第5回：単位円上の Poisson 問題で曲線メッシュを比較する
Introduction
前回は，UFL の grad と dx の中で何が起きているかを整理した．
弱形式は物理座標で書いているように見えるが，実際の要素積分は参照要素上へ引き戻され，
勾配には \(DF_K^{-\top}\)，測度には \(|\det DF_K|\) が現れる，という話である．
今回は，第1期の締めくくりとして，単位円上の Poisson 問題を使い， 直線メッシュと曲線メッシュを比較する． ここでの目的は，高次有限要素法の性能を示すことではない． 有限要素次数は常に \(k=1\) に固定し，幾何次数 \(q_{\rm geo}\) だけを変える．
つまり，今回見たいのは次の点である．
\[ \text{同じ } P^1 \text{ 有限要素空間を使ったとき，幾何表現を変えると何が変わるのか．} \]この実験により，曲線幾何が主に幾何誤差を改善すること，そして \(P^1\) 近似の主要な有限要素誤差次数は自動的には変わらないことを確認する．
計算する問題
領域は単位円とする．
\[ \Omega=\{(x_1,x_2)^\top\in\mathbb R^2:x_1^2+x_2^2<1\}. \]Poisson 問題は，
\[ -\Delta u=f\quad\text{in }\Omega,\qquad u=0\quad\text{on }\partial\Omega \]である．
今回は，厳密解を
\[ u(x_1,x_2)=1-x_1^2-x_2^2 \]と選ぶ． この関数は単位円の境界上で
\[ u=0\quad\text{on }\partial\Omega \]を満たす．また，
\[ -\Delta u=4 \]である．したがって，右辺は
\[ f=4 \]とすればよい．
このように厳密解が分かっている問題を使うことで，数値解 \(u_h\) の誤差を直接計算できる．
比較するメッシュ
今回比較するのは，幾何次数
\[ q_{\rm geo}=1,2,3 \]である．有限要素次数はすべての場合で
\[ k=1 \]に固定する．
つまり，未知関数の近似空間はすべて一次 Lagrange 要素であり，変えているのはメッシュの幾何表現だけである．
- \(q_{\rm geo}=1\)：直線辺をもつ polygonal mesh
- \(q_{\rm geo}=2\)：二次の曲線幾何をもつ mesh
- \(q_{\rm geo}=3\)：三次の曲線幾何をもつ mesh
\(q_{\rm geo}=1\) では，円周は直線の chord で近似される． そのため，計算領域は厳密な単位円ではなく，多角形領域になる． 一方，\(q_{\rm geo}=2,3\) では，境界辺の途中に幾何ノードが入り，境界を曲線的に表現できる．
この違いが，面積誤差や境界誤差に大きく反映される．
Gmsh/FEniCSx での基本方針
Gmsh 側では，メッシュ生成後に geometry order を指定する． 概念的には次の部分である．
gmsh.model.mesh.generate(2)
gmsh.model.mesh.setOrder(geometry_order)
ここで geometry_order が \(q_{\rm geo}\) に対応する．
一方，FEniCSx 側では，有限要素空間を一次 Lagrange 空間として定義する．
V = fem.functionspace(mesh, ("Lagrange", 1))
この 1 が有限要素次数 \(k=1\) に対応する．
したがって，たとえ geometry_order が 2 や 3 であっても，未知関数の近似空間は一次 Lagrange のままである．
UFL では，Poisson 問題の弱形式を次のように書く．
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
a = ufl.inner(ufl.grad(u), ufl.grad(v)) * dx
L = f * v * dx
ここで ufl.grad の中では物理座標での勾配が扱われ，
dx の中では曲線幾何に対応する積分測度が扱われる．
実際の要素積分では，参照要素上に引き戻され，\(DF_K^{-\top}\) と \(|\det DF_K|\) が内部的に使われる．
幾何誤差
まず，メッシュの幾何表現そのものを調べるため，二つの幾何誤差を計算する．
一つ目は，面積誤差である．
\[ E_{\rm area}:=\bigl||\Omega_h|-\pi\bigr|. \]ここで \(|\Omega_h|\) は，計算メッシュが表す領域の面積である． 単位円の面積は \(\pi\) なので，この差を見れば，領域の面積がどれだけ正確に表されているかが分かる．
二つ目は，境界半径誤差である．
\[ E_{\rm bdry}:= \max_{x_q\in\partial\Omega_h} \bigl||x_q|-1\bigr|. \]これは，境界上のサンプル点 \(x_q\) が，単位円の半径 1 からどれだけずれているかを測る量である． 厳密な最大値ではなく，境界 coordinate map から得たサンプル点に対する診断量として用いる．
有限要素誤差
次に，有限要素解 \(u_h\) の相対誤差を計算する． \(H^1\)-seminorm 相対誤差は，
\[ E_{H^1} := \left( \sum_{K\in\mathcal K_h} \int_K |\nabla(u-u_h)|^2\,dx \right)^{1/2} /|u|_{H^1(\Omega)} \]で定義する．
また，\(L^2\)-norm 相対誤差は，
\[ E_{L^2} := \left( \sum_{K\in\mathcal K_h} \int_K |u-u_h|^2\,dx \right)^{1/2} /\|u\|_{L^2(\Omega)} \]で定義する．
観測次数は，連続する二つのメッシュに対して
\[ {\rm rate} = \frac{\log(E_{\rm old}/E_{\rm new})} {\log(h_{\rm old}/h_{\rm new})} \]で計算する．
一次 Lagrange 要素を用いているので，滑らかな解に対しては標準的に
\[ E_{H^1}=O(h),\qquad E_{L^2}=O(h^2) \]が期待される．
幾何誤差の結果
まず，単位円の幾何誤差を比較する．
| \(q_{\rm geo}\) | level | \(h\) | \(E_{\rm area}\) | \(E_{\rm bdry}\) |
|---|---|---|---|---|
| 1 | 0 | \(4.000\times10^{-1}\) | \(8.013\times10^{-2}\) | \(1.916\times10^{-2}\) |
| 1 | 1 | \(2.000\times10^{-1}\) | \(2.015\times10^{-2}\) | \(4.802\times10^{-3}\) |
| 1 | 2 | \(1.000\times10^{-1}\) | \(5.205\times10^{-3}\) | \(1.240\times10^{-3}\) |
| 1 | 3 | \(5.000\times10^{-2}\) | \(1.302\times10^{-3}\) | \(3.100\times10^{-4}\) |
| 2 | 0 | \(4.000\times10^{-1}\) | \(1.549\times10^{-4}\) | \(4.596\times10^{-5}\) |
| 2 | 1 | \(2.000\times10^{-1}\) | \(9.717\times10^{-6}\) | \(2.887\times10^{-6}\) |
| 2 | 2 | \(1.000\times10^{-1}\) | \(6.473\times10^{-7}\) | \(1.924\times10^{-7}\) |
| 2 | 3 | \(5.000\times10^{-2}\) | \(4.047\times10^{-8}\) | \(1.203\times10^{-8}\) |
| 3 | 0 | \(4.000\times10^{-1}\) | \(2.282\times10^{-5}\) | \(1.184\times10^{-5}\) |
| 3 | 1 | \(2.000\times10^{-1}\) | \(1.437\times10^{-6}\) | \(7.451\times10^{-7}\) |
| 3 | 2 | \(1.000\times10^{-1}\) | \(9.587\times10^{-8}\) | \(4.968\times10^{-8}\) |
| 3 | 3 | \(5.000\times10^{-2}\) | \(5.995\times10^{-9}\) | \(3.107\times10^{-9}\) |
\(q_{\rm geo}=1\) では，境界は直線 chord で表されるため，面積誤差も境界誤差も比較的大きい． 一方，\(q_{\rm geo}=2,3\) では，これらの誤差が大きく小さくなる． これは，計算が単なる多角形メッシュではなく，曲線幾何をもつメッシュ上で行われていることを示している．
有限要素誤差の結果
次に，Poisson 問題の相対有限要素誤差を見る．
| \(q_{\rm geo}\) | level | \(h\) | \(E_{H^1}\) | rate | \(E_{L^2}\) | rate |
|---|---|---|---|---|---|---|
| 1 | 0 | \(4.000\times10^{-1}\) | \(1.465\times10^{-1}\) | -- | \(6.316\times10^{-2}\) | -- |
| 1 | 1 | \(2.000\times10^{-1}\) | \(7.746\times10^{-2}\) | 0.92 | \(1.674\times10^{-2}\) | 1.92 |
| 1 | 2 | \(1.000\times10^{-1}\) | \(4.044\times10^{-2}\) | 0.94 | \(4.426\times10^{-3}\) | 1.92 |
| 1 | 3 | \(5.000\times10^{-2}\) | \(2.031\times10^{-2}\) | 0.99 | \(1.111\times10^{-3}\) | 1.99 |
| 2 | 0 | \(4.000\times10^{-1}\) | \(1.480\times10^{-1}\) | -- | \(2.598\times10^{-2}\) | -- |
| 2 | 1 | \(2.000\times10^{-1}\) | \(7.800\times10^{-2}\) | 0.92 | \(6.940\times10^{-3}\) | 1.90 |
| 2 | 2 | \(1.000\times10^{-1}\) | \(4.057\times10^{-2}\) | 0.94 | \(1.837\times10^{-3}\) | 1.92 |
| 2 | 3 | \(5.000\times10^{-2}\) | \(2.036\times10^{-2}\) | 0.99 | \(4.600\times10^{-4}\) | 2.00 |
| 3 | 0 | \(4.000\times10^{-1}\) | \(1.426\times10^{-1}\) | -- | \(2.447\times10^{-2}\) | -- |
| 3 | 1 | \(2.000\times10^{-1}\) | \(7.665\times10^{-2}\) | 0.90 | \(6.737\times10^{-3}\) | 1.86 |
| 3 | 2 | \(1.000\times10^{-1}\) | \(4.023\times10^{-2}\) | 0.93 | \(1.811\times10^{-3}\) | 1.89 |
| 3 | 3 | \(5.000\times10^{-2}\) | \(2.027\times10^{-2}\) | 0.99 | \(4.565\times10^{-4}\) | 1.99 |
すべての幾何次数で，\(H^1\)-seminorm 誤差はおおむね一次，\(L^2\)-norm 誤差はおおむね二次に近づいている． これは，理論的に期待される \(P^1\) 有限要素法の収束挙動と整合している．
一方で，\(q_{\rm geo}=1,2,3\) の \(H^1\) 誤差はかなり似ている． これは不思議ではない． なぜなら，今回変えているのは幾何次数だけであり，有限要素次数は常に \(k=1\) だからである．
ただし，\(L^2\) 誤差については，\(q_{\rm geo}=1\) と \(q_{\rm geo}=2,3\) の間にある程度の差が見られる． これは，\(L^2\) 誤差が \(H^1\) 誤差よりも幾何表現の影響を受けやすいことを示唆している． 特に，直線辺による polygonal representation では，領域そのものが単位円からずれているため， 幾何誤差が \(L^2\) 誤差にも反映されやすい．
一方で，\(q_{\rm geo}=2\) と \(q_{\rm geo}=3\) の差はかなり小さい． この実験では有限要素空間を \(P^1\) に固定しているため， 幾何表現を二次以上に改善した後は，主要な有限要素誤差は \(P^1\) 近似によって支配されると考えられる．
結果の読み方
今回の結果から，次のことが分かる．
今回の実験では，geometry order を上げることで面積誤差と境界誤差は大きく改善する． また，\(L^2\) 誤差にも \(q_{\rm geo}=1\) と \(q_{\rm geo}=2,3\) の間で一定の改善が見られる． しかし，有限要素空間を \(P^1\) に固定しているため，\(H^1\) 誤差の主要な振る舞いは大きく変わらない． つまり，曲線幾何は幾何誤差を減らし，\(L^2\) 誤差にも影響を与えるが， 有限要素近似次数そのものを高次化するわけではない．曲線幾何は，主に「領域をどれだけ正確に表すか」に効いている．有限要素近似そのものを高次化したいなら，finite element degree \(k\) も上げる必要がある．
特に重要なのは，
\[ \text{曲線幾何を高次にすること} \neq \text{有限要素近似を高次にすること} \]である．
この区別が曖昧だと，geometry order を上げても \(H^1\) 誤差が劇的に改善しないことを不自然に感じるかもしれない． しかし，今回の実験では，未知関数の近似空間はあくまで \(P^1\) に固定されている． したがって，主要な \(H^1\) 誤差が \(P^1\) 近似に支配されるのは自然である．
第1期のまとめ
第1期では，exact-curved FEM の基本的な考え方と，単位円上の Poisson 問題での実装実験を整理した．
- 曲線領域を多角形近似せず，曲線のまま扱うという問題意識
- 要素写像 \(F_K=\Psi_K\circ\Phi_{T_K}\) の分解
- geometry order と finite element degree の違い
- UFL の
gradとdxの中での Jacobian 変換 - 単位円上の Poisson 問題による数値確認
これで，二次元の \(P^1\) Lagrange 要素に対する exact-curved FEM の原型が見えてきた． 次の段階では，高次 Lagrange 要素や三次元への拡張を考えることになる．
参考文献
この第5回では，投稿中論文の数値実験部分をもとに，単位円上の Poisson 問題で直線メッシュと曲線メッシュを比較した． 背景となる文献として，曲線要素上の補間理論，有限要素法の標準的文献，および筆者の関連研究を挙げておく．
- C. Bernardi, Optimal finite-element interpolation on curved domains, SIAM Journal on Numerical Analysis, 26(5), 1212--1240, 1989.
- P. G. Ciarlet, The Finite Element Method for Elliptic Problems, North-Holland, Amsterdam, 1978.
- P. G. Ciarlet and P.-A. Raviart, Interpolation theory over curved elements, with applications to finite element methods, Computer Methods in Applied Mechanics and Engineering, 1, 217--249, 1972.
- H. Ishizaka, K. Kobayashi and T. Tsuchiya, Anisotropic interpolation error estimates using a new geometric parameter, Japan Journal of Industrial and Applied Mathematics, 40, 475--512, 2023.
- H. Ishizaka, Exact-curved Lagrange finite elements for the Poisson problem in two dimensions, submitted.
