2次元有限要素法

正方領域$(0,1)\times (0,1)$に定義されたPoisson方程式の境界値問題を考えます。

$$-\Delta u = f \mbox{ in } \Omega, \quad u=0 \mbox{ on }\partial \Omega$$

上記の式の両側にテスト関数$v\in H^1_0(\Omega)$をかけて、積分すると、以下の変分式が得られます。

$$ (\nabla u, \nabla v)_{\Omega} = (f,v)_{\Omega},\quad \forall v \in H^1_0(\Omega) $$

ここで、 $$ (f, g)_{\Omega}:=\int_{\Omega} f\cdot g d\Omega \:. $$

有限要素法は領域のメッシュ分割を用いて立てられる有限次元の関数空間$V^h$(有限要素空間という)で上記の変分式を解いています。 即ち、以下の式を満たす$u_h\in V^h$を求める。

$$ (\nabla u_h, \nabla v_h)_{\Omega} = (f,v_h)_{\Omega},\quad \forall v_h \in V^h $$

関数空間$V^h$の作り方については、基本的に区分的な多項式が成すベクトル空間です。メッシュ分割の形、多項式の次数と関数の連続性の設定によって、多くのバリエーションがあります。例えば、以下の特徴を持っている「一次適合有限要素関数空間」はよく使用されています。

  • 領域を三角メッシュに分割する
  • 三角メッシュの上に定義される区分的な1次多項式であり、全体としては連続である。
  • 境界で値が0となる。

上記の関数空間$V^h$は$H_0^1(\Omega)$の部分空間となっているので、適合有限要素関数空間(conforming finite element space)と呼ばれます。 $V^h\not\subset H_0^1(\Omega)$の場合、非適合有限要素空間(non-conforming finite element space )と呼ばれます。

以下のコードでは、一次適合有限要素空間を利用してPoisson方程式の近似解を解いています。

実際の計算には、$f=2x(1-x)+2y(1-y)$を使っています。この場合、厳密解は $$u=x(1−x)y(1−y)$$ です。

劉 雪峰
最終更新:2019年12月10日
最初作成:2015年6月4日

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

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

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

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

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

In [2]:
N=8
mesh = UnitSquareMesh(N, N)
In [3]:
%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)
Out[3]:
[<matplotlib.lines.Line2D at 0x7f78e8b8be90>,
 <matplotlib.lines.Line2D at 0x7f78e8b8bf50>]

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

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

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

In [7]:
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)_{\Omega} = (f,v)_{\Omega} \quad \forall v \in V$$

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

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

a = inner(grad(u), grad(v))*dx
b = f*v*dx
--- Instant: compiling ---

4. 有限要素解の計算

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

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

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

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

In [10]:
%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[10]:
<matplotlib.collections.TriMesh at 0x7f78ed5e3890>
In [12]:
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[12]:
<mpl_toolkits.mplot3d.art3d.Poly3DCollection at 0x7f78e7c9bdd0>

6. 有限要素解の誤差

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

$$ \|Err\|_{H^1} = \sqrt{ \int_\Omega \nabla Err \cdot \nabla Err~ dx} \sqrt{ \int_\Omega |\nabla Err|^2~ dx} $$
In [14]:
exact_solution = Expression("x[0]*(1- x[0])*x[1]*(1-x[1])",degree=6)
exact_solution_grad = Expression(["(1-2*x[0])*x[1]*(1-x[1])","(1-2*x[1])*x[0]*(1-x[0])"],degree=6)

fem_error = fem_u - exact_solution
fem_error_grad = grad(fem_u) - exact_solution_grad

l2_error = sqrt(assemble(fem_error*fem_error*dx))
print "Error in L2 norm:      %f"% l2_error

h1_error = sqrt(assemble( inner( fem_error_grad , fem_error_grad )*dx))
print "Error in H1-semi norm: %f"% h1_error

# print "Error in H1-semi norm: %f"% norm(fem_error,"H1",mesh)
--- Instant: compiling ---
Error in L2 norm:      0.001441
Error in H1-semi norm: 0.030161

上記のコードでは、厳密解の勾配をあらかじめ用意するのは必要です。FEniCSで厳密解の勾配を計算するために、厳密解を有限要素空間への射影または補間を利用して数値的に計算することができます。

In [28]:
#The interpolation is used.
V4 = FunctionSpace(mesh, "Lagrange", 4);
exact_solution = Expression("x[0]*(1- x[0])*x[1]*(1-x[1])",degree=4)
interp_exact_solution = interpolate(exact_solution, V4);
interp_fem_u = interpolate(fem_u, V4);
fem_error = interp_fem_u - interp_exact_solution

l2_error = sqrt(assemble(fem_error*fem_error*dx))
h1_error = sqrt(assemble( inner( fem_error_grad , fem_error_grad )*dx))

print("Error estimated by function interpolation.")
print("Error in L2 norm (by interpolation):      %f"% l2_error)
print("Error in H1-semi norm (by interpolation): %f"% h1_error)

print("\n---\n")

fem_error = project(fem_u - exact_solution, V4);
print("Error estimated by function projection.")
print( "Error in L2 norm (by projection):     %f"% norm(fem_error,"L2",mesh))
print( "Error in H1-semi norm  (by projection):  %f"% norm(fem_error,"H1",mesh))
Error estimated by function interpolation.
Error in L2 norm (by interpolation):      0.001441
Error in H1-semi norm (by interpolation): 0.030161

---

Error estimated by function projection.
Error in L2 norm (by projection):     0.001441
Error in H1-semi norm  (by projection):  0.030196

7. 計算コードのまとめ

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

#------------ mesh triangulation and FEM space --------------#
N=16
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]) ")
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])")
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")
--- Instant: compiling ---

TypeErrorTraceback (most recent call last)
<ipython-input-29-1b03ca79396c> in <module>()
     18 u = TrialFunction(V)
     19 v = TestFunction(V)
---> 20 f = Expression("2*x[0]*(1- x[0]) + 2*x[1]*(1- x[1]) ")
     21 a = inner(grad(u), grad(v))*dx
     22 b = f*v*dx

/usr/lib/python2.7/dist-packages/dolfin/functions/expression.pyc in __init__(self, cppcode, *args, **kwargs)
    155         # Select an appropriate element if not specified
    156         if element is None:
--> 157             element = _auto_select_element_from_shape(value_shape, degree, cell)
    158         else:
    159             # Check that we have an element

/usr/lib/python2.7/dist-packages/dolfin/functions/expression.pyc in _auto_select_element_from_shape(shape, degree, cell)
    886     "Automatically select an appropriate element from cppcode."
    887     if not isinstance(degree, int):
--> 888         raise TypeError("Expression needs an integer argument 'degree' if no 'element' is provided.")
    889 
    890     # Default element, change to quadrature when working

TypeError: Expression needs an integer argument 'degree' if no 'element' is provided.
In [ ]: