２次元混合有限要素法¶
正方領域$(0,1)\times (0,1)$に定義されたPoisson方程式の境界値問題を考えます。
$$-\Delta u = f \mbox{ in } \Omega, \quad u=0 \mbox{ on }\partial \Omega$$
実際の計算には、$f=2x(1-x)+2y(1-y)$を使っています。この場合、厳密解は $$u=x(1−x)y(1−y)$$ です。
劉 雪峰
最終更新：２０１７年９月４日
最初作成：２０１５年６月４日
最終更新：２０１７年９月４日
最初作成：２０１５年６月４日
1. 準備（パッケージの導入）¶
DOLFINのパッケージを導入します。
- (DOLFINはFEniCSプロジェクトの核部分であり、数行のコードで有限要素法計算ができます。
- FEniCS Project: https://fenicsproject.org/ )
In [1]:
from fenics import *
2. メッシュ分割と有限要素空間の設定¶
まず、メッシュ分割を行います。
In [2]:
N=8
mesh = UnitSquareMesh(N, N)
１次混合有限要素法の関数Raviart-Thomas空間$RT$と区分的な定数区間$DG$を作ります。空間$V$は$RT$と$DG$の直和空間です。
In [31]:
degree=1
RT = FiniteElement("RT", mesh.ufl_cell(), degree)
DG = FiniteElement("DG", mesh.ufl_cell(), degree-1)
V = FunctionSpace( mesh , MixedElement ([RT, DG]) )
混合法をPoisson問題の境界値問題を考える時、Dirichlet境界値条件は自然条件となって、境界における解の値の設定は不要です。
3. 変分式を定義する¶
有限要素法は関数空間$V$の中で以下の変分式を解きます。
Find $(p,u)$ in $RT \times DG$, s.t.,
$$(p,q) + (div\: q,u) =0 \quad \forall q \in RT$$
$$(div\: p + f ,v) =0 \quad \forall v \in DG$$
ただし、$f=2x(1-x)+2y(1-y)$。
上記の２つの変分式を一つにまとめてます。
Find $(p,u)$ in $RT \times DG$, s.t., $$(p,q) + (div\: q,u) + (div\: p, v) =(- f ,v) \quad \forall (q,v) \in RT \times DG$$
以下のコードは変分式の両辺の内積を計算します。
In [11]:
(p,u) = TrialFunctions(V)
(q,v) = TestFunctions(V)
f = Expression("2*x[0]*(1- x[0]) + 2*x[1]*(1- x[1])", degree=4) # C 言語で書かれた数式
a = (inner(p,q)+ div(p)*v + div(q)*u)*dx
b = -f*v*dx
4. 有限要素解の計算¶
変分式はa=bを解いて、有限要素解が得られます。計算の結果はfem_uに格納されます。
In [12]:
fem_sol = Function(V)
solve(a == b, fem_sol)
(p,u) = fem_sol.split()
solution_u = u
5. 有限要素解の描画：カーラーマップと３Dの図¶
有限要素解の色マップを描画します。分割の各点での値を再計算することが必要です。fem_uのベクトルの自由度の順番とmeshの節点リストの順番は互いに独立であることを注意しなさい。
In [13]:
%matplotlib inline
import matplotlib.pyplot as plt
import matplotlib.tri as tri
nodes = mesh.coordinates()
solution_vec = u.compute_vertex_values(mesh)
plt.gca().set_aspect('equal')
my_triangle = tri.Triangulation(nodes[:, 0], nodes[:, 1], mesh.cells())
plt.tripcolor(my_triangle, solution_vec, shading='gouraud')
Out[13]:
In [14]:
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot_trisurf(nodes[:,0],nodes[:,1],solution_vec, cmap=plt.cm.Spectral, linewidth=0.1, antialiased=True)
Out[14]:
6. 有限要素解の誤差¶
以下は$(u-u_h)$の$L^2$ノルムと$(p-\nabla u)$の$L^2$ノルムを計算します。
In [28]:
V4 = FunctionSpace(mesh, "CG", 4);
exact_solution = Expression("x[0]*(1- x[0])*x[1]*(1-x[1])", degree=4, domain=mesh)
fem_error = project(solution_u - exact_solution, V4);
print "L2 norm of (u_h-u): %f"% norm(fem_error,"L2",mesh=mesh)
DG2=VectorFunctionSpace(mesh, "DG",3, dim=2)
fem_grad_error = project(p - grad(exact_solution), DG2);
print "L2 norm of (p-grad(u)): %f"% norm(fem_grad_error,"L2",mesh=mesh)
#Another way to calculate the norm
print "Another way to calculate the norm"
u=TrialFunction(DG2)
v=TestFunction(DG2)
a2 = inner(u,v)*dx
b2 = inner( p - grad(exact_solution), v)*dx
L2_projection = Function(DG2)
solve(a2 == b2, L2_projection)
print "L2-norm of (p-grad(u)): %f"% norm(L2_projection,"L2",mesh=mesh)
7. 計算コードのまとめ¶
In [33]:
from dolfin import *
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
%matplotlib inline
#------------ mesh triangulation and FEM space --------------#
N=16*6
mesh = UnitSquareMesh(N, N)
degree=1
RT = FiniteElement("RT", mesh.ufl_cell(), degree)
DG = FiniteElement("DG", mesh.ufl_cell(), degree-1)
V = FunctionSpace( mesh , MixedElement ([RT, DG]) )
#--------------- Variational equation ----------------------------#
(p,u) = TrialFunctions(V)
(q,v) = TestFunctions(V)
f = Expression("2*x[0]*(1- x[0]) + 2*x[1]*(1- x[1])", degree=4, domain=mesh) # C 言語で書かれた数式
a = (inner(p, q)+ div(p)*v + div(q)*u) *dx
b = -f*v*dx
#--------------- Compute solution ------------------------------#
fem_sol = Function(V)
solve(a == b, fem_sol)
(p,u) = fem_sol.split()
solution_u = u
#---------------- Draw solution -------------------------------------#
nodes = mesh.coordinates()
solution_vec = solution_u.compute_vertex_values(mesh)
fig = plt.figure(1)
ax = fig.gca(projection='3d')
ax.plot_trisurf(nodes[:,0],nodes[:,1],solution_vec, cmap=plt.cm.Spectral, linewidth=0.1, antialiased=True)
fig = plt.figure(2)
plt.gca().set_aspect('equal')
my_triangle = tri.Triangulation(nodes[:, 0], nodes[:, 1], mesh.cells())
plt.tripcolor(my_triangle, solution_vec, shading='gouraud')
#------------- Error computation ------------------------------------#
V4 = FunctionSpace(mesh, "CG", 4);
exact_solution = Expression("x[0]*(1- x[0])*x[1]*(1-x[1])", degree=4, domain=mesh)
fem_error = project(solution_u - exact_solution, V4);
print "L2 norm of (u_h-u): %f"% norm(fem_error,"L2",mesh=mesh)
DG2=VectorFunctionSpace(mesh, "DG",3,2)
fem_grad_error = project(p - grad(exact_solution), DG2);
print "L2 norm of (p-grad(u)): %f"% norm(fem_grad_error,"L2",mesh=mesh)
