2次元有限要素法 (簡易版)

正方領域$(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)$$ です。

劉 雪峰
最終更新:2017年9月4日
最初作成:2015年6月4日

注意

該当コードはFEniCS 2016.01以降のバージョンで実行できます。2016.01より古いFEniCSで実行する場合、不具合のことがあります。

1. 準備(パッケージの導入)

DOLFINのパッケージを導入します。

  • (DOLFINはFEniCSプロジェクトの核部分であり、数行のコードで有限要素法計算ができます。
  • FEniCS Project: https://fenicsproject.org/ ) ``
In [1]:
from fenics import *

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

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

In [2]:
N=8
mesh = UnitSquareMesh(N, N)

一次適合有限要素法の関数空間$V$を作ります。空間$V$は区分一次連続関数によって構成される線形関数空間です。

In [3]:
V = FunctionSpace(mesh, "Lagrange", 1)

関数boundaryで指定される境界でDirichlet境界値条件の設定を行います。

In [4]:
def boundary(x):
    return near(x[0], 0) or near(x[0],1.0) or near(x[1], 0) or near(x[1],1.0)
u0 = Constant(0.0)
bc = DirichletBC(V, u0, boundary)

3. 変分式を定義する

有限要素法は関数空間$V$の中で以下の変分式を解きます。

$$(\nabla u, \nabla v) = (f,v) \quad \forall v \in V$$

That is,

$ \int_\Omega \nabla u \cdot \nabla v dxdy $

ただし、$f=2x(1-x)+2y(1-y)$。

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

In [9]:
u = TrialFunction(V)
v = TestFunction(V)
f = Expression("2*x[0]*(1- x[0]) + 2*x[1]*(1- x[1])", degree=4)  # C 言語で書かれた数式

a = inner(grad(u), grad(v))*dx
b = f*v*dx

4. 有限要素解の計算

変分式はa=bを解いて、有限要素解が得られます。計算の結果はfem_uに格納されます。

In [10]:
fem_u = Function(V)
solve(a == b, fem_u, bc)

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

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

In [12]:
%matplotlib inline
import matplotlib.pyplot as plt
import matplotlib.tri as tri



nodes = mesh.coordinates()
solution_vec = fem_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[12]:
<matplotlib.collections.TriMesh at 0x7fcfa16f8150>
In [13]:
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[13]:
<mpl_toolkits.mplot3d.art3d.Poly3DCollection at 0x7fcfa15856d0>

6. 有限要素解の誤差

以下は$(u-u_h)$の$L_2$ノルムと$H^1_0$ノルムを計算します。 $(u-u_h)$のノルムを直接に計算できないので、誤差式の4次Lagrange有限要素空間への射影を使って、ノルムを計算します。

In [16]:
V4 = FunctionSpace(mesh, "CG", 4);
exact_solution = Expression("x[0]*(1- x[0])*x[1]*(1-x[1])",degree=4)
fem_error = project(fem_u - exact_solution, V4);
print "Error in L2 norm:      %f"% norm(fem_error,"L2",mesh=mesh)
print "Error in H1-semi norm: %f"% norm(fem_error,"H1",mesh=mesh)
Error in L2 norm:      0.001441
Error in H1-semi norm: 0.030196

7. 計算コードのまとめ

In [1]:
from dolfin import *
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import matplotlib.tri as tri
%matplotlib inline

#------------ mesh triangulation and FEM space --------------#
N=64
mesh = UnitSquareMesh(N, N)
V = FunctionSpace(mesh, "Lagrange", 1)

#---------------   boundary condition  ----------------------------#
def boundary(x):
    return near(x[0], 0) or near(x[0],1.0) or near(x[1], 0) or near(x[1],1.0)

u0 = Constant(0.0)
bc = DirichletBC(V, u0, boundary)

u = TrialFunction(V)
v = TestFunction(V)
f = Expression("2*x[0]*(1- x[0]) + 2*x[1]*(1- x[1]) ", degree=4)
a = inner(grad(u), grad(v))*dx
b = f*v*dx

#---------------   Compute solution   ------------------------------#
fem_u = Function(V)
solve(a == b, fem_u, bc)

#---------------- Draw solution   -------------------------------------#
nodes = mesh.coordinates()
solution_vec = fem_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 ------------------------------------#
exact_solution = Expression("x[0]*(1- x[0])*x[1]*(1-x[1])", degree=4)
mesh = refine(mesh)
V4 = FunctionSpace(mesh, "DG", 4);
fem_error = project(fem_u - exact_solution, V4);
print "Error in L2 norm: %f"% norm(fem_error,"L2")
print "Error in H1 norm: %f"% norm(fem_error,"H1")
/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.')
Error in L2 norm: 0.000023
Error in H1 norm: 0.003803
In [ ]: