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日

注意

該当コードはFEniCS 2016.01以降のバージョンで実行できます。2016.01より古いFEniCSで実行する場合、不具合のことがあります。

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

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

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

1. メッシュ分割

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

メッシュの確認を確認します。

2. 有限要素空間と変分式

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

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境界値条件の設定を行います。

有限要素法は関数空間$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)$。

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

3. 有限要素解を求める

a,bを使って、有限要素解を計算できます。計算の結果はfem_uに格納されます。

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

有限要素解を直接に求める。

剛性行列(stiffness matrix)$A$とロードベクトル(load vector)の値を確認して、連立方程式を解いてから解を求めることもできます。

3Dの図を描画します。

有限要素解の誤差を計算します。

計算コードのまとめ