［境界要素法プログラムを設計する（境界積分ｻﾌﾞﾙｰﾃｨﾝ：実装編）］ [ネコ騙し数学]
［境界要素法プログラムを設計する（境界積分ｻﾌﾞﾙｰﾃｨﾝ：実装編）］
１．積分パラメータ関数の実装
前々回、メインの中に与えた境界積分サブルーティン（関数）は、以下でした。
REM (x0，y0)，(x1，y1)，(x2，y2)から、要素kの境界積分値を計算
REM ψj1に関する結果をB1，qj1に関する結果をH1で表す
REM ψj2に関する結果をB2，qj2に関する結果をH2で表す
Lk = Length(x1，y1，x2，y2)
Beta_k = Angle(x2 - x1, y2 - y1, Lk)
r1 = Length(x0，y0，x1，y1)
r2 = Length(x0，y0，x2，y2)
Gamm_1 = Angle(x1 - x0, y1 - y0, r1)
Gamm_2 = Angle(x2 - x0, y2 - y0, r2)
h = Normal_Length(x2 - x1, y2 - y1, x0，y0)
s = Normal_Foot(x2 - x1, y2 - y1, x0，y0, h)
B1 = B_Inte_B1(Lk, r1, r2, Beta_k, Gamm_1, Gamm_2, h, s)
B2 = B_Inte_B2(Lk, r1, r2, Beta_k, Gamm_1, Gamm_2, h, s)
H1 = B_Inte_H1(Lk, r1, r2, Beta_k, Gamm_1, Gamm_2, h, s)
H2 = B_Inte_H2(Lk, r1, r2, Beta_k, Gamm_1, Gamm_2, h, s)
最初にLength関数ですが、明らかに次でOKです。
External Function Length(x1，y1，x2，y2)
Length = Sqr((x2 - x1)^2 + (y2 - y1)^2)
End Function
前々回に予告したように、関数とサブルーティンはローカル変数を持てる、外部手続きで統一します。Normal_LengthとNormal_Footの計算に疑問の余地はないのですが、ちょっと長そうなので後にします。
次にBeta_kの扱いですが、前回の検討からGamm_1，Gamm_2の特殊ケースとして扱って良いとわかりました。という訳でGamm_1とGamm_2の計算です。これらは単独に計算しては駄目で、Gamm_2－Gamm_1の形で計算する必要がありました。それをGam_12で表します。
まず特異点(x0，y0)が要素節点(x1，y1)と(x2，y2)に一致しないケースから考え始めると、わかりやすいと思います。以下、前回を参照しつつ・・・。
計算開始の前提として、図-1のr1，r2および(x0，y0)，(x1，y1)と(x2，y2)は引数渡しされるとします。また特異点と要素節点の一致／不一致判定のために、0距離判定値Epsも必要です。
External Function Angle(x0, y0, x1，y1，x2，y2, r1, r2, Eps)
If r1 < Eps Then
x01 = 1
y01 = 0
x02 = x2 - x1
y02 = y2 - y1
r1 = 1
ElseIf r2 < Eps Then
x01 = x1 - x2
y01 = y1 - y2
x02 = 1
y02 = 0
r2 = 1
Else
x01 = x1 - x0
y01 = y1 - y0
x02 = x2 - x0
y02 = y2 - y0
End If
Gum_12 = Acos((x01 * x02 + y01 * y02) / r1 / r2)
Sign = Sgn(x01 * y02 - x02 * y01)
Gum_12 = Gum_12 * Sign
Angle = Gum_12
End Function
r1 < Eps（r1＝0）の時は、(x01，y01)＝(1，0)とするのでした。この時(x0，y0)は(x1，y1)に一致するとみなせるので、(x02，y02)＝(x2 - x1，y2 - y1)です。そしてr1＝1とします。
r2 < Eps（r2＝0）の時は、(x02，y02)＝(1，0)とするのでした。この時(x0，y0)は(x2，y2)に一致するとみなせるので、(x01，y01)＝(x1 - x2，y1 - y2)です。そしてr2＝1とします。
以上より、積分パラメータを与える関数呼び出しは、
Lk = Length(x1，y1，x2，y2)
r1 = Length(x0，y0，x1，y1)
r2 = Length(x0，y0，x2，y2)
Gam_12 = Angle(x0, y0, x1，y1，x2，y2, r1, r2, Eps)
h = Normal_Length(x0, y0, x1，y1，x2，y2)
s = Normal_Foot(x0, y0, x1，y1，x2，y2)
になります（Normal_Length，Normal_Footの引数が前回のままでは不味いと気付いた(^^;)）。
Normal_Lengthの実装は以下です。Normal_Lengthは図-1の垂線距離hを求めます。
a = x2 - x1
b = y2 - y1
とすると、(x1，y1)と(x2，y2)を通る直線は、y＝b/a*xに平行です。a＝0の場合も考慮して、これをbx－ay＝0の形にし、(x1，y1)を通る事を考慮すると要素kの直線の方程式は、b(x－x1)－a(y－y1)＝0になります（なんか受験を思い出すなぁ～(^^)）。直線の式からその垂線ベクトルの方位は、明らかに(b，－a)です。従って点(x0，y0)からの垂線は、tを任意の実数として、
と書けます。(1)がb(x－x1)－a(y－y1)＝0に載る条件は、(1)を代入し、
(2)をtについて解くと、
従ってhは、tが(3)の値の時、
tを利用すれば垂線の足もわかるので、図-1のsも出せます。垂線の足の座標は、tを(3)の値にした(1)です。よってsは、
から、
です。
こうなるとNormal_LengthとNormal_Footで同じパラメータtを利用したいですよね。こうしましょう。
h = Normal_Length(x1，y1，x2，y2, t)
s = Normal_Foot(x0, y0, x1，y1，x2，y2, t)
External Function Normal_Param(x0, y0, x1，y1，x2，y2)
a = x2 - x1
b = y2 - y1
t = (b * (x1 - x0) - a * (y1 - y0)) / (a^2 + b^2)
Normal_Param = t
End Function
External Function Normal_Length(x1，y1，x2，y2, t)
a = x2 - x1
b = y2 - y1
h = Abs(t) * Sqr(a^2 + b^2)
Normal_Length = h
End Function
External Function Normal_Foot(x0, y0, x1，y1，x2，y2, t)
a = x2 - x1
b = y2 - y1
s = Sqr((t * b + x0 - x1)^2 + (-t * a + y0 - y1)^2)
Normal_Foot = s
End Function
２．境界積分関数の実装
以上で境界積分に必要な積分パラメータは全て揃いました。
B1 = B_Inte_B1(Lk, r1, r2, Gam_12, h, s, Eps)
B2 = B_Inte_B2(Lk, r1, r2, Gam_12, h, s, Eps)
H1 = B_Inte_H1(Lk, r1, r2, Gam_12, h, s, Eps)
H2 = B_Inte_H2(Lk, r1, r2, Gam_12, h, s, Eps)
で良いはずです。引数の最後にEpsが入るのは、これらでも特異点と要素節点の一致／不一致判定が必要だからです。
後はB_Inte_？の中に、［内点方程式の積分公式（ラプラス型）］と［内点方程式から境界方程式へ（ラプラス型）］で導出した式を並べるだけですが、関数の内部構造はExternal Function Angleと同一構造になるとわかります。なのでIfブロックの形を決め、「関数名 = 値」の行を書いてから、計算式を並べる事をお奨めします。
External Function B_Inte_B1(Lk, r1, r2, Gam_12, h, s, Eps)
If r1 < Eps Then
B1 = Gam_12 / 2 /Pai
ElseIf r2 < Eps Then
B1 = 0
Else
B1 = Gum_12 - 1 / Lk * (s * Gum_12 + h * (Log(r2) - Log(r1)))
B1 = B1 / 2 /Pai
End If
B_Inte_B1 = B1
End Function
External Function B_Inte_B2(Lk, r1, r2, Gam_12, h, s, Eps)
If r1 < Eps Then
B2 = 0
ElseIf r2 < Eps Then
B2 = Gam_12 / 2 /Pai
Else
B2 = 1 / Lk * (s * Gum_12 + h * (Log(r2) - Log(r1)))
B2 = B1 / 2 /Pai
End If
B_Inte_B2 = B2
End Function
External Function B_Inte_H1(Lk, r1, r2, Gam_12, h, s, Eps)
If r1 < Eps Then
H1 = Lk / 4 / Pai * (Log(Lk) - 3 / 2)
ElseIf r2 < Eps Then
H1 = Lk / 4 / Pai * (Log(Lk) - 1 / 2)
Else
H1 = -1 / 4 / Pai / Lk * (r2^2 * Log(r2) - r1^2 * Log(r1) - 1 / 2 *((Lk - s)^2 - s^2))
H1 = H1 + 1 / 2 / Pai * (1 - s / Lk) * (h * Gum_12 + (Lk - s) * Log(r2) + s * Log(r1) - Lk)
End If
B_Inte_H1 = H1
End Function
External Function B_Inte_H2(Lk, r1, r2, Gam_12, h, s, Eps)
If r1 < Eps Then
H2 = Lk / 4 / Pai * (Log(Lk) - 1 / 2)
ElseIf r2 < Eps Then
H2 = Lk / 4 / Pai * (Log(Lk) - 3 / 2)
Else
H2 = 1 / 4 / Pai / Lk * (r2^2 * Log(r2) - r1^2 * Log(r1) - 1 / 2 *((Lk - s)^2 - s^2))
H2 = H2 + 1 / 2 / Pai * s / Lk * (h * Gum_12 + (Lk - s) * Log(r2) + s * Log(r1) - Lk)
End If
B_Inte_H2 = H2
End Function
以上の結果は、[計算部]と[外部手続き]のところへは、まだはめ込みません。[計算部]では、まだやる事があるからです（もぉ～、境界要素法って面倒くさいなぁ～(^^;)）。
