２次元有限要素法 (詳細版)

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

劉　雪峰
最終更新：２０１７年９月４日
最初作成：２０１５年６月４日

注意

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

準備（パッケージの導入）

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

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

from fenics import *

1. メッシュ分割

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

N=3
mesh = UnitSquareMesh(N, N)

メッシュの確認を確認します。

%matplotlib inline

import matplotlib.pyplot as plt
import matplotlib.tri as tri

nodes = mesh.coordinates();
elements = tri.Triangulation(nodes[:, 0], nodes[:, 1], mesh.cells());
plt.triplot(elements)

/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.lines.Line2D at 0x7f76fa0bac90>,
 <matplotlib.lines.Line2D at 0x7f76fa0bae50>]

#-----  メッシュの確認 ------------
print "要素のリストの数"
print len( mesh.cells() )
print "要素のリスト"
print mesh.cells()
print "分割の節点のリスト"
print mesh.coordinates()

要素のリストの数
18
要素のリスト
[[ 0  1  5]
 [ 0  4  5]
 [ 1  2  6]
 [ 1  5  6]
 [ 2  3  7]
 [ 2  6  7]
 [ 4  5  9]
 [ 4  8  9]
 [ 5  6 10]
 [ 5  9 10]
 [ 6  7 11]
 [ 6 10 11]
 [ 8  9 13]
 [ 8 12 13]
 [ 9 10 14]
 [ 9 13 14]
 [10 11 15]
 [10 14 15]]
分割の節点のリスト
[[ 0.          0.        ]
 [ 0.33333333  0.        ]
 [ 0.66666667  0.        ]
 [ 1.          0.        ]
 [ 0.          0.33333333]
 [ 0.33333333  0.33333333]
 [ 0.66666667  0.33333333]
 [ 1.          0.33333333]
 [ 0.          0.66666667]
 [ 0.33333333  0.66666667]
 [ 0.66666667  0.66666667]
 [ 1.          0.66666667]
 [ 0.          1.        ]
 [ 0.33333333  1.        ]
 [ 0.66666667  1.        ]
 [ 1.          1.        ]]

2. 有限要素空間と変分式

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

#import logging
#ffc_logger = logging.getLogger('FFC')
#ffc_logger.setLevel(logging.ERROR)

V = FunctionSpace(mesh, "Lagrange", 1)

Remark: To disable the debug output:

import logging
ffc_logger = logging.getLogger('FFC')
print 'Log level was:', logging.getLevelName(ffc_logger.getEffectiveLevel())
ffc_logger.setLevel(logging.ERROR)
print 'Log level is now:', logging.getLevelName(ffc_logger.getEffectiveLevel())

関数boundaryによって境界部分を定義して、Dirichlet境界値条件の設定を行います。

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)

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

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

即ち、

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

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

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

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

3. 有限要素解を求める

a,bを使って、有限要素解を計算できます。計算の結果はfem_uに格納されます。

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

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

nodes = mesh.coordinates()
solution_vec = fem_u.compute_vertex_values(mesh)
## Or, follow the way below.
#solution_vec = fem_u.vector()[vertex_to_dof_map(V)]
print solution_vec

[ 0.          0.          0.          0.          0.          0.04526749
  0.04526749  0.          0.          0.04526749  0.04526749  0.          0.
  0.          0.          0.        ]

有限要素解を直接に求める。

剛性行列(stiffness matrix)$A$とロードベクトル(load vector)の値を確認して、連立方程式を解いてから解を求めることもできます。

import numpy as np

s_A,s_b=assemble_system(a,b,bc)

MatA = as_backend_type(s_A)
Vecb = as_backend_type(s_b)


print "剛性行列"
#The DOF order in matrix is different from the one in mesh.
dof_map = vertex_to_dof_map(V) 
mat_A = MatA.array()[np.ix_(dof_map,dof_map)]
print mat_A

print "ロードベクトル"
vec_b = Vecb.array()[np.ix_(dof_map)]
print vec_b

print "連立一次方程式を解いて、解を表示します。"
solution = np.linalg.solve(mat_A, vec_b)
print solution

剛性行列
[[ 2.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
 [ 0.  3.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
 [ 0.  0.  3.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
 [ 0.  0.  0.  1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  3.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  4. -1.  0.  0. -1.  0.  0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0. -1.  4.  0.  0.  0. -1.  0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.  0.  3.  0.  0.  0.  0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.  0.  0.  3.  0.  0.  0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0. -1.  0.  0.  0.  4. -1.  0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  0. -1.  0.  0. -1.  4.  0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  3.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  3.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  3.  0.]
 [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  2.]]
ロードベクトル
[ 0.          0.          0.          0.          0.          0.09053498
  0.09053498  0.          0.          0.09053498  0.09053498  0.          0.
  0.          0.          0.        ]
連立一次方程式を解いて、解を表示します。
[ 0.          0.          0.          0.          0.          0.04526749
  0.04526749  0.          0.          0.04526749  0.04526749  0.          0.
  0.          0.          0.        ]

import matplotlib.pyplot as plt
import matplotlib.tri as tri

plt.gca().set_aspect('equal')
my_triangle = tri.Triangulation(nodes[:, 0], nodes[:, 1], mesh.cells())
plt.tripcolor(my_triangle, solution_vec, shading='gouraud')

<matplotlib.collections.TriMesh at 0x7f76f9ff20d0>

3Dの図を描画します。

from mpl_toolkits.mplot3d import Axes3D
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)

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

有限要素解の誤差を計算します。

N=8
mesh2 = UnitSquareMesh(N, N)

V4 = FunctionSpace(mesh2, "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=mesh2)
print "Error in H2-semi norm: %f"% norm(fem_error,"H1",mesh=mesh2)

fem_error = fem_u - exact_solution
L2_error = sqrt(assemble(fem_error*fem_error*dx ))
#H1_error = assemble( inner(nabla_grad(fem_error), nabla_grad(fem_error))*dx(mesh) )
#The above way to evalute H1_error does not work.

print "Error in L2 norm (by integral directly): %f"% L2_error

##print H1_error

Error in L2 norm:      0.009658
Error in H2-semi norm: 0.068945
Error in L2 norm (by integral directly): 0.009153

計算コードのまとめ

from dolfin import *
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import matplotlib.tri as tri

%matplotlib inline


N=8*8
mesh = UnitSquareMesh(N, N)
V = FunctionSpace(mesh, "Lagrange", 1)

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")

Error in L2 norm: 0.000023
Error in H1 norm: 0.003803