# 2次元Poisson方程式に対するIGA

$%2018/12/06
\newcommand\setN[0]{\mathbb{N}}
\newcommand\setZ[0]{\mathbb{Z}}
\newcommand\setQ[0]{\mathbb{Q}}
\newcommand\setR[0]{\mathbb{R}}
\newcommand\setC[0]{\mathbb{C}}
\newcommand\setH[0]{\mathbb{H}}
\newcommand\setE[0]{\mathbb{E}}
\newcommand\pare[1]{{(#1)}}
\newcommand\Pare[1]{\left(#1\right)}
\newcommand\curl[1]{\{#1\}}
\newcommand\Curl[1]{\left\{#1\right\}}
\newcommand\squa[1]{[#1]}
\newcommand\Squa[1]{\left[#1\right]}
\newcommand\abs[1]{\lvert#1\rvert}
\newcommand\Abs[1]{\left\lvert#1\right\rvert}
\newcommand\norm[1]{\|#1\|}
\newcommand\Norm[1]{\left\|#1\right\|}
\newcommand\floor[1]{\lfloor#1\rfloor}
\newcommand\Floor[1]{\left\lfloor#1\right\rfloor}
\newcommand\ceil[1]{\lceil#1\rceil}
\newcommand\Ceil[1]{\left\lceil#1\right\rceil}
\newcommand\angl[1]{\langle#1\rangle}
\newcommand\Angl[1]{\left\langle#1\right\rangle}
\newcommand\transpose[1]{\,{\vphantom{#1}}^t\!#1}
\newcommand\sfrac[2]{#1/#2}
\newcommand\od[2]{\frac{d#1}{d#2}}
\newcommand\pd[2]{\frac{\partial#1}{\partial#2}}
\newcommand\sod[2]{\sfrac{d#1}{d#2}}
\newcommand\spd[2]{\sfrac{\partial#1}{\partial#2}}
\newcommand\tensor[2]{\boldsymbol{\mathrm{T}}^#1_#2}
\DeclareMathOperator{\tr}{tr}
\DeclareMathOperator{\Re}{Re}
\DeclareMathOperator{\Im}{Im}
\DeclareMathOperator{\supp}{supp}
\DeclareMathOperator{\sgn}{sgn}
\DeclareMathOperator{\Diff}{Diff}
\DeclareMathOperator{\Sym}{Sym}
\DeclareMathOperator{\span}{span}
\DeclareMathOperator{\grad}{grad}
\DeclareMathOperator{\rot}{rot}
\DeclareMathOperator{\div}{div}$

これは[数値計算Advent Calendar 2018](https://qiita.com/advent-calendar/2018/numerical_analysis)の12月12日の記事ですが, 公開日は2018年12月26日です.
遅くなってしまい大変申し訳ありません.

みなさん, [IGA(Isogeometric Analysis)](https://en.wikipedia.org/wiki/Isogeometric_analysis)はご存知でしょうか.
IGAは日本語の情報が少なく, たぶん知らない方も多いなので, 今回の記事ではざっくりと解説することを目標にします.

## tl;dr
ざっくり, 本当にざっくり言えばIGAは「数値計算・ミーツ・形状表現」です.
ここでの数値計算とはFEM(Finite Element Method, 有限要素法)でも使われているGalerkin法のことで, 形状表現とはNURBS(Non-uniform Rational B-spline)のことです.

通常のFEMでは, 近似解を構成するために解析対象の形状をメッシュ分割します.
それに対してIGAでは三角形分割のようなメッシュ分割は行いません.
これを可能にしているのがNURBSで, 近似解構成のための基底関数が形状表現から自然に誘導される事を利用しています.

IGAで解ける問題は幾つかありますが, 本記事では2次元Poisson方程式を例として解説します.

<!-- IGAは基本的にはFEMと同じで, Galerkin法を使って偏微分方程式の数値解を得る手法です.
通常のFEMと大きく異なるのは, 物体形状を単体分割すること無しに解を離散化する手法にあります.
具体的には, この手法ではNURBSから誘導して基底関数を構成することで離散化してます.

つまりIGAは, NURBSで表された形状の上の偏微分方程式の近似解を, Galerkin法で求める手法です.
-->

<!-- ## はじめに

## NURBS
NURBSはNon-uniform Raitonal B-splineの略です.
詳しくは次のpdfを参照して下さい.
https://drive.google.com/file/d/1JK2tBEXQavGOq2nlkpJyoeQC0QoIUOQ3/view
-->

## 2次元Poisson方程式
領域$\Omega$がNURBS多様体で表されているとします.
次の2次元Poisson方程式を考えます.

\begin{align}
\pd{^2u}{x^2}+\pd{^2u}{y^2}+f&=0 \\
u|_{\Gamma_D}&=g \\
\nabla u\cdot \boldsymbol{n}|_{\Gamma_N}&=0
\end{align}

ここで, $f$は与えられた$\Omega$上の関数で, $u$が求めるべき未知関数です.
$\Gamma_D$はDirichlet境界, $\Gamma_N$はNeumann境界で$\overline{\Gamma_D\cup\Gamma_N}=\partial \Omega, \Gamma_D\cap\Gamma_N=\varnothing$を満たします.


とくに$f=0$であればLaplace方程式と呼ばれ, その解$u$は調和関数と呼ばれます.

定常の熱伝導方程式と思えば, $u$が温度, $f$が$\Omega$上の熱源, $g$が境界$\Gamma_D$上の温度, $\nabla u\cdot \boldsymbol{n}|_{\Gamma_N}=0$は$\Gamma_N$上での断熱条件に相当します.

### 数値解析の方針
IGAの手順を書けば次のようになります:

1. 偏微分方程式を弱形式に書き直す
1. 領域をNURBSで表す
1. 解の精度のためにNURBS多様体の細分(refinement)を取る
1. NURBSから誘導される基底関数で弱形式の解空間を離散化して近似解を求める(Galerkin法)

通常のFEMとは基底関数の作り方が違うだけなので, 計算手順も基本的に同じになります.

### 1. 弱形式の導出
Galerkin法のために弱形式に書き換えます.
$\Gamma_D$上で$0$を取るテスト関数$w$を掛けて$\Omega$で積分します.

\begin{align}
0
&=\int_\Omega w\Pare{\pd{^2u}{{x_1}^2}+\pd{^2u}{{x_2}^2}+f} d{x_1}d{x_2} \\
&=\int_\Omega \Pare{\pd{}{{x_1}}\Pare{w\pd{u}{{x_1}}}+\pd{}{{x_2}}\Pare{w\pd{u}{{x_2}}}}d{x_1}d{x_2} \\
&\qquad -\int_\Omega\Pare{\pd{w}{{x_1}}\pd{u}{{x_1}}+\pd{w}{{x_2}}\pd{u}{{x_2}}}d{x_1}d{x_2} \\
&\qquad +\int_\Omega wf \ d{x_1}d{x_2} \\
&=\int_{\partial \Omega} \Pare{{n_1}\Pare{\pd{u}{{x_1}}}+{n_2}\Pare{\pd{u}{{x_2}}}}w \ ds \\
&\qquad -\int_\Omega\Pare{\pd{w}{{x_1}}\pd{u}{{x_1}}+\pd{w}{{x_2}}\pd{u}{{x_2}}}d{x_1}d{x_2} \\
&\qquad +\int_\Omega wf \ d{x_1}d{x_2} \\
&=-\angl{w,u}+\squa{w,f} \\
\end{align}

ここで$\boldsymbol{n}=(n_1,n_2)$は外向き単位法線ベクトルで, 記号$\angl{\cdot, \cdot}, \squa{\cdot, \cdot}$は次のように定めました.

\begin{align}
\angl{u,v}&=\int_\Omega\Pare{\pd{v}{{x_1}}\pd{u}{{x_1}}+\pd{v}{{x_2}}\pd{u}{{x_2}}}d{x_1}d{x_2} \\
\squa{u,v}&=\int_\Omega uv \ dx_1dx_2
\end{align}

つまり, 任意のテスト関数$w$に対して

\begin{align}
\angl{w,u}=\squa{w,f}
\end{align}

が充たされるような$u$を求めれば良いということになります.
これを弱形式と言います.
これに対して元の偏微分方程式は強形式と呼ばれ, 弱形式と(数学的な厳密さを言わなければ)等価です.

### 2. NURBSによる領域の形状表現
NURBSの詳細は[NURBS多様体による形状表現pdf](https://drive.google.com/file/d/1JK2tBEXQavGOq2nlkpJyoeQC0QoIUOQ3)に書いたのでそちらを参照してください.
今回使う性質で重要なものは次の2点です.

* 矩形領域$D=[a_1, b_1]\times[a_2,b_2]\subset\setR^2$から$\setR^2$への写像$\boldsymbol{p}$で, $\Omega$をパラメータ付けます.
* そのパラメータ付けは$n$個の基底関数$B_i:D\to\setR$と制御点$\boldsymbol{a}_i\in\setR^2$の線型結合で表されます. $\boldsymbol{p}(t_1,t_2)=\sum_{i_1,i_2} B_{i_1i_2}(t_1, t_2)\boldsymbol{a}_{i_1i_2}$といった具合です.

領域$\Omega$の具体例として次の図のような円環領域の半分を考えましょう.


この形状$\Omega$は次のようにして作れます.

* ノット列 $\boldsymbol{k}_1=(0,0,1,1), \boldsymbol{k}_2=(0,0,0,1,1,2,2,2)$
* 次数 $p_1=1, p_2=2$
* 重み $w_{i_1i_2}$ (図の通り)
* 制御点 $\boldsymbol{a}_{i_1i_2}$ (図の通り)


図には書いてませんが, 内半径1, 外半径2としています.

$D$は矩形領域で$D=[0,1]\times[0,2]$です.

### 3. NURBS多様体の細分
IGAにおける基底関数という用語には次の2つの意味があるため注意が必要です.

* NURBSの形状表現のための基底関数 $B_{i_1i_2}:D\to\setR$
* Galerkin法で使う基底関数 $N_i:\Omega\to\setR$

Galerkin法で使う基底関数$N_i$はNURBS基底関数$B_{i_1i_2}$から誘導されるため実質的にこれらは同じで, 混乱の恐れはありません.
($B_{i_1i_2}$は$\boldsymbol{p}(t_1,t_2)=\sum_{i_1,i_2} B_{i_1i_2}(t_1, t_2)\boldsymbol{a}_{i_1i_2}$のために2つの添字が走りますが, Galerkin法での基底関数は添字一つで十分なので束ねて$N_i$と書いています)


ただしNURBSによる形状表現が決まった段階で基底関数(と制御点)の数は決まってしまい, このままでは近似解を得るためには不十分という問題があります.

つまり近似解の精度向上には基底関数を増やす必要があり, ここで細分(refinement)と呼ばれる操作が役立ちます.
細分はパラメータ付け$\boldsymbol{p}:D\to\Omega$を変えることなく, 基底関数と制御点の数を増やす操作のことです.
(こちらも詳しくは[NURBS多様体による形状表現pdf](https://drive.google.com/file/d/1JK2tBEXQavGOq2nlkpJyoeQC0QoIUOQ3)を参照して下さい)
$\Omega$に対して細分を行って制御点を増やした例が次の図です.



この細分の操作は, 有限要素法で言えばメッシュを細かくする操作に相当します.

### 4. Galerkin法による近似解構成
理論的背景は他の文献を当ってもらうとして, ここでは解法のみ簡単に述べます.

$w$をテスト関数とした弱形式は次のように表されるのでした.

\begin{align}
\angl{w,u}=\squa{w,f} \\
\end{align}

NURBS基底関数$N_i \ (i\in I=\{1,\dots,n\})$を用いて, 厳密解$u$を$\tilde{u}$で近似します:

\begin{align}
\tilde{u}
&=\sum_{i\in I}[u]_iN_i=\tilde{v}+\tilde{g} \\
\tilde{v}
&=\sum_{i\in I_D}[v]_iN_i \\
\tilde{g}
&=\sum_{i\in I}[g]_iN_i \\
\end{align}

ここで$[u]_i, [v]_i, [g]_i\in \setR$は係数です.
$I_D$はDirichlet条件から決まる添字の集合で, 次を満たします.

\begin{align}
i\in I_D\Rightarrow N_i|_{\Gamma_D}=0
\end{align}

$\tilde{g}$は$g$を$\Gamma_D$上で近似できるように係数$[g]_i$を決定して構成されます. ($[g]_i$の決定は難しくなく, とくに今回の例では定数のみを考えるのでより簡単です)
$[u]_i=[v]_i+[g]_i$なので, あとは$[v]_i$を求めればよいことになります.

$w$も同様に$\tilde{w}=N_j \ (j \in I_D)$のようにして離散化します.
これらを先の弱形式に代入すれば, $[v]_i$を未知変数とする連立方程式($j\in I_D$)に帰着されます.

\begin{align}
\sum_{i\in I_D}[v]_i\angl{N_j,N_i}+\sum_{i\in I}[g]_i\angl{N_j,N_i}=\squa{N_j,f} \\
\end{align}

積分はNURBS多様体で入る座標$(t_1, t_2)$で行うため, 次のようにして変数を取り替えます.

\begin{align}
\angl{N_j,N_i}
&=\int_\Omega\Pare{\pd{N_j}{{x_1}}\pd{N_i}{{x_1}}+\pd{N_j}{{x_2}}\pd{N_i}{{x_2}}}d{x_1}d{x_2} \\
&=\int_D\Pare{\sum_{m,k,l}\pd{t_k}{x_m}\pd{t_l}{x_m}\pd{N_j}{{t_k}}\pd{N_i}{{t_l}}}\det_{k,l}\Pare{\pd{x_k}{t_l}}d{t_1}d{t_2} \\
\squa{N_j,f}
&=\int_\Omega N_j f \ dx_1dx_2 \\
&=\int_D N_j f \det_{k,l}\Pare{\pd{x_k}{t_l}}d{t_1}d{t_2}
\end{align}

### 計算例 A-1
簡単のために$f=0$とし, 次のように$\Gamma_D, g$を定めます.


この条件で数値計算した結果が次のグラフです.


厳密解$u$が調和関数になることを思い出せば, これは$\log(z)$の実部と同じ形になることが分かります.
つまり$u$は対数関数のグラフを原点中心に回転させたものに一致し, 次の図はその断面です.


### 計算例 A-2
簡単のために$f=0$とし, 次のように$\Gamma_D, g$を定めます.


この条件で数値計算した結果が次のグラフです.


厳密解$u$が調和関数になることを思い出せば, これは$\log(z)$の虚部と同じ形([常螺旋面](https://en.wikipedia.org/wiki/Helicoid))になることが分かります.

### 計算例 A-3
全てDirichlet条件として$f=4$とすれば次のようになります.


この条件で解いた結果が次のグラフです.


この計算結果では, 中央付近に不自然な山があるように見えます.
これは基底関数の数が十分でないために起こっており, もう少し細かい細分を取ると改善できます.


新しい基底関数で計算した結果が次の図です.
不自然な山は見える範囲では現れず, 近似解の精度が上がっている事が確認できます.


これまでの例と異なり, この条件に対する厳密解$u$を陽に書き下すのは難しそうです.
前述のようにこれも物理的解釈が出来て, 「$\Omega$全体に熱源$f$があって境界での温度が$g$で固定された定常状態」というものです.

### 計算例 B-1
$f=4$に対しても厳密解($u=1-{x_1}^2-{x_2}^2$)が得られるような形状(単位円板)・境界条件として次を与えましょう.


この領域$\Omega$のNURBSによる形状表現は次のようになります.

* ノット列 $\boldsymbol{k}_1=(0,0,0,1,1,1), \boldsymbol{k}_2=(0,0,0,1,1,1)$
* 次数 $p_1=2, p_2=2$
* 重み $w_{i_1i_2}$ (図の通り)
* 制御点 $\boldsymbol{a}_{i_1i_2}$ (図の通り)


基底関数$N_i$は次のようになります.


この条件で解いた結果が次のグラフです.


厳密解$u=1-{x_1}^2-{x_2}^2$と近似解$\tilde{u}$の差を取れば次のようになります.
ただし差を大きく見せるため, 高さ方向に200倍して表示しています


細分を取って基底関数の数を増やすと次のようになります.


これらの基底関数を使って近似解$\tilde{u}$を構成すると次のようになりました.


$\tilde{u}$の見た目だけでは細分前後の変化が分かり難いですが, 同じように差$u-\tilde{u}$を取ると誤差が減っていることが確認できます.


## まとめ
[isoは「等しい」](https://ja.wiktionary.org/wiki/iso-), geometricは「幾何的に」という意味なので, IGA(isogeometric analysis)は直訳すれば「幾何的に等しい解析」ということになります.

通常の有限要素法では, 領域$\Omega$を三角形分割する必要があるため, その意味で形状に対する近似が入ります.
これに対してIGAでは, 今回見たように$\Omega$がNURBSで表現されていれば形状に対する近似無しに数値解を得ることが出来ます.
これが「幾何的に等しい解析」の意味する所です.

CADなどで使われる形状の多くはNURBSで表現可能であり, その意味でNURBSは非常に柔軟に様々な形状を表すことが可能です.
IGAではNURBSで表現可能な形状しか扱えませんが, これはデメリットでは無く, 寧ろCADで作成した形状からシームレスに数値計算を行えるためにメリットと考える方が適切です.

## コード
今回の数値計算のためのコードはJuliaで書きました.
全然綺麗には書けてませんが, 一応githubに上げておきます.
https://github.com/hyrodium/IGA-PoissonEquation

## 参考文献
数学的にはGalerkin法を使っているため, 有限要素法に関する文献も有用です.

* [有限要素法概説](http://www.saiensu.co.jp/?page=book_details&ISBN=ISBN4-7819-0911-6)
Poisson方程式の式変形とかはこの書籍を参照しました.
* [Isogeometric Analysis: Toward Integration of CAD and FEA](https://www.wiley.com/en-sg/Isogeometric+Analysis%3A+Toward+Integration+of+CAD+and+FEA-p-9780470748732)
IGAに関してまず参照される書籍がこちらです.
* [Geometric Modeling with Splines An Introduction](https://www.crcpress.com/Geometric-Modeling-with-Splines-An-Introduction/Cohen-Riesenfeld-Elber/p/book/9781568811376)
NURBSに関する書籍ではこれが良いです.
* [NURBS多様体による形状表現](https://drive.google.com/file/d/1JK2tBEXQavGOq2nlkpJyoeQC0QoIUOQ3)
私が書いたNURBS入門です.

## おわりに
QiitaではなくHackMDで書いたため, コメントなどが出来ません.
質問・間違いなどありましたら[筆者のTwitter](https://twitter.com/Hyrodium)までご連絡ください.