Efficient finite element quadrature on immersed voxel data

技術要約：埋め込みボクセルデータにおける効率的な有限要素求積法

問題提起
高解像度のコンピュータ断層撮影（CT）スキャンは、形状と材料特性に関する詳細なボクセル単位の情報を提供する。すべてのボクセルを有限要素として扱う直接的な数値シミュレーションも可能であるが、これは膨大な数の未知数（数十億個）を生じさせることが多い。自由度を減らすためにメッシュを粗視化する場合、正確な数値求積が大きな課題となる。単一の有限要素が複数のボクセル（CTのグレー値による不均質な材料特性など）を包含する場合、剛性行列の被積分関数は不連続になる。標準的なガウス・ルジャンドル求積法は、滑らかな多項式被積分関数を前提としているため、このような文脈では大きな積分誤差が生じる。

標準的な解決策である「合成求積（composed integration）」は、要素の積分をボクセル境界で細分化し、実質的にすべてのボクセル上で積分を行うものである。これは正確ではあるが、膨大な数の積分点が必要となるため計算コストが高い。これに代わる「事前積分（pre-integration）」は、ボクセル固有のサブ剛性行列を一度に計算し、後に組み立てる手法である。しかし、事前積分には決定的な限界がある。要素あたりのボクセル数が変化する場合や、多項式次数が変化する場合、あるいは最も重要な点として、適応的メッシュ細分化（adaptive mesh refinement）が行われる場合に、容易に適応できない。適応的なスキーム（マルチレベル $hp$ 法など）では、基底関数が細分化レベルによって変化するため、事前計算されたサブ行列が無効になってしまう。

手法
著者らは、これらの制限に対処するために、**ボクセルベースのモーメント適合（voxel-based moment fitting）**スキームを提案している。この手法は、モーメント適合法を埋め込みボクセル幾何学へと拡張したものである。核心となるアイデアは、各有限要素に対してカスタムの積分則を導出し、不連続な材料分布を積分点そのものではなく、積分重みに直接吸収させることである。

1. 数学的定式化: 要素剛性積分を、不連続な材料関数 $\alpha(x)$ と滑らかな部分 $h(x)$ に分割する。積を直接積分する代わりに、この手法は、固定された積分点 $x_i$ に対するカスタム重み $w_j$ を決定するための線形方程式系を解く。
2. ラグランジュ基底: 積分点で定義されたラグランジュ基底関数を選択することで、モーメント方程式の系は対角化される。重みはモーメントとして計算される：$w_j = \int_{\Omega_e} \alpha(x) l_j(x) d\Omega$。
3. ボクセル積分: モーメント方程式の右辺は、ボクセルレベルでの合成求積を用いて評価される。決定的な点として、著者らは事前積分スキーム（アルゴリズム2）を導入しており、そこでは幾何学的なモーメント（ボクセル上の $\int l_j(x) dx$）が行列 $S$ に事前計算される。最終的な要素の重みは、要素固有の材料値 $\alpha_k$ と $S$ との単純な縮退（contraction）によって得られる。これにより、各要素に対してボクセル上の積分を再評価する手間を回避できる。
4. 適応性: 重みは材料分布と現在の要素の特定の基底関数に基づいて計算されるため、本手法はグローバルなメッシュトポロジーに依存しない。これにより、マルチレベル $hp$ 法のような、細分化レベルによって基底関数が変化する手法ともシームレスに機能する。

主な貢献

• 基底に依存しない求積: 本手法は、任意の多項式次数および任意のボクセルレイアウトを持つ要素に対して、剛性行列の厳密な積分を可能にする。これは、事前積分技術とは異なり、事前計算されたサブ行列を必要としない。これにより、これまでボクセルベースの事前積分と互換性がなかった適応的細分化戦略（マルチレベル $hp$ など）が可能になる。
• 計算効率: 本スキームは、ボクセルごとの材料情報をカスタムの求積重みに吸収する。これにより、グローバル剛性行列の組み立て時における、全ボクセルにわたる高コストなループというボトルネックを排除できる。
• 事前積分モーメント行列: 事前計算された行列 $S$ の導入により、重みを決定するための計算コストが大幅に削減され、材料の解像度と有限要素基底がデカップリング（分離）される。

結果
論文では、以下の3つの数値例を通じて本手法を検証している。

1. スケーリング解析（不均質な単位立方体）: $768^3$ 個のボクセルを持つ立方体を用いた、合成求積、事前積分、およびボクセルモーメント適合の比較。
   • 精度: ボクセルモーメント適合は、合成求積と同一のひずみエネルギーを与え、数学的な等価性を確認している。
   • 性能: 要素あたりのボクセル数が（$p=1$ において）辺の長さに対して $\approx 6$ 個を超える場合、ボクセルモーメント適合は事前積分よりも優れている。高いボクセル数（例：辺あたり384ボクセル）において、ボクセルモーメント適合は、合成求積と比較して99.6%、事前積分と比較して**89.7%**の実行時間の短縮を達成した。
   • 損益分岐点: 効率のクロスオーバーは多項式次数 $p$ に依存する。$p$ が高くなると、損益分岐点はより多くのボクセルを含む要素へとシフトする。
2. 適応的L字型領域（2D）: マルチレベル $hp$ 適応化を用いた、再入角の特異点を持つ平面応力問題への適用。
   • 本スキームは、特異点およびボクセル化された境界の幾何学的な階段状構造を正常に解決する。
   • この研究では、離散化誤差と幾何学的誤差の相互作用を定量化しており、要素がボクセルサイズ以下に細分化された後は、さらなる細分化は物理的な領域ではなくボクセル化された幾何学を解決するだけになるため、適応的細分化を停止すべきであることを示している。
3. 相場破壊（Phase-Field Fracture, 3D）: 穴のある切り欠きプレートと空間的に変化するヤング率を用いたシミュレーション。
   • 本手法は、力学的平衡と相場発展の結合を扱う。劣化関数 $g(s)$ が被積分関数の次数を上昇させる（これにより剛性のアンダーインテグレーションを招く）が、結果として限界荷重や亀裂経路はほとんど影響を受けないことが示された。
   • 実現可能性: ボクセルモーメント適合を用いた50ステップのフルシミュレーションには30時間を要した。著者らは、合成求積を使用した場合、この問題は計算上不可能になる（118日を要すると推定）と述べている。

意義
本論文は、ボクセルベースのモーメント適合が、CT由来のデータに対する高忠実度シミュレーションのための堅牢かつ効率的な経路を提供すると主張している。その主な意義は、材料の解像度を有限要素基底からデカップリングすることにあり、これにより、従来の事前積分技術では不可能であった適応的細分化戦略を可能にしている。求積の計算コストを合成積分と比較して数桁削減することで、本手法は大規模な3次元不均質破壊シミュレーションを標準的なハードウェアで実行可能なものにしている。著者らは、本手法は重みの計算以外に標準的な有限要素コードを変更する必要がなく、既存のソルバーへの実用的な追加要素であると強調している。