2次元混合有限要素法

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

劉 雪峰
最終更新:2017年9月4日
最初作成:2015年6月4日

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

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

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

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

1次混合有限要素法の関数Raviart-Thomas空間$RT$と区分的な定数区間$DG$を作ります。空間$V$は$RT$と$DG$の直和空間です。

混合法をPoisson問題の境界値問題を考える時、Dirichlet境界値条件は自然条件となって、境界における解の値の設定は不要です。

3. 変分式を定義する

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

Find $(p,u)$ in $RT \times DG$, s.t.,

$$(p,q) + (div\: q,u) =0 \quad \forall q \in RT$$$$(div\: p + f ,v) =0 \quad \forall v \in DG$$

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

上記の2つの変分式を一つにまとめてます。

Find $(p,u)$ in $RT \times DG$, s.t., $$(p,q) + (div\: q,u) + (div\: p, v) =(- f ,v) \quad \forall (q,v) \in RT \times DG$$

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

4. 有限要素解の計算

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

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

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

6. 有限要素解の誤差

以下は$(u-u_h)$の$L^2$ノルムと$(p-\nabla u)$の$L^2$ノルムを計算します。

7. 計算コードのまとめ