いつになったらサンプルコードまでたどり着けるんだろうか...
境界の積分
ノイマン境界条件を処理するために境界でも積分しなければ足りません。
線分要素の局所座標系
この前と同じように局所座標系に変換していきます。元の要素上の点$\mathbf{x} = (x,y)$と局所座標系における値$\zeta$には以下のような関係が成り立ちます。
\mathbf{x} = \zeta (\mathbf{x_2} - \mathbf{x_1})
つまり、その点は線分の端から端までのうちどの辺りに存在するかを$[0,1]$で表したものとなります。
内挿関数
その局所座標系における内挿関数は以下のように表されます。
\begin{align}
\psi_1 &= 2(x-\frac{1}{2})(x-1) = 2x^2-3x+1\\
\psi_2&=-4x(x-1) = -4x^2+4x\\
\psi_3 &= 2x(x-\frac{1}{2}) = 2x^2-x
\end{align}
積分
\begin{align}
q_i^l \int_{\Gamma_{N_l}} \psi_i^l \psi_j^l \mathrm{d\Gamma}
&= q_i^l \int_0^1
\begin{pmatrix}4 \zeta^{4} - 12 \zeta^{3} + 13 \zeta^{2} - 6 \zeta + 1 & - 8 \zeta^{4} + 20 \zeta^{3} - 16 \zeta^{2} + 4 \zeta & 4 \zeta^{4} - 8 \zeta^{3} + 5 \zeta^{2} - \zeta\\- 8 \zeta^{4} + 20 \zeta^{3} - 16 \zeta^{2} + 4 \zeta & 16 \zeta^{4} - 32 \zeta^{3} + 16 \zeta^{2} & - 8 \zeta^{4} + 12 \zeta^{3} - 4 \zeta^{2}\\4 \zeta^{4} - 8 \zeta^{3} + 5 \zeta^{2} - \zeta & - 8 \zeta^{4} + 12 \zeta^{3} - 4 \zeta^{2} & 4 \zeta^{4} - 4 \zeta^{3} + \zeta^{2}\end{pmatrix} \mathrm{d}\zeta\\
&=\begin{pmatrix}
q_1 & q_2 & q_3
\end{pmatrix}
\begin{pmatrix}\frac{2}{15} & \frac{1}{15} & - \frac{1}{30}\\\frac{1}{15} & \frac{8}{15} & \frac{1}{15}\\- \frac{1}{30} & \frac{1}{15} & \frac{2}{15}\end{pmatrix}
\\
i &= 1,2,3
\end{align}
以後、これを$\mathbf{q^l}$と置きます。
全体の行列作成
もう一度要素に分割した際に導出した方程式を見てみましょう。
\begin{align}
\sum_{e=1} \left[
w_j^e u_i^e K_{i,j}^e
+w_j^e {f_j^e}
\right]
=
\sum_{l=1} w_j^e q_{j}^l\\
\end{align}
そして、方程式を立てる元となっている要素分割の際の図も見てみましょう。
節点に各要素内の番号と全体の番号が割り当てられていて、要素内の番号（図ではe=3の要素のみ表示）は要素内に1,2,...6といったように割り当てられています。全体の節点の番号は1...17まであります。また、境界要素については、$l=1,2,...5$まであり、各要素内では1,2,3と番号が割り当てられています(煩雑になるため図には書いていません)。これらの$e=1,2,...6$の６つの要素で計算される$\mathbf{K^e},\mathbf{f_e}$、$l=1,2,...5$の５つの境界要素で計算される$\mathbf{q_l}$をまとめて一つの大きな方程式を作りたいのです。つまり、
\mathbf{w}^T\mathbf{K}\mathbf{u} = \mathbf{w}^T (-\mathbf{f} + \mathbf{q})
としたいのです。
この大きな行列を$\mathbf{K}$の作り方を説明します。
$e=3$のとき、$\mathbf{K_e}$の(1,1)要素であれば、要素内の節点番号1の節点は全体番号では7であるため、$\mathbf{K}$の(7,7)要素に$\mathbf{K_e}$の(1,1)要素を足します。$\mathbf{K_e}$の(1,2)要素であれば、要素内の節点番号1の節点は全体番号では7、要素内の節点番号2の節点は全体番号では8であるため、$\mathbf{K}$の(7,8)要素に$\mathbf{K_e}$の(1,2)要素を足します。
これを(1,1)要素から(6,6)要素まで繰り返し行い、更に、$e=1$から$e=6$でも同じことをすることで$\mathbf{K}$が完成します。
$\mathbf{f}$と$\mathbf{q}$に関しても同様に行います。
また、任意の$\mathbf{w}$に対して
\mathbf{w}^T\mathbf{K}\mathbf{u} = \mathbf{w}^T (-\mathbf{f} + \mathbf{q})
が成り立つため
\mathbf{K}\mathbf{u} = -\mathbf{f} + \mathbf{q}
となり、この連立方程式を解くことで$\mathbf{u}$を求めることができます.
次回、プログラム書きます...
pythonかCかfortran03かRustで書こうかなと...
ようやくqiitaらしく...
参考文献
偏微分方程式の数値解法 / 神谷紀生, 北栄輔著 東京 : 共立出版, 1998.3 (三角形一次要素を使った解法がのっています)