👨‍👩‍👧 ♑️ ♌️ 有限要素プログラミング 🎓 🛀 👩‍🔬

この記事は、拡散反応方程式系の有限要素法の独自の実装（ソルバージョーカーFEM ）に特化しています。

通常、既製のソリューションを使用することをお勧めしますが、問題に特定の機能がある場合は、単純なライブラリに基づいて、タスクを解決するのが簡単です。

はじめに

有限要素法[1-4]は、数理物理学、連続体力学の問題を解決するための効果的な方法のグループです。 FEMプログラミングの質問は、ソース[5–9]に当てられています。

問題の声明

システムのロビン境界値問題 $N$ 限られた領域での拡散反応方程式 $\オメガ\サブセット\ mathbb R ^ 2$ ボーダー付き $\ガンマ$ ：

- a_{i} D e l t a u_{i} （ x ） + s u m_{k = 1}^{N} q_{i k} （ x ） u_{k} （ x ） = g_{i} （ x ） 、 x i n O m e g a 、

$-a_i \ Delta u_i（x）+ \ sum_ {k = 1} ^ N q_ {ik}（x）u_k（x）= g_i（x）、\; \; x \ in \ Omega、$

a_{i} f r a c p a r t i a l u_{i} （ x ） p a r t i a l m a t h b f n + b_{i} （ x ） u_{i} （ x ） = b_{i} （ x ） w_{i} （ x ） 、 x i n \ガ ン マ 。

$a_i \ frac {\ partial u_i（x）} {\ partial \ mathbf n} + b_i（x）u_i（x）= b_i（x）w_i（x）、\; \; x \ in \ガンマ。$

ここに $i = 1、2、\ ldots、N$ 、 $\ mathbf n$ 外部法線のベクトルは $\ガンマ$ 。

弱い表現

それぞれを掛ける $N$ 任意のテスト関数の方程式 $v$ 、取得したアイデンティティを統合します $\オメガ$ 、グリーンの式を適用し、境界条件を考慮します。ゲット

a_{i} i n t l i m i t s_{} O m e g a n a b l a u_{i} n a b l a v \、 d x + i n t l i m i t s_{} G a m m a b_{i} u_{i} v \、 d G a m m a + s u m_{k = 1}^{N} i n t l i m i t s_{} O m e g a q_{i k} u_{k} v \、 d x = i n t l i m i t s_{} O m e g a g_{i} v \、 d x + i n t l i m i t s_{} G a m m a b_{i} w_{i} v \、 d G a m m a 、 i = 1 、 l d o t s 、 N .

$a_i \ int \ limits_ \ Omega \ nabla u_i \ nabla v \、dx + \ int \ limits_ \ Gamma b_i u_i v \、d \ Gamma + \ sum_ {k = 1} ^ N \ int \ limits_ \ Omega q_ { ik} u_k v \、dx = \ int \ limits_ \ Omega g_i v \、dx + \ int \ limits_ \ Gamma b_i w_i v \、d \ Gamma、\; \; i = 1、\ ldots、N.$

機能的なソボレフ空間を紹介します $V = H ^ 1（\オメガ）$ 。機能 $u_1、\ ldots、u_N \ in V$ これらの恒等式が任意のテスト関数に当てはまる場合、元の境界値問題の弱解と呼ばれます $v \ in V$ 。

弱い表現は次のように書くことができます：関数を見つける $u_1、\ ldots、u_N \ in V$ そのような

A_{i} （ u_{i} 、 v ） + s u m_{k = 1}^{N} B_{i k} （ u_{k} 、 v ） = F_{i} （ v ） f o r a l l v i n V 、 i = 1 、 l d o t s 、 N 、

$A_i（u_i、v）+ \ sum_ {k = 1} ^ N B_ {ik}（u_k、v）= F_i（v）\; \; \ forall v \ in V、\; i = 1、\ ldots、N、$

どこで $A_i、B_ {ik}$ -バイリニア形式、 $F_i$ -線形形式：

A_{i} （ u 、 v ） = a_{i} i n t l i m i t s_{} O m e g a n a b l a u n a b l a v \、 d x + i n t l i m i t s_{} G a m m a b_{i} u v \、 d G a m m a 、 q u a d B_{i k} （ u 、 v ） = s u m_{k = 1}^{N} i n t l i m i t s_{} O m e g a q_{i k} u v \、 d x 、

$A_i（u、v）= a_i \ int \ limits_ \ Omega \ nabla u \ nabla v \、dx + \ int \ limits_ \ Gamma b_i uv \、d \ Gamma、\ quad B_ {ik}（u、v） = \ sum_ {k = 1} ^ N \ int \ limits_ \ Omega q_ {ik} uv \、dx、$

F_{i} （ v ） = i n t l i m i t s_{} O m e g a g_{i} v \、 d x + i n t l i m i t s_{} G a m m a b_{i} w_{i} v \、 d G a m m a 。

$F_i（v）= \ int \ limits_ \ Omega g_i v \、dx + \ int \ limits_ \ Gamma b_i w_i v \、d \ Gamma。$

FEMアルゴリズム

宇宙で紹介します $V$ 有限次元部分空間 $V_m$ 。ガラーキン法によると、離散問題に到達します：関数を見つける $u_1、\ ldots、u_N \ in V_m$ そのような

A_{i} （ u_{i} 、 v ） + s u m_{k = 1}^{N} B_{i k} （ u_{k} 、 v ） = F_{i} （ v ） f o r a l l v V_{m} で 、 i = 1 、 l d o t s 、 N 。

$A_i（u_i、v）+ \ sum_ {k = 1} ^ N B_ {ik}（u_k、v）= F_i（v）\; \; \ forall v \ V_mで、\; i = 1、\ ldots、N。$

させる $\ phi_1、\ ldots、\ phi_m$ -基礎 $V_m$ 。入れたら $u_i（x）= \ displaystyle \ sum_ {s = 1} ^ m c_ {is} \ phi_s（x）$ 、 $v（x）= \ phi_j（x）$ 、それから問題はSLAE $Nm$ 不明 $c_ {ks}$ ：

s u m_{s = 1}^{m} A_{i} （ p h i_{s} 、 p h i_{j} ） c_{i s} + s u m_{k = 1}^{N} s u m_{s = 1}^{m} B_{i k} （ p h i_{s} 、 p h i_{j} ） c_{k s} = F_{i} （ p h i_{j} ） 、 i = 1 、 l d o t s 、 N 、 j = 1 、 l d o t s 、 m 、

$\ sum_ {s = 1} ^ m A_i（\ phi_s、\ phi_j）c_ {is} + \ sum_ {k = 1} ^ N \ sum_ {s = 1} ^ m B_ {ik}（\ phi_s、\ phi_j）c_ {ks} = F_i（\ phi_j）、\; \; i = 1、\ ldots、N、\; j = 1、\ ldots、m、$

または

\ sum_ {s = 1 c_ {is} + \ sum_ {k = 1} ^ N \ sum_ {s = 1} ^ m \ left（\ int_ \ Omega q_ {ik} \ phi_s \ phi_j \、dx \ right）c_ {ks} =

$\ sum_ {s = 1 c_ {is} + \ sum_ {k = 1} ^ N \ sum_ {s = 1} ^ m \ left（\ int_ \ Omega q_ {ik} \ phi_s \ phi_j \、dx \ right）c_ {ks} =$

= i n t l i m i t s_{} O m e g a g_{i} p h i_{j} \、 d x + i n t l i m i t s_{} G a m m a b_{i} w_{i} p h i_{j} \、 d G a m m a 、 i = 1 、 l d o t s 、 N 、 j = 1 、 l d o t s 、 m

$= \ int \ limits_ \ Omega g_i \ phi_j \、dx + \ int \ limits_ \ Gamma b_i w_i \ phi_j \、d \ Gamma、\; \; i = 1、\ ldots、N、\; j = 1、\ ldots、m$

部分空間の最も単純な例 $V_m$ -連続的な区分的線形関数の空間。面積 $\オメガ$ 有限要素（通常は三角形）に分割され、各三角形の関数 $V_m$ 線形。 FEMでは、基底関数は、ほとんどの領域でゼロになるように選択され、SLAE行列には多くのゼロがあります。基本機能として $\ phi_j$ 1つのグリッドノードで1、残りのノードで0に等しい「ピラミッド」を取得します。各グリッドノードには独自の基底関数があり、ディメンション $m$ 部分空間はグリッドノードの数に等しくなります。

ことに注意してください $s$ そして $j$ ノードがエッジで接続されていない場合、 $\ varphi_s \ varphi_j \ equiv 0$ 、 $\ nabla \ varphi_s \ cdot \ nabla \ varphi_j \ equiv 0$ 、したがって、合計すると $s$ インデックスのみを反復処理するのに十分 $s$ その間 $s$ そして $j$ 番目のノードはエッジを $s = j$ 。

したがって、有限要素法による境界値問題の解は、対応するSLAEの構築と解に還元されます。

区分線形関数

機能を設定するには $u \ in V_m$ 、グリッドのノードでその値を指定する必要があります。座標を示します $i$ 番目のノード $（X_i、Y_i）$ 、および関数の値 $u$ で $i$ 番目のノード-スルー $U_i = u（X_i、Y_i）$ 。

関数の値を見つけるには $u$ 任意の点で $（x、y）$ との三角形に属する $i$ 番目 $j$ そして $k$ 頂点、重心座標を導入します $L_1、L_2、L_3$ 0から1の値を取り、デカルト座標に関連付けられています $x、y$ 次の関係（[3、セクション4.7.1]、[4、p。35]を参照）：

b e g i n g a t h e r e d x = L_{1} X_{i} + L_{2} X_{j} + L_{3} X_{k} 、 y = L_{1} Y_{i} + L_{2} Y_{j} + L_{3} Y_{k} 、 1 = L_{1} + L_{2} + L_{3} 。 e n d g a t h e r e d

$\ begin {gathered} x = L_1 X_i + L_2 X_j + L_3 X_k、\\ y = L_1 Y_i + L_2 Y_j + L_3 Y_k、\\ 1 = L_1 + L_2 + L_3。 \ end {gathered}$

ここから

b e g i n g a t h e r e d L_{1} （ x 、 y ） = f r a c 1 2 A （ a_{i} + b_{i} x + c_{i} y ） 、 a_{i} = X_{j} Y_{k} - X_{k} Y_{j} 、 b_{i} = Y_{j} - Y_{k} 、 c_{i} = X_{k} - X_{j} 、 L_{2} （ x 、 y ） = f r a c 1 2 A （ a_{j} + b_{j} x + c_{j} y ） 、 a_{j} = X_{k} Y_{i} - X_{i} Y_{k} 、 b_{j} = Y_{k} - Y_{i} 、 c_{j} = X_{i} - X_{k} 、 L_{3} （ x 、 y ） = f r a c 1 2 A （ a_{k} + b_{k} x + c_{k} y ） 、 a_{k} = X_{i} Y_{j} - X_{j} Y_{i} 、 b_{k} = Y_{i} - Y_{j} 、 c_{k} = X_{j} - X_{i} 、 e n d g a t h e r e d

$\ begin {gathered} L_1（x、y）= \ frac {1} {2A}（a_i + b_i x + c_i y）、\; a_i = X_j Y_k-X_k Y_j、\; b_i = Y_j-Y_k、\; c_i = X_k-X_j、\\ L_2（x、y）= \ frac {1} {2A}（a_j + b_j x + c_j y）、\; a_j = X_k Y_i-X_i Y_k、\; b_j = Y_k-Y_i、\; c_j = X_i-X_k、\\ L_3（x、y）= \ frac {1} {2A}（a_k + b_k x + c_k y）、\; a_k = X_i Y_j-X_j Y_i、\; b_k = Y_i-Y_j、\; c_k = X_j-X_i、\\ \ end {gathered}$

どこで $A = \ dfrac {1} {2} \ begin {vmatrix} 1＆X_i＆Y_i \\ 1＆X_j＆Y_j \\ 1＆X_k＆Y_k \ end {vmatrix}$ -記号の付いた三角形の領域。

関数の値は、式によって見つけることができます $u（x、y）= U_i L_1（x、y）+ U_j L_2（x、y）+ U_k L_3（x、y）$ 。

数値積分

境界エッジ上の積分を $i$ そして $j$ セグメント上の積分に終わる $[-1、1]$ ：

i n t l i m i t s_{i}^{j} f （ x 、 y ） \、 d l = f r a c | i \、 j | 2 i n t l i m i t s_{- 1}^{1} f （ x （ t ） 、 y （ t ） ） \、 d t 。

$\ int \ limits_i ^ jf（x、y）\、dl = \ frac {| i \、j |} {2} \ int \ limits _ {-1} ^ 1 f（x（t）、y（t））\、dt。$

ここに $| i \、j |$ -リブの長さ $x（t）= \ dfrac {X_i + X_j} {2} + \ dfrac {X_j-X_i} {2} t$ 、
$y（t）= \ dfrac {Y_i + Y_j} {2} + \ dfrac {Y_j-Y_i} {2} t$ 。

最後の積分を計算するには、ガウス求積式を適用します。

i n t l i m i t s_{- 1}^{1} v a r p h i （ s ） \、 d s \お よ そ s u m_{j = 1}^{n} w_{j} v a r p h i （ x i_{j} ）

$\ int \ limits _ {-1} ^ 1 \ varphi（s）\、ds \およそ\ sum_ {j = 1} ^ n w_j \ varphi（\ xi_j）$

ポイント $\ xi_j \ in [-1、1]$ および重み係数 $w_j$ テーブルから取られた[3、秒。 5.9.2]。

三角形の積分を $i$ 番目 $j$ そして $k$ 三角形上の積分への頂点 $D = \ {（L_1、L_2）\ in [0,1] \ times [0,1] \ colon L_1 + L_2 \ leq 1 \}$ ：

i n t l i m i t s_{i \、 j \、 k} f （ x 、 y ） \、 d x d y = 2 | i \、 j \、 k | i n t l i m i t s_{D} f （ x （ L_{1} 、 L_{2} ） 、 y （ L_{1} 、 L_{2} ） ） \、 d L_{1} d L_{2} 。

$\ int \ limits_ {i \、j \、k} f（x、y）\、dxdy = 2 | i \、j \、k | \ int \ limits_D f（x（L_1、L_2）、y（L_1、L_2））\、dL_1dL_2。$

ここに $| i \、j \、k | = | A |$ 三角形の面積です $x（L_1、L_2）= L_1X_i + L_2X_j +（1-L_1-L_2）X_k$ 、
$y（L_1、L_2）= L_1Y_i + L_2Y_j +（1-L_1-L_2）Y_k$ 。

最後の積分を計算するために、求積式[3、Sec。 5.11]。

ソフトウェア実装

一般的な考慮事項

このサブジェクト領域では、テストが難しい[10]ため、コードを非常に明確に記述することが重要です。 OOPを使用すると、コードを読みやすくすることができますOOPを使用するには、一連の概念を選択し、これらの用語でアルゴリズムを策定する必要があります。これにより、まずアルゴリズムに注意を払い、次に実装の詳細にのみ注意を払うことができます。

三角測量

グリッドを構築するには、 FreeFem ++パッケージが使用され、三角測量は.mesh形式で保存されます。

ライブラリには、ファイルから読み取られるグリッドに関する情報を含むMeshクラスが導入されています。

SLAUの構築

アルゴリズムを記述するために、基底関数と被積分関数（関数の積）のクラスが導入され、これらは基底クラスIntegrandおよびBoundaryIntegrandから継承されます。これらのクラスには、それぞれフィールドsupport （media）およびboundary_support （境界メディア）があります。キャリアは、関数が0に等しくない三角形のリストです。境界キャリアは、関数が0に等しくない境界エッジのリストですValueメソッドは、被積分関数のクラスに対しても定義され、特定のポイントで関数の値を計算します。

IntegrateおよびBoundaryIntegrate関数は、被積分関数を入力として受け取り、ホストリストを反復処理し、 Integrand基本クラス（ BoundaryIntegrand ）で定義された特定の三角形（または境界エッジ）で積分関数を呼び出し、仮想関数Valueを呼び出します。

したがって、SLAU構築コードの形式は次のとおりです。

 for (int i = 0; i < data.N; ++i) { for (int j = 0; j < data.mesh.nodes_num; ++j) { for (auto s : data.mesh.adjacent_nodes[j]) { sys.AddCoeff(i, j, i, s, data.a[i] * Integrate(grad(phi[s]) * grad(phi[j])) + BoundaryIntegrate(data.b[i] * phi[s] * phi[j]) ); } for (int k = 0; k < data.N; ++k) { for (auto s : data.mesh.adjacent_nodes[j]) { sys.AddCoeff(i, j, k, s, Integrate(data.q[i][k] * phi[s] * phi[j]) ); } } sys.AddRhs(i, j, Integrate(data.g[i] * phi[j]) + BoundaryIntegrate(data.b[i] * data.w[i] * phi[j]) ); } }

SLAUソリューション

SLAEを解決するには、 MTL4ライブラリに実装されている反復法が使用されます。

結果へのアクセス

特定のポイントで解の値を見つけるには、このポイントを含む三角形を短時間で見つける必要があります。これにはハッシュテーブルを使用できます。

領域を囲む長方形を破る $\オメガ$ に $nm$ ブロックどこ $n = [W / a]$ 、 $m = [H / b]$ 、 $W、H$ -囲む長方形の長さと幅、 $a、b$ -三角形分割三角形を囲む長方形の平均寸法。

それぞれの $nm$ ブロックは、交差する三角形のリストを塗りつぶします。これを行うには、すべての三角形をループして、それらを囲む長方形を見つけます。
交差するブロックを簡単に見つけることができます。

特定のポイントに対してクエリを実行すると、それが属するブロックを見つけ、このブロックに対応する三角形のリストを調べて、特定のポイントがどのブロックに属しているかを確認します。

テスト中

プログラムのテストは、既知の正確なソリューションを使用してタスクで実行されました。ダブルメッシュリファインメントでは、誤差は約4倍または2倍減少します。これは、メソッドの精度が2次または1次であることを示しています。 1番目への収束の順序の減少は、関数の異なる値が与えられる境界部分の接合部での境界条件の不一致による可能性があります $w_i$ 。

おわりに

開発されたライブラリの大きなソルバーに対する利点は、そのシンプルさと、タスクの特定の機能（境界条件での区分的に一定の係数）を考慮できることです。将来、3次元問題用の同様のライブラリを実装する予定です。このライブラリは、複雑な熱伝達のモデルで境界係数の最適制御の問題を解決するために使用されます。

ソース

アレクシーエフG.V. 数理物理学の数値的方法。有限要素法の紹介。ウラジオストク：ダルナウカ、1999年。
ダウトフR.Z.、カルチェフスキーM.M. 有限要素法の理論の紹介。カザン：KSU、2004年。
Zienkiewicz OC、Taylor RL、Zhu JZ有限要素法：その基礎と基礎。エルゼビア、2005年。
Segerlind L.有限要素法の適用。 M。：ミール、1979年。
ポアソン方程式の例の有限要素法 // Habrahabr、07.27.2015。
180行未満のコードでFEM電卓を作成する //
Habrahabr、2015年12月1日。
Freefem ++の使用方法 // Habrahabr、2013年6月6日。
Gockenbach MS有限要素法の理解と実装。 SIAM、2006年。
Smith IM、Griffiths DV、Margetts L.有限要素法のプログラミング。ワイリー、2014年。
科学的なアプリケーションでユニットテストが機能しない理由 // Habrahabr、2010年4月26日。
グレンキンG.V. 複雑な熱伝達モデルの境界最適制御問題を解決するためのアルゴリズム // Dalnevost。仲間。ジャーナル 2016.V. 16. No. 1. P. 24-38。

有限要素プログラミング