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

上記の関数空間$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のパッケージを導入します。

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

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

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

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

3. 変分式を定義する

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

$$(\nabla u, \nabla v)_{\Omega} = (f,v)_{\Omega} \quad \forall v \in V$$

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

4. 有限要素解の計算

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

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

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

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} $$

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

7. 計算コードのまとめ