第4回：UFL の grad と dx の中で何が起きているか
Introduction
前回は，Gmsh の geometry_order と有限要素次数 \(k\) は別物であることを説明した．
幾何次数 \(q_{\rm geo}\) はメッシュの形をどの次数で表すかを決める量であり，有限要素次数 \(k\) は未知関数をどの次数の多項式で近似するかを決める量である．
今回は，FEniCSx/UFL で弱形式を書くときに使う grad と dx の意味を整理する．
特に，曲線要素写像
を使うとき，勾配変換に現れる \(DF_K^{-\top}\) と，積分測度に現れる \(|\det DF_K|\) が，実装上どこで扱われているのかを見る．
結論から言うと，FEniCSx/UFL では，ユーザーが明示的に Jacobian を書かなくても，grad と dx の内部で，参照要素から物理要素への変換が処理される．
物理要素上の弱形式
Poisson 問題の標準的な弱形式を考える．有限要素解 \(u_h\) と試験関数 \(v_h\) に対して，剛性行列に対応する双線形形式は
\[ a(u_h,v_h)=\int_{\Omega_h}\nabla u_h\cdot\nabla v_h\,dx \]である．要素ごとに書けば，
\[ a_K(u_h,v_h)=\int_K\nabla u_h\cdot\nabla v_h\,dx \]である．ここで \(K\) は物理要素であり，曲線要素の場合は直線三角形とは限らない． この式は，数学的には物理座標 \(x\) 上で書かれている．しかし，実際の有限要素計算では，各要素積分は参照三角形 \(\widehat T\) 上の数値積分として評価される．
参照要素への引き戻し
曲線要素 \(K\) は，要素写像
\[ F_K:\widehat T\to K \]によって参照三角形 \(\widehat T\) から作られる．参照座標を \(\widehat x\)，物理座標を \(x\) とすれば，
\[ x=F_K(\widehat x) \]である．物理要素上の積分は，この写像を使って参照要素上の積分へ変換される．すなわち，
\[ \int_K g(x)\,dx = \int_{\widehat T} g(F_K(\widehat x))\,|\det DF_K(\widehat x)|\,d\widehat x \]となる．ここで \(DF_K(\widehat x)\) は要素写像の Jacobian matrix であり，\(|\det DF_K(\widehat x)|\) は面積要素の変換を表す． したがって，物理要素上で \(dx\) と書いていても，実装では参照要素上で \(|\det DF_K|\,d\widehat x\) として評価される．
勾配はどう変換されるか
次に，基底関数の勾配を考える．参照三角形上の基底関数を \(\widehat\phi_i\) とし，物理要素上の基底関数を
\[ \phi_{K,i}=\widehat\phi_i\circ F_K^{-1} \]と定義する．つまり，
\[ \phi_{K,i}(F_K(\widehat x))=\widehat\phi_i(\widehat x) \]である．このとき，連鎖律により，物理座標での勾配は
\[ \nabla_x\phi_{K,i}(x) = DF_K(\widehat x)^{-\top} \nabla_{\widehat x}\widehat\phi_i(\widehat x), \qquad x=F_K(\widehat x) \]
となる．つまり，物理勾配 \(\nabla_x\) を計算するためには，参照勾配 \(\nabla_{\widehat x}\) に \(DF_K^{-\top}\) を掛ける必要がある．
これが，grad の中で起きている変換である．
局所剛性行列の形
以上を使うと，局所剛性行列は次のように書ける．
\[ A^K_{ij} = \int_K \nabla_x\phi_{K,j}\cdot\nabla_x\phi_{K,i}\,dx. \]これを参照要素上に引き戻すと，
\[ A^K_{ij} = \int_{\widehat T} \left( DF_K(\widehat x)^{-\top} \nabla_{\widehat x}\widehat\phi_j(\widehat x) \right) \cdot \left( DF_K(\widehat x)^{-\top} \nabla_{\widehat x}\widehat\phi_i(\widehat x) \right) |\det DF_K(\widehat x)|\,d\widehat x. \]この式を見ると，二つの重要な要素が現れている．
- \(DF_K^{-\top}\)：勾配の変換
- \(|\det DF_K|\)：積分測度の変換
曲線要素では \(DF_K\) が要素内で一定とは限らない．そのため，数値積分点ごとに \(DF_K\)，\(DF_K^{-\top}\)，\(|\det DF_K|\) が評価される．
UFL ではどう書くか
FEniCSx/UFL では，ユーザーは上のような Jacobian を直接書かずに，弱形式を物理座標で書く． たとえば，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(u) は物理座標での勾配 \(\nabla_x u_h\) を表す．
しかし，実際の assembly では，参照要素上の勾配から物理勾配へ変換される．すなわち，内部的には
が使われる．また，* dx は物理要素上の積分を表す．しかし，実際には
として参照要素上で評価される．したがって，UFL の grad と dx は，単なる記号ではなく，要素写像を通じた変換を背後に持っている．
コードに Jacobian が見えない理由
自分で有限要素コードを書く場合，局所行列を作るところで
J = ...
detJ = ...
invJT = ...
のような変数を明示的に計算することが多い．しかし，FEniCSx/UFL では，通常このような変数をユーザーが直接書かない．
その理由は，FEniCSx の form compiler が，mesh に保存された coordinate map を使って，必要な Jacobian 変換を内部的に行うからである． 曲線メッシュの場合でも，import された mesh の coordinate map が曲線幾何を表していれば，その map を通じて \(DF_K\)，\(DF_K^{-\top}\)，\(|\det DF_K|\) が評価される．
つまり，コード上で Jacobian が見えないからといって，Jacobian を使っていないわけではない．むしろ，grad と dx の評価時に内部で使われている．
曲線メッシュで特に重要な点
直線三角形の場合，\(F_K\) はアフィン写像であり，\(DF_K\) は要素内で定数である．そのため，勾配変換や面積要素の変換も比較的単純である．
一方，曲線要素の場合，\(F_K\) は一般に非アフィンである．そのため，
\[ DF_K(\widehat x),\qquad DF_K(\widehat x)^{-\top},\qquad |\det DF_K(\widehat x)| \]は参照要素内の点 \(\widehat x\) に依存する．このため，曲線要素では，quadrature point ごとに幾何情報が評価される． この違いが，曲線メッシュを扱う有限要素実装の本質である．
論文中の \(F_K=\Psi_K\circ\Phi_{T_K}\) との関係
理論では，要素写像を
\[ F_K=\Psi_K\circ\Phi_{T_K} \]と分解している．したがって，Jacobian は
\[ DF_K(\widehat x) = D\Psi_K(\Phi_{T_K}(\widehat x))\,D\Phi_{T_K} \]と書ける．この式は，勾配変換や積分測度に入る \(DF_K\) の中に，affine core のスケーリングと curved correction の両方が含まれていることを意味する．すなわち，
- \(D\Phi_{T_K}\)：アフィンなスケール，向き，異方性
- \(D\Psi_K\)：曲線幾何による補正
が合成されて，実際の要素 Jacobian \(DF_K\) になる． 実装では，この分解をユーザーが明示的に書いていなくても，mesh の coordinate map が \(F_K\) として使われる．解析では，この \(F_K\) を \(\Psi_K\circ\Phi_{T_K}\) と見て，幾何の影響を分けて評価する．
まとめ
今回は，FEniCSx/UFL の grad と dx の中で何が起きているかを整理した．
UFL では弱形式を物理座標で書くが，実際の要素積分は参照要素上で行われる．
その際，
\[ \nabla_x\phi_{K,i} = DF_K^{-\top}\nabla_{\widehat x}\widehat\phi_i \]によって勾配が変換され，
\[ dx = |\det DF_K|\,d\widehat x \]によって積分測度が変換される．
FEniCSx では，これらの処理は form compiler によって内部的に行われる．したがって，コードに Jacobian が明示的に出てこなくても，曲線要素の幾何情報は grad と dx の評価に反映されている．
次回は，実際に単位円上の Poisson 問題を設定し，\(q_{\rm geo}=1,2,3\) の曲線幾何表現を比較する． 特に，幾何誤差 \(E_{\rm area}\)，\(E_{\rm bdry}\) と有限要素誤差 \(E_{H^1}\)，\(E_{L^2}\) の違いを見る．
参考文献
この第4回では，有限要素法における参照要素への引き戻し，勾配変換，および積分測度の変換を説明した．背景として，有限要素法の標準的な組み立てと曲線要素に関する文献を挙げておく．
- P. G. Ciarlet, The Finite Element Method for Elliptic Problems, North-Holland, Amsterdam, 1978.
- C. Bernardi, Optimal finite-element interpolation on curved domains, SIAM Journal on Numerical Analysis, 26(5), 1212--1240, 1989.
- 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, Exact-curved Lagrange finite elements for the Poisson problem in two dimensions, submitted.
