第1回：有限要素法で「丸い領域」を丸いまま扱いたい

Introduction

有限要素法で円領域や曲線境界をもつ領域を扱うとき，多くの場合，まず領域を三角形や四面体で近似する．二次元の単位円であれば，円周を直線の辺でつないだ多角形領域として近似し，その上で有限要素法を適用する．これは実装が簡単であり，標準的な方法でもある．

しかし，このとき計算している問題は，厳密には元の円領域上の問題ではなく，多角形近似された領域上の問題である．つまり，有限要素近似による誤差だけでなく，領域を多角形で近似したことによる幾何誤差も含まれる．

そこで自然に出てくる問いは次である．

有限要素法で，曲がった領域を本当に曲がったまま扱うことはできないのか？

このシリーズでは，この問いを出発点として，exact-curved finite element method の理論と実装を少しずつ整理していく． 第1期では，二次元の Poisson 問題を題材に，曲線境界をもつ領域上の Lagrange 有限要素法を考える．

この記事は，投稿中の論文 Exact-curved Lagrange finite elements for the Poisson problem in two dimensions の数値実験部分を，実装ノートとして補足するシリーズの第1回である．論文本文では理論と結果を中心にまとめたため，ブログでは，Gmsh/FEniCSx で曲線幾何をどう扱うか，geometry order と有限要素次数の違いは何か，UFL の grad と dx の内部で何が起きているかを，できるだけ丁寧に記録していく．

直線メッシュと曲線メッシュ

たとえば，単位円 \[ \Omega=\{(x_1,x_2)^{\top}\in\mathbb R^2:x_1^2+x_2^2<1\} \] を考える．通常の三角形メッシュでは，円周は直線の辺で近似される．そのため，計算領域は円ではなく，多角形になる．

一方で，曲線メッシュを用いると，境界辺そのものを曲線として表すことができる．このとき，境界が円周に沿うように要素写像を作れば，計算領域はより正確に元の曲線領域を表すことができる．

ここで重要なのは，曲線メッシュを使うことと，高次有限要素を使うことは同じではない，という点である．メッシュの幾何をどれだけ正確に表すかという問題と，未知関数をどの次数の多項式で近似するかという問題は，別の問題である．

geometry order と finite element degree は違う

今回の数値実験では，有限要素空間は常に一次 Lagrange 要素，つまり \(P^1\) 要素に固定する．その上で，幾何表現の次数だけを変える．

記号で書けば，有限要素次数は \[ k=1 \] に固定し，幾何次数を \[ q_{\rm geo}=1,2,3 \] と変える．

\(q_{\rm geo}=1\) の場合，境界は直線辺で表される．これは通常の polygonal mesh に対応する．一方，\(q_{\rm geo}=2,3\) とすると，境界や要素の幾何表現が曲線的になり，単位円をより正確に表すことができる．

ただし，未知関数の近似空間は \(P^1\) のままである．したがって，幾何表現を高次にしても，有限要素近似そのものが自動的に高次になるわけではない．この点は，数値結果を読むときに非常に重要である．

理論で使う要素写像

曲線要素上の有限要素法には，Ciarlet--Raviart，Zlámal，Scott，Lenoir らによる古典的研究がある．また，Bernardi は曲線領域に正確に適合する有限要素補間理論を展開している．このシリーズで扱う exact-curved FEM は，これらの流れを背景にしつつ，要素写像を affine core と curved correction に分けて見直す立場から出発する．

このシリーズの中心になるのは，要素写像の分解である．曲がった要素 \(K\) を扱うために，参照三角形 \(\widehat T\) から物理要素 \(K\) への写像を \[ F_K=\Psi_K\circ\Phi_{T_K} \] と書く．

ここで，\(\Phi_{T_K}\) は参照三角形から直線三角形である affine core \(T_K\) への写像である．一方，\(\Psi_K\) は affine core \(T_K\) を曲要素 \(K\) へ移す曲線補正である．

つまり，

  • \(\Phi_{T_K}\)：三角形のスケール，向き，異方性を表す部分
  • \(\Psi_K\)：曲線境界や曲がった幾何を表す部分

と分けて考える． この分解により，有限要素補間の誤差を，まず affine core 上で評価し，その後，曲線写像 \(\Psi_K\) を通じて物理要素へ移すことができる．

実装では何が起きているのか

実装では，解析的な \(\Psi_K\) を自分で完全に書き下しているわけではない．Gmsh で曲線メッシュを作り，それを FEniCSx に読み込むと，メッシュの中に coordinate map が保存される．この coordinate map が，数値積分や基底関数の評価の段階で，要素写像 \(F_K\) の役割を果たす．

FEniCSx/UFL では，弱形式を物理座標で書く．たとえば Poisson 問題なら， \[ \int_{\Omega_h}\nabla u_h\cdot\nabla v_h\,dx \] のように書く．しかし，実際の要素積分は参照要素上で行われる．そのとき，勾配には \(DF_K^{-\top}\)，積分測度には \(|\det DF_K|\) が内部的に現れる．

つまり，コード上で明示的に Jacobian を書いていなくても，UFL の grad と dx の内部で，参照要素から物理要素への変換が処理されている．この点は，後の回で詳しく見る．

このシリーズで扱うこと

第1期では，次の内容を順に整理する予定である．

  1. 有限要素法で「丸い領域」を丸いまま扱いたい
  2. 曲線要素写像 \(F_K=\Psi_K\circ\Phi_{T_K}\) の意味
  3. Gmsh の geometry order と finite element degree の違い
  4. FEniCSx/UFL における Jacobian，grad，dx の意味
  5. 単位円上の Poisson 問題で，直線メッシュと曲線メッシュを比較する

その後，高次 Lagrange 要素，三次元の曲四面体，非適合要素，混合型要素，さらに Stokes 問題への拡張を考えていく．特に，Crouzeix--Raviart 要素や Raviart--Thomas 要素を曲単体上でどう定義するかは，今後の大きなテーマになる．

まとめ

この第1回では，曲線領域を有限要素法で扱うときの基本的な問題意識を整理した．通常の polygonal mesh では，円領域は多角形として近似される．一方，exact-curved FEM の立場では，曲線境界を曲線として扱い，物理領域により忠実な有限要素空間を構成することを目指す．

ただし，曲線幾何を高次にすることと，未知関数の近似次数を高次にすることは同じではない．この区別が，今回の数値実験を理解する鍵である．

次回は，要素写像 \[ F_K=\Psi_K\circ\Phi_{T_K} \] の意味をもう少し丁寧に見ていく．特に，affine core と curved correction を分けることで，何が見通しやすくなるのかを説明する．

参考文献

このシリーズの理論的背景として，曲線要素上の有限要素補間，曲線境界をもつ領域での有限要素法，異方性補間誤差評価，および筆者の関連研究を挙げておく．第1期では，特に Bernardi の exact curved element の考え方と，Ciarlet--Raviart，Zlámal，Scott，Lenoir らの曲線要素に関する古典的研究を意識している．また，本シリーズは，筆者の異方性メッシュ上の補間誤差解析，新しい幾何パラメータ，および投稿中の exact-curved Lagrange FEM 論文へ接続する研究ノートでもある．

  1. C. Bernardi, Optimal finite-element interpolation on curved domains, SIAM Journal on Numerical Analysis, 26(5), 1212--1240, 1989.
  2. 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.
  3. M. Zlámal, Curved elements in the finite element method. I, SIAM Journal on Numerical Analysis, 10(1), 229--240, 1973.
  4. R. Scott, Finite Element Techniques for Curved Boundaries, Ph.D. thesis, Massachusetts Institute of Technology, Cambridge, MA, 1973.
  5. M. Lenoir, Optimal isoparametric finite elements and error estimates for domains involving curved boundaries, SIAM Journal on Numerical Analysis, 23(3), 562--580, 1986.
  6. C. M. Elliott and T. Ranner, Finite element analysis for a coupled bulk--surface partial differential equation, IMA Journal of Numerical Analysis, 33(2), 377--402, 2013.
  7. G. Dziuk and C. M. Elliott, Finite element methods for surface PDEs, Acta Numerica, 22, 289--396, 2013.
  8. 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.
  9. H. Ishizaka, Exact-curved Lagrange finite elements for the Poisson problem in two dimensions, submitted.

ページの先頭へ