Лабораторная работа №3¶
ДВ Метод конечных элементов, ММФ, БГУ
преподаватель Лаврова О.А., октябрь 2017
Задача Дирихле для уравнения Пуассона¶
Постановка задачи¶
Для заданного точного решения $u(x,y)$ построить задачу Дирихле для уравнения Пуассона в единичном квадрате $[0,1] \times [0,1]$
\begin{equation*} \left\{ \begin{array}{ll} -\Delta u=f, & \quad (x,y) \in \Omega=(0,1)^2\\ u=g, & \quad (x,y) \in \partial \Omega. \end{array} \right. \end{equation*}Решить построенную задачу средствами библиотеки FEniCS и сравнить визуально численное решение, полученное методом конечных элементов, с точным решением $u(x,y)$.
Например, полагаем $u(x,y) = x^3+y^3$. Тогда задача Дирихле примет вид
\begin{equation*} \left\{ \begin{array}{ll} -\Delta u=-6 x - 6 y, & \quad (x,y) \in \Omega=(0,1)^2\\ u(x,y)=y^3, & \quad (x,y) \in \partial \Omega \cap \{(x,y) : x=0\} \\ u(x,y)=1+y^3, & \quad (x,y) \in \partial \Omega \cap \{(x,y) : x=1\}\\ u(x,y)=x^3, & \quad (x,y) \in \partial \Omega \cap \{(x,y) : y=0\}\\ u(x,y)=x^3+1. & \quad (x,y) \in \partial \Omega \cap \{(x,y) : y=1\}\\ \end{array} \right. \end{equation*}Реализация¶
Будет осуществлена поэтапно:
- Построение сетки
- Определение конечно-элементного пространства
- Определение граничных условий
- Определение вариационной задачи
- Решение системы линейных алгебраических уравнений
- Постпроцессинг (визуализация, вычисление ошибки)
Подключаем библиотеки¶
from fenics import *
import fenics
import matplotlib.pyplot as plt
%matplotlib inline
import numpy as np
Получение справки по функциям библиотеки fenics¶
help(fenics.plot)
или с помощью pydoc
!pydoc fenics.mesh
1 Этап. Построение сетки¶
mesh = UnitSquareMesh(20, 8)
plot(mesh)
2 Этап. Определение конечно-элементного пространства¶
V1 = FunctionSpace(mesh, 'P', 1) # линейные элементы на треугольниках
V2 = FunctionSpace(mesh, 'P', 2) # квадратичные элементы на треугольниках
V3 = FunctionSpace(mesh, 'P', 3) # кубические элементы на треугольниках
Важно: In FEniCS we DO NOT specify the boundary conditions as part of the function space, so it is sufficient to work with one common space $V$ for both the trial and test functions in the program [From FEniCS Tutorial 2016]
u1 = TrialFunction(V1)
v1 = TestFunction(V1)
u2 = TrialFunction(V2)
v2 = TestFunction(V2)
2.1. Определение точного решения исходной задачи¶
Точное решение будет использовано на этапе постпроцессинга
u_exact = Expression('pow(x[0],3) + pow(x[1],3)', degree=2)
Для визуализации необходимо описание функции (тип Expression) преобразовать в функцию
u_exact_plt = interpolate(u_exact,V3)
Важно: if an Expression is used to represent an exact solution which is used to evaluate the accuracy of a computed solution, a higher degree than for the space $V$ must be used for the expression (one or two degree higher) [From FEniCS Tutorial 2016]
3 Этап. Определение граничных условий¶
3.1. Описание границы¶
tol = 1E-14
def boundary(x):
return abs(x[0]) < tol or abs(x[1]) < tol \
or abs(x[0] - 1) < tol or abs(x[1] - 1) < tol
3.2. Задание функции из условия Дирихле¶
g_D1 = Expression('pow(x[0],3) + pow(x[1],3)', degree=1)
g_D2 = Expression('pow(x[0],3) + pow(x[1],3)', degree=2)
Важно: To obtain optimal (order of) accuracy in computations, it is usually a good choice to use the same degree as for the space $V$ that is used for the trial and test functions [From FEniCS Tutorial 2016]
3.3. Граничное условие¶
bc1 = DirichletBC(V1, g_D1, boundary)
bc2 = DirichletBC(V2, g_D2, boundary)
4 Этап. Определение вариационной задачи¶
Определение функции правой части¶
f = Expression('-6*x[0]-6*x[1]', degree=1)
a1 = dot(grad(u1), grad(v1))*dx # линейные элементы
L1 = f*v1*dx
a2 = dot(grad(u2), grad(v2))*dx # квадратичные элементы
L2 = f*v2*dx
Важно: The language used to express weak forms is called UFL (Unified Form Language) and is an integral part of FEniCS [From FEniCS Tutorial 2016]
5 Этап. Решение системы линейных алгебраических уравнений¶
u1 = Function(V1) # переопределение функции u1
solve(a1 == L1, u1, bc1)
u2 = Function(V2) # переопределение функции u2
solve(a2 == L2, u2, bc2)
6 Этап. Постпроцессинг¶
6.1. Визуализация¶
plot(u1, interactive = True)
#colorbar()
plot(u1-u_exact_plt)
interactive()
Важно: The call to interactive() is usually placed at the end of a program that creates plot [From FEniCS Tutorial 2016]
Важно: The built-in plotting in FEniCS may not work as expected when running inside a FEniCS Docker container. For Docker users who need plotting, it is recommended to work within a Jupyter/FEniCS notebook (command fenicsproject notebook) [From FEniCS Tutorial 2016]
#from IPython.display import HTML
#HTML(X3DOM().html(u1))
6.2. Вычисление ошибки¶
1) Вычисление ошибки в норме $L_2$¶
error_L2_1 = errornorm(u_exact,u1,'L2')
error_L2_2 = errornorm(u_exact,u2,'L2')
print(error_L2_1,error_L2_2)
Замечание: Обратите внимание на уменьшение ошибки за счет перехода из пространства $P1$-элементов в пространство $P2$-элементов на одинаковой сетке.
2) Вычисление ошибки в максимальной норме¶
Данный код реализует вычисление ошибки в узлах сетки (не в степенях свободы конечно-элементого решения)
vertex_values_u_exact = u_exact.compute_vertex_values(mesh)
vertex_values_u1 = u1.compute_vertex_values(mesh)
vertex_values_u2 = u2.compute_vertex_values(mesh)
error_max_1 = np.max(np.abs(vertex_values_u_exact - vertex_values_u1))
error_max_2 = np.max(np.abs(vertex_values_u_exact - vertex_values_u2))
print(error_max_1,error_max_2)
Важно: We have here used the maximum and absolute value functions from numpy, because these are much more efficient for large arrays (a factor of 30) than Python's built-in max and abs functions [From FEniCS Tutorial 2016]