２次元混合有限要素法

正方領域$(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/ )

from fenics import *

2. メッシュ分割と有限要素空間の設定

まず、メッシュ分割を行います。

N=8
mesh = UnitSquareMesh(N, N)

１次混合有限要素法の関数Raviart-Thomas空間$RT$と区分的な定数区間$DG$を作ります。空間$V$は$RT$と$DG$の直和空間です。

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$$

以下のコードは変分式の両辺の内積を計算します。

(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に格納されます。

fem_sol = Function(V)
solve(a == b, fem_sol)

(p,u) = fem_sol.split()
solution_u = u

Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.

5. 有限要素解の描画：カーラーマップと３Dの図

有限要素解の色マップを描画します。分割の各点での値を再計算することが必要です。fem_uのベクトルの自由度の順番とmeshの節点リストの順番は互いに独立であることを注意しなさい。

%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')

/usr/lib/python2.7/dist-packages/matplotlib/font_manager.py:273: UserWarning: Matplotlib is building the font cache using fc-list. This may take a moment.
  warnings.warn('Matplotlib is building the font cache using fc-list. This may take a moment.')

<matplotlib.collections.TriMesh at 0x7f6509827c90>

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)

<mpl_toolkits.mplot3d.art3d.Poly3DCollection at 0x7f650984c750>

6. 有限要素解の誤差

以下は$(u-u_h)$の$L^2$ノルムと$(p-\nabla u)$の$L^2$ノルムを計算します。

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)

L2 norm of (u_h-u):  0.003936
L2 norm of (p-grad(u)): 0.018379
Another way to calculate the norm
L2-norm of (p-grad(u)): 0.018379

7. 計算コードのまとめ

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)

L2 norm of (u_h-u):      0.000313
L2 norm of (p-grad(u)):      0.001553