２次元有限要素法 (詳細版)¶
正方領域$(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/ )
In [1]:
from fenics import *
1. メッシュ分割¶
まず、粗いメッシュ分割を行います。
In [2]:
N=3
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]:
In [4]:
#----- メッシュの確認 ------------
print "要素のリストの数"
print len( mesh.cells() )
print "要素のリスト"
print mesh.cells()
print "分割の節点のリスト"
print mesh.coordinates()
2. 有限要素空間と変分式¶
一次適合有限要素法の関数空間$V$を作ります。空間$V$は区分一次連続関数によって構成される線形関数空間です。
In [5]:
#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境界値条件の設定を行います。
In [6]:
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)$。
以下のコードは変分式の両辺の内積を計算します。
In [7]:
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に格納されます。
In [8]:
# Compute solution
fem_u = Function(V)
solve(a == b, fem_u, bc)
有限要素解の色マップを描画します。分割の各点での値を再計算することが必要です。fem_uのベクトルの自由度の順番とmeshの節点リストの順番は違いますので、注意が必要です。
In [9]:
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
有限要素解を直接に求める。¶
剛性行列(stiffness matrix)$A$とロードベクトル(load vector)の値を確認して、連立方程式を解いてから解を求めることもできます。
In [11]:
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
In [12]:
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')
Out[12]:
3Dの図を描画します。
In [14]:
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)
Out[14]:
有限要素解の誤差を計算します。
In [16]:
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
計算コードのまとめ¶
In [2]:
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")
