CAE-AI Lab 近畿大学 English

3Dで学ぶ電磁気学

G-8 高次基底と適応解析

G-7では、大域的に許容なcurl側とdiv側の解を使い、真値を知らずに誤差を挟む方法を学びました。最終回では、最低次基底を適合性の土台として保ちながら、滑らかな領域では高次基底を加え、特異点や材料界面ではメッシュを細分する考え方を学びます。

NGSolveの高次要素は、最低次空間を残したまま、辺・面・要素内部の基底を階層的に追加します。特異点や材料界面にはh細分、滑らかな領域にはp次数上昇を使い、必要なら両方を組み合わせます。この回では、基底関数を一つずつ見てから、誤差情報を SOLVE・ESTIMATE・MARK・REFINE へ渡します。

研究室所蔵の亀有資料「Electromagnetic Analyses Using Higher Order Hierarchic Finite Elements」と「Adaptive Meshing with Bossavit's Error Estimator」にあるNGSolve実装を照合しています。ここで扱う局所回復型指標は適応位置を選ぶための実用量であり、大域的な許容性を確認しない限り、G-7の保証付き上下界とは呼びません。

まず30秒でつかむ

有限要素解を良くする方法は、要素を細かくするh細分だけではありません。解が滑らかな場所では、同じ要素へ高次の辺・面・内部モードを足すp次数上昇が効きます。一方、角の特異点、薄いギャップ、材料界面では局所的なh細分が必要です。最低次基底は捨てず、接続と適合性を保つ土台として残します。

\[ V_h^k(1)\subset V_h^k(2)\subset\cdots\subset V_h^k(p), \qquad \eta^2=\sum_{K\in\mathcal T_h}\eta_K^2, \qquad \sum_{K\in\mathcal M}\eta_K^2\ge\theta\eta^2 \]

最初の包含関係が階層基底、中央が要素ごとの誤差情報、最後がDorfler markingです。可視化では、単にメッシュが細かくなる様子ではなく、どの自由度を足し、どの要素を選び、目的量がどう改善したかを数値と一緒に追います。

式を読む順番

1. 適合空間を選ぶ。 スカラーなら H1、循環をもつ場なら HCurl、法線磁束を保つ場なら HDiv を選びます。

2. 自由度の所属を読む。 最低次の頂点・辺・面自由度に、高次の辺・面・要素内部モードがどう加わるかを確認します。

3. 誤差情報を要素へ分ける。 全体誤差 $\eta$ を要素寄与 $\eta_K$ へ分解し、局所回復型か保証付き相補型かを明記します。

4. MARKで候補を選ぶ。 大きな寄与を持つ要素集合 $\mathcal M$ を、Dorfler条件で機械的に選びます。

5. 滑らかさからhかpを決める。 特異・不連続ならh、滑らかならp、判断が難しければhp候補を比較します。

6. 同じ診断量で再計算を評価する。 SOLVE・ESTIMATE・MARK・REFINE を繰り返し、エネルギー誤差と目的量誤差を混同しません。

この回の到達点

1. 最低次基底と高次の辺・面・要素内部基底の役割を区別できる。

2. NGSolveの H1、HCurl、HDiv で、order と自由度の所属を確認できる。

3. curl側とdiv側の場を別の適合空間へ局所回復し、要素指標を作れる。

4. h細分とp次数上昇を、解の滑らかさに応じて選べる。

5. エネルギー誤差、局所平均で定義した目的量誤差、要素品質を別の診断量として監視できる。

修士課程の電磁気学演習

自由度を増やすなら、hとpのどちらへ使うか

電磁気の問い: 同じ自由度予算で、滑らかな磁界と角特異性を持つ磁界を改善する。要素を細分するh法と、辺・面・内部へ高次モードを足すp法を、何を根拠に選べばよいだろうか。

観察ポイント: 滑らかな解と特異な解について、p次数を上げた場合と局所h細分した場合の誤差低下を順位付けする。エネルギー誤差と目的量誤差で選択が変わる場合も考える。

操作と観測: 可視化G-8-4で解の滑らかさを切り替え、可視化G-8-5でp次数・自由度・誤差を同時に比較する。可視化G-8-3では同じ誤差指標を使ってh細分とp高次化を競わせる。

判定基準と微分幾何による説明

高次化しても離散de Rham複体と接線・法線の適合性は保ち、階層空間は $V_h^k(1)\subset V_h^k(2)\subset\cdots\subset V_h^k(p)$ と広がる。滑らかな領域では階層的pモードが効率よく誤差を落とし、角、薄いギャップ、材料界面の低正則性には局所h細分が必要になる。

診断: エネルギー誤差、局所目的量、最小要素品質、$\boldsymbol C\boldsymbol G=\boldsymbol0$ を別々に監視する。pを上げても誤差が停滞するなら特異性または幾何近似を疑い、h細分で品質が悪化するなら細分規則を見直す。

研究へ: NGSolveで `HCurl`・`HDiv` の階層自由度を確認し、`SOLVE`・`ESTIMATE`・`MARK`・`REFINE` にh/p選択と目的量監視を組み込む。

階層基底で近似空間を広げる

最低次要素は、向き、接続、境界条件、大域循環を堅牢に表す土台です。したがって「高次要素へ替える」とは最低次基底を消すことではありません。階層基底では

\[ V_h^k(1)\subset V_h^k(2)\subset\cdots\subset V_h^k(p) \]

となり、$p=1$ の空間を保ったまま新しいモードを足します。滑らかな解に対する代表的な近似評価

\[ \|u-u_h\|_{H^1(\Omega)} \le C h^p |u|_{H^{p+1}(\Omega)} \]

を見ると、同じ $h$ でも $p$ を上げる意味が分かります。反対に、角の特異点や材料係数の不連続では高階微分が滑らかでないため、p次数だけを上げても効率が落ちます。この違いが、h細分とp高次化を使い分ける理由です。

三角形上の完全多項式空間の次元は

\[ N_p=\dim\mathbb P_p(K) =\dfrac{(p+1)(p+2)}{2} \]

です。$p=1$ なら3自由度、$p=4$ なら15自由度です。自由度は増えますが、一つの要素内で滑らかな変化を表せるため、「粗い形状を低次要素だけで何段も割る」より少ない自由度で高精度になる場合があります。

NGSolveでは基底を階層的に足す

参照三角形

\[ K=\{(x,y):x\ge0,\ y\ge0,\ x+y\le1\}, \qquad \lambda_1=1-x-y,\quad\lambda_2=x,\quad\lambda_3=y \]

を考えます。$H^1$ の最低次頂点基底は $\phi_1=\lambda_1$ です。辺 $E_{12}$ に所属する代表的な高次バブルは、Legendre多項式 $P_n$ を使って

\[ \phi_{12}^{(p)} =\lambda_1\lambda_2 P_{p-2}(\lambda_2-\lambda_1), \qquad p\ge2 \]

と書けます。これは対象辺以外でゼロになり、既存の頂点基底を変えずに辺上の表現力を増やします。$H^1$ の階層分解は概念的に

\[ W_{hp}=W_{p=1} +\sum_E W_E+\sum_F W_F+\sum_K W_K \]

です。三次元の $H(\mathrm{curl})$ でも、最低次辺基底に高次の辺・面・要素内部基底を加えます。

\[ V_{hp}=V_{\mathrm{l.o.}} +\sum_E V_E+\sum_F V_F+\sum_K V_K \]

辺 $E_{12}$ の最低次Nedelec基底の一例は

\[ \vec N_{12} =\lambda_1\nabla\lambda_2-\lambda_2\nabla\lambda_1 \]

です。高次 $H^1$ 辺基底の勾配は同じ辺に所属する $H(\mathrm{curl})$ モードとなり、局所de Rham列を保ちます。

\[ \nabla W_E\subset V_E, \qquad \nabla W_F\subset V_F, \qquad \nabla W_K\subset V_K \]

三次元のde Rham列は

\[ H^1\xrightarrow{\nabla} H(\mathrm{curl})\xrightarrow{\nabla\times} H(\mathrm{div})\xrightarrow{\nabla\cdot}L^2 \]

と並びます。電磁気量とNGSolve空間の対応は次のように整理できます。

NGSolve空間主な電磁気量要素間で共有する量高次化で加わる基底
H1スカラーポテンシャル値辺・面・要素内部のバブル
HCurl磁気ベクトルポテンシャル、電界接線成分高次辺、面、要素内部モード
HDiv磁束密度、電流密度法線成分高次面、要素内部モード
L2体積密度、残差共有しない要素内部多項式

NGSolveでは、まず空間と次数を明示します。次の4行は空間の選び方を並べた例であり、同じ p を指定すれば自動的に一つの完全列になる、という意味ではありません。方程式とNGSolveの次数規約に合わせて組合せを決めます。

p = 3
fes_phi = H1(mesh, order=p)
fes_a = HCurl(mesh, order=p)
fes_b = HDiv(mesh, order=p, RT=True)
fes_rho = L2(mesh, order=max(p - 1, 0))

RT=True はRaviart--Thomas型の $H(\mathrm{div})$ 空間を選びます。GetDofNrs を使えば、特定の辺に所属する基底を一つだけ励起し、基底そのものとcurlを表示できます。

gfu = GridFunction(fes_a)
edge_dofs = fes_a.GetDofNrs(NodeId(EDGE, edge_nr))
gfu.vec[:] = 0
gfu.vec[edge_dofs[k]] = 1
Draw(gfu, mesh, "edge_basis")
Draw(curl(gfu), mesh, "curl_basis")

可視化 G-8-1 NGSolveの階層基底を一つずつ見る

最初のイメージ。 最低次基底を粗い輪郭、高次基底をその輪郭へ局所的な細部を足す筆だと考えます。次数を上げても最低次の辺・面自由度を捨てるのではなく、その上へ高次辺モードや要素内部モードを積み重ねます。どの基底が要素境界の情報を隣要素と共有し、どれが要素内部だけで消えるかを見ることが、適合性と静的縮約の理解につながります。

一要素の形だけでは「適合」の意味は見えません。共有辺 $e=K^-\cap K^+$ に共通の向きで単位接線 $\vec t$ と単位法線 $\vec n$ を固定すると、隣接要素間で一致させる量は

\[ u^-|_e=u^+|_e\quad(H^1),\qquad (\vec v^--\vec v^+)\cdot\vec t=0\quad(H(\mathrm{curl})),\qquad (\vec q^--\vec q^+)\cdot\vec n=0\quad(H(\mathrm{div})) \]

です。ベクトルの全成分が一致する必要はありません。$H(\mathrm{curl})$ では法線成分、$H(\mathrm{div})$ では接線成分が跳んでも、それぞれの適合性は壊れません。

Canvasの高次辺モードには、共有辺座標 $s\in[-1,1]$ 上の正規化代表トレース

\[ g_{p}(s)=(1-s^2)P_{p-2}(s)\qquad(p\ge2) \]

を使います。内部モードは $b_K=\lambda_1\lambda_2\lambda_3q_{p-3}$ とし、適合トレースは境界で零です。これは支持・トレース・階層包含を比較する診断モデルであり、NGSolve内部の正規化、符号規約、基底番号を再現するものではありません。

基底の種類とトレース
まず共有辺のどの量が左右で一致するか図を動かして確かめてください。

動かして確かめる順序。 初期状態では $p=1$ の最低次基底に固定されています。次数を動かすと $p=3$ の高次辺モードが開き、左要素の青実線と右要素の橙破線の適合トレースが共有辺上で重なります。HCurl と HDiv を切り替え、重なる成分が接線から法線へ変わる一方、非拘束成分は一致しなくてもよいことを確認してください。次に内部モードへ切り替え、境界トレースが零になることを確かめます。次数を上げても下段の $V_1$ が残り、追加モードだけが積み重なることも読み取ってください。

$H(\mathrm{curl})$ と二次元 $H(\mathrm{div})$ の対応には、90度回転 $\mathcal R(a,b)=(b,-a)$ を使っています。辺の単位接線を $\vec t$、単位法線を $\vec n=\mathcal R\vec t$ とすれば

\[ (\mathcal R\vec v)\cdot\vec n =\vec v\cdot\vec t \]

となり、接線トレースの連続性が法線トレースの連続性へ移ります。

高次化は無料ではありません。要素内部自由度を $i$、骨格上の自由度を $b$ と分けると、静的縮約では

\[ \widehat{\boldsymbol K}_{bb} =\boldsymbol K_{bb} -\boldsymbol K_{bi}\boldsymbol K_{ii}^{-1}\boldsymbol K_{ib} \]

を解き、内部自由度を要素ごとに消去できます。最低次ブロックを大域的な粗視化へ、高次バブルを局所平滑化へ分けられることも、階層基底を使う実装上の利点です。

局所補間で誤差の居場所を探す

スカラーポテンシャル側の有限要素解 $\vec H_h$ はcurl-freeですが、$\mu\vec H_h$ が選んだdiv適合空間へ入るとは限りません。そこで局所補間 $I_{\mathrm{div}}$ を使い、回復磁束密度を

\[ \vec B_h^{\,*}=I_{\mathrm{div}}(\mu\vec H_h) \]

と作ります。要素 $T$ の指標は

\[ \eta_{T,\mathrm{div}}^2 =\int_T\mu^{-1} \left|\vec B_h^{\,*}-\mu\vec H_h\right|^2\,dV \]

です。ベクトルポテンシャル側では、$\vec B_h=\nabla\times\vec A_h$ から回復磁界

\[ \vec H_h^{\,*}=I_{\mathrm{curl}}(\mu^{-1}\vec B_h) \]

を作り、

\[ \eta_{T,\mathrm{curl}}^2 =\int_T\mu \left|\vec H_h^{\,*}-\mu^{-1}\vec B_h\right|^2\,dV \]

を評価します。NGSolveの GridFunction.Set による局所射影と共有自由度の平均化は、この回復操作を実装する一つの方法です。ただし、回復場が境界条件、独立循環、総磁束まで満たすかは別に検査します。

SOLVE・ESTIMATE・MARK・REFINE

要素指標を全体へ集めると、計算領域全体の指標は

\[ \eta^2=\sum_{T\in\mathcal T_h}\eta_T^2 \]

です。単純な最大値基準では、最大指標に対する割合 $r$ を使って

\[ \mathcal M =\left\{T\in\mathcal T_h: \eta_T>r\max_{K\in\mathcal T_h}\eta_K\right\}, \qquad 0\lt r\lt1 \]

をMARKします。より安定したDörfler markingでは、指標の一定割合を担う最小集合を選びます。

\[ \sum_{T\in\mathcal M}\eta_T^2 \ge\theta \sum_{T\in\mathcal T_h}\eta_T^2, \qquad 0\lt\theta\lt1 \]

その後、MARKした要素だけを細分して再び解きます。

\[ \mathrm{SOLVE} \longrightarrow\mathrm{ESTIMATE} \longrightarrow\mathrm{MARK} \longrightarrow\mathrm{REFINE} \longrightarrow\mathrm{SOLVE} \]

可視化G-8-2では、材料界面の角近傍に大きな指標をもつ教材モデルを使います。細分段階、要素次数 $p$、MARK比率を変え、赤い要素へ自由度が集まる様子を確認してください。

可視化 G-8-2 局所指標に応じたh細分

要素数と推定誤差

h細分とp次数上昇

h細分は要素径 $h_T$ を小さくし、特異点や材料界面を局所的に解像します。一方、高次階層要素は既存の低次基底を残したまま次数を増やします。

\[ V_h^k(p)\subset V_h^k(p+1) \]

辺・面・要素内部の次数を独立に選ぶ場合は、隣接自由度の整合と局所de Rham列を壊さないことが前提です。

\[ p_e\le p_f\le p_K, \qquad V_h^0\xrightarrow{\boldsymbol G}V_h^1 \xrightarrow{\boldsymbol C}V_h^2 \xrightarrow{\boldsymbol D}V_h^3 \]
\[ \boldsymbol C\boldsymbol G=\boldsymbol0, \qquad \boldsymbol D\boldsymbol C=\boldsymbol0 \]

角、薄いギャップ、材料係数の不連続ではh細分が有効です。滑らかな領域ではp次数を上げる方が、同じ自由度数で速く収束することがあります。実務では滑らかさ指標 $s_T$ を併用し、$s_T$ が低い要素はh、高い要素はpへ振り分けます。

\[ s_T\lt s_{\mathrm{crit}}\Rightarrow h\text{-refinement}, \qquad s_T\ge s_{\mathrm{crit}}\Rightarrow p\text{-enrichment} \]

局所目的量とメッシュ品質を別に監視する

エネルギーノルム用の指標が小さくても、狭い領域の磁束密度や局所反力が高精度とは限りません。ただし $H(\mathrm{div})$ や $H(\mathrm{curl})$ の場では、点評価 $B_z(\vec x_0)$ は一般には連続な汎関数ではなく、厳密な誤差評価の目的量として常に定義できるわけではありません。小領域 $D_\varepsilon$ の平均やコイル積分など、有限要素空間上で有界な目的汎関数 $J_\varepsilon(u)$ を定めてgoal-oriented adaptivityを行います。

\[ J_\varepsilon(u) =\dfrac{1}{|D_\varepsilon|} \int_{D_\varepsilon}B_z(\vec x)\,dV, \qquad \left|J_\varepsilon(u)-J_\varepsilon(u_h)\right| \]

解が十分滑らかで点値が定義できる場合には、$D_\varepsilon$ を小さくした極限を点値の近似として監視できます。ただし、そのときも離散化誤差と平均化半径 $\varepsilon$ による誤差を分けて扱います。

随伴解 $z$ を使うと、残差が目的量へ与える寄与を要素ごとに分けられます。

\[ J_\varepsilon(u)-J_\varepsilon(u_h) \simeq \sum_{T\in\mathcal T_h}R_T(u_h)(z-z_h) \]

また、極端に小さい角をもつ要素では、細分しても条件数が悪化し、局所指標が減らないことがあります。三角形なら内接円半径 $r_{\mathrm{in},T}$ と外接円半径 $r_{\mathrm{circ},T}$ から

\[ q_T=\dfrac{2r_{\mathrm{in},T}}{r_{\mathrm{circ},T}}, \qquad 0\lt q_T\le1 \]

を監視できます。誤差指標、保存則、目的量、要素品質、反復解法残差を分けて記録することが、安定した適応解析に必要です。

保存則・検算とF/G各回の接続

メッシュが変わるたびに、接続行列の恒等式と相補場の挟み込みを再検査します。

\[ \boldsymbol C_h\boldsymbol G_h=\boldsymbol0, \qquad \boldsymbol D_h\boldsymbol C_h=\boldsymbol0 \quad\text{at every refinement level} \]
\[ R_{\mathrm{div},h}^{(m+1)} \ge R_{\mathrm{div},h}^{(m)}, \qquad R_{\mathrm{curl},h}^{(m+1)} \le R_{\mathrm{curl},h}^{(m)} \]

非入れ子空間、反復解法誤差、境界・位相制約の不整合があると単調性は崩れるため、期待ではなく数値で確認します。G-6の総合課題へ戻り、同じ電磁石を粗いメッシュから始め、保存則、上下界、目的量を監視しながら細分してください。

可視化 G-8-3 同じ節点数でh細分とp高次化を比べる

比較領域を一つの参照三角形 $K$ に揃え、面積座標 $\lambda_1+\lambda_2+\lambda_3=1$ 上の三次曲面

\[ u^*(\lambda_1,\lambda_2,\lambda_3) =27\lambda_1\lambda_2\lambda_3 \]

を共通の近似対象にします。比較レベルを $q$ とすると、どちらも同じ三角格子点

\[ N(q)=\dfrac{(q+1)(q+2)}{2} \]

を使います。ただし、h側は $q^2$ 個の小三角形ごとに一次、p側は一つの三角形全体で完全 $q$ 次です。

\[ V_h(q)=\{v:\ v|_T\in\mathbb P_1(T),\ T\in\mathcal T_q\}, \qquad V_p(q)=\mathbb P_q(K) \]

Canvasの灰色格子は固定した基準曲面 $u^*$、青と紫の面は二つの近似、赤い鉛直線は同じ標本点での差です。表示するRMSは固定した三角格子標本上の無次元診断値であり、$H^1$ 誤差評価や一般問題の収束保証そのものではありません。

まず、要素数と多項式次数のどちらが三次曲面の再現性を決めるか図を動かして確かめてください。
操作すると、同じ節点数の二つの近似を同時表示します。

動かして確かめる順序。 図を動かして $q=3$ の比較を開き、灰色の基準曲面と青いh-P1面、紫のp面の隙間を見ます。その後 $q=1\sim5$ を動かし、節点数が同じでも「小要素ごとの一次」と「一要素全体の高次」では作れる関数空間が異なることを、赤線とRMSで説明してください。

可視化 G-8-4 滑らかな解と特異な解でp高次化の効き方を比べる

要素内で十分滑らかな解なら、次数 $p$ の補間誤差は高次化によって速く減ります。代表的な評価は

\[ \|u-I_pu\|_{H^1(K)} \le C\,h_K^p|u|_{H^{p+1}(K)} \]

です。しかし $u(x)=\sqrt{x}$ のように端点で高階導関数が特異になると、右辺の滑らかさ仮定が満たされません。

\[ u'(x)=\dfrac{1}{2\sqrt{x}}, \qquad u'\notin L^\infty(0,1) \]

$e^x$ と $\sqrt{x}$ を切り替え、同じ次数 $p$ で補間誤差の減り方を比較してください。滑らかな領域ではp高次化、特異点近傍ではh細分という選択を、曲線の差として判断できます。

可視化 G-8-5 p次数・自由度・誤差を同時に見る

対象を単一四面体上のスカラー $H^1$ 完全 $p$ 次階層要素に固定します。次数を上げたとき、頂点・辺・面・要素内部の自由度がどのように積み重なり、滑らかな解を仮定した誤差モデルがどう低下するかを同じ $p$ 軸で比べます。

\[ \begin{aligned} N_p &=\underbrace{4}_{\text{頂点}} +\underbrace{6(p-1)}_{\text{辺}} +\underbrace{2(p-1)(p-2)}_{\text{面}} +\underbrace{\frac{(p-1)(p-2)(p-3)}{6}}_{\text{要素内部}}\\ &=\binom{p+3}{3}. \end{aligned} \]
\[ \frac{E_p}{E_1}=10^{-\alpha(p-1)},\qquad \alpha=0.60 \]

ここで $E_p/E_1$ は、要素内で解析的に滑らかな解を想定した比較用モデルです。実際の収束率は解・要素形状・ノルムに依存し、特異点や材料界面がある問題で指数低下を保証する式ではありません。

まず自由度の増え方と誤差の減り方を図を動かして確かめてください。
操作すると p=1 の固定基準と比較します。

動かして確かめる順序。 先に増加・低下の組合せを確認し、固定基準 $p=1$ と $p=5$ を比較します。その後スライダーを動かし、低次の4頂点自由度を残したまま辺・面・内部自由度が追加されること、自由度曲線は三次元では $O(p^3)$、滑らかな解のモデル誤差は片対数図で直線になることを説明してください。

理解の確認

1. 階層基底でp次数を上げても、最低次基底が消えないことを式で説明してください。

2. HCurl の辺自由度と HDiv の面自由度は、どのトレースを要素間で共有するか。

3. GetDofNrs(NodeId(EDGE, edge_nr)) で確認できるものは何か。

4. 局所回復型指標が、そのまま保証付き誤差上界ではない理由を説明してください。

5. 材料界面の角近傍では、p次数上昇だけよりh細分が有効になりやすいのはなぜか。

6. 磁気エネルギーが収束しても一点の $B_z$ が揺れる可能性があるのはなぜか。

7. 適応細分の各段階で $\boldsymbol C_h\boldsymbol G_h=\boldsymbol0$ を検査する意味を説明してください。

参考資料