はじめに
この記事では2次元のポアソン方程式を有限要素法で解くためのシンプルなPythonプログラムを紹介します。
前記事では本記事で出てくる操作の意味について解説しました。
実装編では、その理論を実際のコードへ落とし込みます。局所剛性行列の計算、大域剛性行列の構築、境界条件の適用といった処理を順に実装し、最終的に数値解を求めていきます。
前回記事
https://qiita.com/fried_nasubi/items/5f30d583a2f6345b918d
局所剛性行列の計算
前記事では、有限要素法の実装において各三角形ごとに局所剛性行列を計算し、それらを足し合わせることで大域剛性行列を構成することを説明しました。
ここでは、三角形の頂点座標から局所剛性行列を計算する方法を導出します。
基底関数の勾配
基底関数は
\phi_{i1}
=
\frac{1}{2A_i}
\left\{
b_{i1}x
+
c_{i1}y
+
(x_{i2}y_{i3}-x_{i3}y_{i2})
\right\}
のように表されるのでした。
したがって、その勾配は
\nabla\phi_{ij}
=
\frac{1}{2A_i}
\begin{pmatrix}
b_{ij}\\
c_{ij}
\end{pmatrix}
です。
この結果、2つの基底関数の勾配の内積は
(\nabla\phi_{ia})
\cdot
(\nabla\phi_{ib})
=
\frac{1}{4A_i^2}
\left(
b_{ia}b_{ib}
+
c_{ia}c_{ib}
\right)
となります。
要素積分
局所剛性行列の成分は
K_{ab}^{(i)}
=
\int_{\Omega_i}
dS\,
(\nabla\phi_{ia})
\cdot
(\nabla\phi_{ib})
で定義されます。
先ほど求めた式を代入すると、
K_{ab}^{(i)}
=
\frac{1}{4A_i^2}
\left(
b_{ia}b_{ib}
+
c_{ia}c_{ib}
\right)
\int_{\Omega_i}dS
となります。
ここで、
\int_{\Omega_i}dS=A_i
であるため、
K_{ab}^{(i)}
=
\frac{1}{4A_i}
\left(
b_{ia}b_{ib}
+
c_{ia}c_{ib}
\right)
が得られます。
局所剛性行列
以上より、第 $i$ 要素の局所剛性行列は
K^{(i)}
=
\frac{1}{4A_i}
\begin{pmatrix}
b_{i1}b_{i1}+c_{i1}c_{i1}
&
b_{i1}b_{i2}+c_{i1}c_{i2}
&
b_{i1}b_{i3}+c_{i1}c_{i3}
\\
b_{i2}b_{i1}+c_{i2}c_{i1}
&
b_{i2}b_{i2}+c_{i2}c_{i2}
&
b_{i2}b_{i3}+c_{i2}c_{i3}
\\
b_{i3}b_{i1}+c_{i3}c_{i1}
&
b_{i3}b_{i2}+c_{i3}c_{i2}
&
b_{i3}b_{i3}+c_{i3}c_{i3}
\end{pmatrix}
となります。
この式をそのまま実装することも可能ですが、行列演算を用いることでより簡潔に記述できます。
行列による表現
まず、三角形の形状を表す行列
J_i
=
\begin{pmatrix}
x_{i2}-x_{i1}
&
y_{i2}-y_{i1}
\\
x_{i3}-x_{i1}
&
y_{i3}-y_{i1}
\end{pmatrix}
を定義します。
これは三角形ごとに定義されるヤコビ行列です。
また、
B
=
\begin{pmatrix}
-1 & 1\\
0 & -1\\
1 & 0
\end{pmatrix}
を定義します。
さらに、
b_{i1}=y_{i2}-y_{i3},
\quad
b_{i2}=y_{i3}-y_{i1},
\quad
b_{i3}=y_{i1}-y_{i2}
c_{i1}=x_{i3}-x_{i2},
\quad
c_{i2}=x_{i1}-x_{i3},
\quad
c_{i3}=x_{i2}-x_{i1}
を用いると、
J_i
=
\begin{pmatrix}
c_{i3} & -b_{i3}\\
-c_{i2} & b_{i2}
\end{pmatrix}
と書くことができます。
ここで、
BJ_i
を計算すると、
BJ_i
=-
\begin{pmatrix}
b_{i1} & c_{i1}\\
b_{i2} & c_{i2}\\
b_{i3} & c_{i3}
\end{pmatrix}
となります。
すなわち、各行には基底関数の勾配を特徴づける係数
(b_{ij},c_{ij})
が並びます。
そのため、
(BJ_i)(BJ_i)^T
の $(a,b)$ 成分は
b_{ia}b_{ib}
+
c_{ia}c_{ib}
となります。
これは先ほど導出した局所剛性行列の分子そのものです。
したがって、
K^{(i)}
=
\frac{1}{4A_i}
(BJ_i)(BJ_i)^T
と書くことができます。
さらに、
2|A_i|
=
|\det J_i|
であるため、
K^{(i)}
=
\frac{1}{2|\det J_i|}
(BJ_i)(BJ_i)^T
が得られます。
実装上の利点
この表現を用いると、三角形の3頂点の座標から
J_i
を構成し、
(BJ_i)(BJ_i)^T
を計算するだけで局所剛性行列が得られます。
そのため、局所剛性行列の各成分を個別に実装する必要がありません。
有限要素法のコードでは、まず各三角形について局所剛性行列を計算し、その結果を対応する節点番号へ加算することで大域剛性行列を構成します。
大域剛性行列の計算
理論編でも説明したように、局所剛性行列から大域剛性行列を作成する際には、
各三角形の頂点がどの節点に対応しているのかが重要になります。
そこで以下のようなリストを作り、これをもとに情報を管理します。
nodes：第$n$成分に第$n$節点の座標を格納したもの。
triangles：三角形の頂点情報を格納したもの。例えば、
triangles[n] = [12, 24, 21]となっていれば、第$n$三角形の第１、２、３頂点はそれぞれ第12、24、21節点であることを示します。
この情報に基づき、
K[np.ix_(nodes_i, nodes_i)] += K_i
と書くことで、
「第$n$三角形の第１頂点は、第12節点に対応する」
↓
「局所剛性行列は、大域剛性行列の第12行、第12列に加算する」
という計算が行われます。
また、大域剛性行列を計算していく行程で同時に荷重ベクトルも計算していきますが、この時点ではノイマン条件に関する計算はしません。あとでまとめて計算します。
ノイマン条件の適用について
ノイマン条件は、弱形式のポアソン方程式の右辺第二項に現れます。
\int_\Gamma dl \, \mathbf{n}\cdot(v\nabla u)
理論編のほうで説明しましたが、ノイマン条件では
\mathbf{n}\cdot \nabla u
つまり、「$\nabla u$、境界線に直交する方向の成分」を直接的に境界値として指定します。のちに扱うディリクレ条件については節点上のみの境界値を決めればいいのですが、ノイマン条件についてはこのように線積分を行うため、節点間の値も補間的に決める必要があります。基底関数を1次関数とする場合、境界値の補間も線形補間とするのが一般的みたいです。
2次元の問題において境界は1次元（曲線）です。そのため、境界上の関数は1変数関数として扱えます。そこで、節点間の境界値を線形補間し、位置を表す変数$\xi$を用いて記述してみましょう。
例えば、第１、２節点上での境界値をそれぞれ$q_1,q_2$とする場合、線形補間は、
\mathbf{n}\cdot \nabla u = \frac{q_2 -q_1}{\xi_2 -\xi_1}(\xi -\xi_0) +q_0
となります。（境界は1次元なので、境界上の位置は単一のパラメータ$\xi$で表現できる。）
大域剛性行列の計算の際には、全ての三角形のすべての頂点に関する基底関数を試験関数として用意し、式を立てています。
ここでまた、計算領域とその分割のイメージ図を掲載します。
ここでは例として、図の第１節点、第２節点、第６節点、第７節点…という境界にノイマン条件を設定している場合を考えましょう。試験関数は$\phi_{1,1}$を選んだ際の計算についてみていきましょう。境界$\Gamma$上で$\phi_{1,1}$が非ゼロとなるのは第１節点と第２節点の間だけなので、
\int_\Gamma dl \, \mathbf{n}\cdot(\phi_{1,1}\nabla u)
= \int_{点１、２間} dl \,\mathbf{n}\cdot(\phi_{1,1}\nabla u)
$\phi_{1,1}$のほうも$\Gamma$上の位置パラメータ$\xi$を用いて表現すると、
\phi_{1,1}(\xi) = \frac{\xi -\xi_1}{\xi_0 -\xi_1}
のようになります。したがって、
\begin{align}
\int_{点１、２間} dl \,\mathbf{n}\cdot(\phi_{1,1}\nabla u) &=
\int_{\xi_0}^{\xi_1} d\xi\,\frac{\xi -\xi_1}{\xi_0 -\xi_1} \left( \frac{q_2 -q_1}{\xi_2 -\xi_1}(\xi -\xi_0) +q_0\right) \\
&= \frac{2q_0 +q_1}{6}(\xi_1 -\xi_0)
\end{align}
です。
試験関数として$\phi_{1,2}$を選んだ際も積分区間は同様です。この場合は、
\phi_{1,2}(\xi) = \frac{\xi -\xi_0}{\xi_1 -\xi_0}
なので、
\int_{点１、２間} dl \,\mathbf{n}\cdot(\phi_{1,2}\nabla u) = \frac{q_0 +2q_1}{6}(\xi_1 -\xi_0)
です。
$\phi_{1,1},\phi_{1,2}$はそれぞれ第１，２節点に対応するため、それぞれの計算結果を荷重ベクトルの第１，２成分に加算します。
同様の作業を次の線分「節点２-節点６」、「節点６-節点７」・・・と続けていき境界の終点まで計算します。第３三角形は頂点のみが境界に接し、辺が境界に接していないため$\phi_{3,2}$に関する計算は存在しません。境界に点だけで接している他の三角形についても同様です。
本記事で紹介する方法では、境界線それぞれにつき適用作業apply_neumannを実施する必要があります。（例えば長方形$ABCD$の辺$AB$、$CD$に対してそれぞれ別に適用作業が必要で、一つにまとめることはできません。）
ディリクレ条件の適用について
ここまでの作業で、行列形式
K \mathbf{u} = \mathbf{f}
はほぼ完成していますが、$u$には既知変数が含まれていて、
u = np.linalg.solve(K, F)
で解くためには既知変数を削除する必要があります。
この既知変数というのはディリクレ条件として直接指定した値です。行、列の削除や移行の考え方については理論編で説明した通りです。
ディリクレ条件の適用の際には大域剛性行列の行、列の削除を行います。プログラム内で利用しているnodesやtrianglesは、$n\times n$ （$n$は総節点数）の行列を前提としています。したがってディリクレ条件を適用した後はこれらを使えなくなってしまいます。したがって、ディリクレ条件の適用作業はノイマン条件の適用が終わった後に実施する必要があります。
こちらについてはノイマン条件とは異なり、境界値は節点上だけで決まっていればいいです。また、連続していないすべての境界線に対する操作をすべてまとめてapply_dirichletで処理すればOKです。
コード
これまでに説明した内容をコードに落とすとこのようになります。
- 節点の定義
- 三角形の定義
- 剛性行列の計算
- ノイマン境界条件の適用
- ディリクレ条件の適用
- 行列形式の問題（連立方程式）を解く
- 結果の描画
という順番です。
このうち、節点の定義と三角形の定義を別の条件に合わせて用意すれば、剛性行列の計算以下は同じものを使って計算できます。
（たとえば歪んだ輪のような形状とかでも計算できます）
この例では長方形を計算領域とし、以下の図ように領域分割を行っています。
import matplotlib.pyplot as plt
import numpy as np
# 節点の座標、各三角形がどの節点からなるかを指定
def generate_mesh(dx, dy, nx, ny):
nodes = []
for j in range(ny):
for i in range(nx):
x = i * dx
y = j * dy
nodes.append([x, y])
triangles = []
for j in range(ny - 1):
for i in range(nx - 1):
# 長方形の4頂点
A = j * nx + i
B = A + 1
D = (j + 1) * nx + i
C = D + 1
# 反時計回り(CCW)
triangles.append([A, B, C])
triangles.append([A, C, D])
return np.array(nodes), np.array(triangles)
# 生成した三角要素の形状、配置を確認する
def check_triangles(nodes, triangles):
color = ['blue', 'green']
for i in range(len(triangles)):
x = nodes[triangles[i],0]
y = nodes[triangles[i],1]
plt.fill(x, y, color = color[i%2])
plt.show()
# ディリクレ条件を課す境界と、境界値を設定
def set_d_bound(xlim, ylim, nodes):
upper = np.where(nodes[:,1] == ylim)[0][1:-1]
bottom = np.where(nodes[:,1] == 0)[0][1:-1]
left = np.where(nodes[:,0] == 0)[0]
right = np.where(nodes[:,0] == xlim)[0]
# ↑重複しないように注意
len_bottom = len(bottom)
len_left = len(left)
x = np.linspace(-1,1,len_left)
up_val = 2.0 *np.ones_like(upper)
bt_val = 2.0 *np.ones_like(upper)
l_val = 2.0 *np.ones_like(left)
r_val = -2*(x**2)
# ↓ノイマン条件を課す境界は出力に含めない
#bound = np.concatenate([upper, bottom, left, right])
#bd_val = np.concatenate([up_val, bt_val, l_val, r_val])
bound = np.concatenate([left, right])
bd_val = np.concatenate([l_val, r_val])
return bound, bd_val
# ノイマン条件を課す境界と、境界値を設定
def set_n_bound(xlim, ylim, nodes, select):
upper = np.sort(np.where(nodes[:,1] == ylim)[0])[1:-1]
bottom = np.sort(np.where(nodes[:,1] == 0)[0])[1:-1]
left = np.sort(np.where(nodes[:,0] == 0)[0])
right = np.sort(np.where(nodes[:,0] == xlim)[0])
len_bottom = len(bottom)
len_left = len(left)
up_val = np.linspace(0,2,len_bottom)
bt_val = np.linspace(0,-2,len_bottom)
l_val = np.linspace(-2,0,len_left)
r_val = np.linspace(-2,0,len_left)
if select == 'upper':
bound = upper
bd_val = up_val
elif select == 'bottom':
bound = bottom
bd_val = bt_val
elif select == 'left':
bound = left
bd_val = l_val
elif select == 'right':
bound = right
bd_val = r_val
return bound, bd_val
# 大域剛性行列を計算する
def calc_stiffness(nodes, triangles):
n = len(nodes)
K = np.zeros((n, n))
F = np.zeros(n)
B = np.array([[-1, 1],
[0, -1],
[1, 0]])
for i in range(len(triangles)):
nodes_i = triangles[i] #i番目の三角形の各頂点のノード番号
xs = nodes[nodes_i,0] #三角形の各頂点のx座標
ys = nodes[nodes_i,1] #同y座標
J_i = np.array([[xs[1] -xs[0], ys[1] -ys[0]],
[xs[2] -xs[0], ys[2] -ys[0]]])
BJ = B@J_i
A_i = abs(np.linalg.det(J_i))/2 #三角形の面積
K_i = (1/(4*A_i))*(BJ @ BJ.T) # 局所剛性行列
f = 10 #ポアソン方程式の右辺のｆ（単純化のため定数とした）
F_i = np.full(3, f*A_i/3)
K[np.ix_(nodes_i, nodes_i)] += K_i
F[np.ix_(nodes_i)] += F_i
return K, F
# ノイマン条件の適用
def apply_neumann(F, bound, bd_val, nodes):
for i in range(len(bound) -1):
xs = nodes[[bound[i],bound[i+1]],0]
ys = nodes[[bound[i],bound[i+1]],1]
L_i = np.sqrt((xs[1] -xs[0])**2 +(ys[1] -ys[0])**2 )
F[bound[i]] += (L_i/6)*(2*bd_val[i] +bd_val[i+1])
F[bound[i+1]] += (L_i/6)*(bd_val[i] +2*bd_val[i+1])
return F
# ディリクレ条件の適用
def apply_dirichlet(K, F, bound, bd_val):
K = np.delete(K, bound, axis=0)
F = np.delete(F, bound)
for i in range(len(bound)):
col = K[:, bound[i]]
F = F -bd_val[i]*col
K = np.delete(K, bound, axis=1)
return K, F
# 剛性行列の非ゼロ成分の分布を確認する
def checK_matrix(K):
plt.spy(K, markersize=1)
plt.show()
dx = 0.01
dy = 0.01
nx = 15
ny = 15
nodes, triangles = generate_mesh(dx, dy, nx, ny)
#check_triangles(nodes, triangles)
K, F = calc_stiffness(nodes, triangles)
# neumann条件の適用
# 境界線1本ずつやる
# dirichletの適用では行・列の削除を伴い、インデックスがずれるため
# neumann条件の適用をdirichlet条件の先に行う
# 左側の辺の境界条件を設定
u_bound, u_bd_val = set_n_bound(dx*(nx -1), dy*(ny -1), nodes, 'upper')
F = apply_neumann(F, u_bound, u_bd_val, nodes)
# 右側の辺の境界条件を設定
b_bound, b_bd_val = set_n_bound(dx*(nx -1), dy*(ny -1), nodes, 'bottom')
F = apply_neumann(F, b_bound, b_bd_val, nodes)
# dirichlet条件の適用
d_bound, d_bd_val = set_d_bound(dx*(nx -1), dy*(ny -1), nodes)
K, F = apply_dirichlet(K, F, d_bound, d_bd_val)
# 境界線が分かれていても全部の境界を一気に計算できる
#checK_matrix(K)
#print(np.linalg.cond(K))
u = np.zeros(nx*ny)
unknown_idx = np.setdiff1d(np.arange(len(u)), d_bound)
# ディリクレ条件を課していない節点の番号をunknown_idxとする
u_unknown = np.linalg.solve(K, F)
# 未知の値を求める
u[unknown_idx] = u_unknown
u[d_bound] = d_bd_val
#dirichlet条件の適用
cmap = plt.pcolormesh(u.reshape(ny, nx))
plt.colorbar(cmap, ax = plt.gca())
plt.show()
ところどころで確認作業を入れています。
まず、計算領域とその内部の分割が思った通りになっているかを確認するために、
check_triangles(nodes, triangles)
を入れています。（冒頭の図）
次に、出来上がった剛性行列の非ゼロ成分の分布を、
checK_matrix(K)
で確認しています。
また、
np.linalg.cond(K)
で、剛性行列の「条件数」を確認しています。この値が無限大（計算上は”非常に大きな値”）になっていると、逆行列が存在せず、問題を解けません。なにかミスをしていたり、境界条件の設定が不適でないか確認できます。
おわりに
以上で、二次元ポアソン方程式に対する有限要素法の基本的な実装は完成です。
もちろん、より複雑な方程式や高次要素への拡張、疎行列を利用した高速化など発展的な話題は数多くあります。しかし、そのような手法も基本的には今回実装した「局所剛性行列を計算し、大域剛性行列を構築する」という考え方の上に成り立っています。
まずは実際にコードを動かし、有限要素法による数値解の計算を試してみてください。