[論文レビュー] Programming of Finite Element Methods in MATLAB

Long Chen|arXiv (Cornell University)|Apr 14, 2018

Matrix Theory and Algorithms参考文献 5被引用数 6

ひとこと要約

本稿では、MATLABを用いたポアソン方程式の線形有限要素法の実装に向けた体系的アプローチを提示する。スパース行列の効率的使用と三角形分割および境界条件のデータ構造に重点を置いている。スパース行列演算を用いて剛性行列と質量行列をアセンブルする方法、さまざまな境界条件の取り扱い方、数値積分則の適用を示し、ハイレベルなMATLABコーディングとスパース行列技術により最適な性能を達成している。

ABSTRACT

We discuss how to implement the linear finite element method for solving the Poisson equation. We begin with the data structure to represent the triangulation and boundary conditions, introduce the sparse matrix, and then discuss the assembling process. We pay special attention to an efficient programming style using sparse matrices in MATLAB.

研究の動機と目的

• ポアソン方程式を解くための、明快かつ効率的な線形有限要素法のMATLAB実装を開発すること。
• グローバルおよびローカルの頂点インデックスを用いて、三角形分割と境界条件を表現する最適なデータ構造を設計すること。
• MATLABにおけるスパース行列の使用により、有限要素アセンブリの計算効率を向上させ、メモリ使用量を削減すること。
• 単体要素における数値積分則を用いた剛性行列と質量行列のアセンブリに関する実用的ガイドラインを提供すること。
• ノイマン問題における適切な解の存在性の問題、特に核空間の取り扱いや適合条件の対処法を検討すること。

提案手法

• 三角形分割を2つの行列で表現する：幾何座標を表す 'e' と、単体のグローバル頂点インデックスを表す 'm' を用い、2次元では反時計回り、3次元では右ねじの法則に従って一貫した向きを保証する。
• 境界条件を密度行列 'lag' で符号化する。値は0（内部）、1（ディリクレ）、2（ノイマン）、3（ロビン）として、境界面をラベル付ける。
• MATLABの組み込みスパース行列形式（CSC）を活用し、大規模なスパースな剛性行列と質量行列を効率的に格納・操作する。
• 要素ごとにループを回し、スパース行列の加算を用いて局所的寄与を加算することで、グローバル行列をアセンブルする。密行列演算を避ける。
• 重心、中点、ガウス積分などの数値積分則を、単体の重心座標で適用し、要素行列を正確かつ効率的に計算する。
• ノイマン問題における剛性行列の核空間を扱うために、解を得た後で u = u - mean(u) を適用し、平均値をゼロに制約する。

実験結果

リサーチクエスチョン

• RQ1 MATLABにおけるスパースデータ構造を用いて、三角形分割と境界条件を効率的に表現する方法は何か？
• RQ2 MATLABのスパース行列演算を用いた有限要素行列のアセンブリにおいて、最適な方法は何か？
• RQ3 異なる数値積分則が、有限要素アセンブリの精度と効率に与える影響は何か？
• RQ4 ノイマン境界条件を伴うポアソン方程式の離散的設定における適切な解の存在性の問題は、どのように解決できるか？
• RQ5 ハイレベルなコーディングスタイルとローレベルなコーディングスタイルの違いが、性能と可読性に与える影響は何か？

主な発見

• MATLABにおけるスパース行列の使用は、特に大規模問題において、メモリ使用量を著しく削減し、行列アセンブリの速度を向上させる。
• 提案された 'e' と 'm' の2つの行列を用いたデータ構造により、三角形分割の局所的細分化および粗化が効率的に行える。
• 2次元では、重心則（N 個の評価）よりも台形則（NT ≈ 2N 個の評価）が関数評価回数が少ないため、より効率的である。
• 三角形では2次積分則（辺の中点を用いる）が存在するが、四面体では存在せず、代替の対称的則が必要となる。
• 離散的適合条件 ∫f dx = ∫g dS は数値的に強制されなければならない。そうでなければ、ノイマン問題は解をもたない可能性がある。
• ノイマン問題の解は定数分だけ一意でない。∫u dx = 0 を強制することで一意性が保たれ、L²誤差の計算が改善される。