［境界要素法プログラムを設計する（[計算部]の完成）］ [ネコ騙し数学]
［境界要素法プログラムを設計する（[計算部]の完成）］
１．Free Termの処理
もう一度［内点方程式から境界方程式へ（ラプラス型）］のFree Termの項を考えます。
図-1でε→0の極限を取った場合、式(1)のψ(ξ，η)の(ξ，η)は、節点j+1に一致します。つまりψ(ξ，η)→ψj+1です。よって、
となります。式(2)の最後の形のψj+1の係数の2項目が、境界積分関数で与えられるGam_12です。
節点jが角点の場合、
になるのでした。kj+1は節点j+1の内角です。ところで今は、角点のない解析領域に対応するプログラムを作成中です。節点が角点でない場合、［境界方程式に関する留意点（ラプラス型）］によればkj+1/2/π→1/2におきかえるのが一般的でした。滑らかな境界点の場合、kj+1＝πと考えられるからです（β(j+1)＝β(j)だから）。
よってメインの「REM 境界要素番号kで境界積分を行うLoop」で全体係数マトリックスA、
を作った後、式(5)のBブロックの対角成分を1/2に上書きする必要があります。上記のLoopの後に、次のコードを追加します。
（「境界要素番号kで境界積分を行うLoop」の後）
REM Bマトリックスの滑らかな境界節点のFree Termを1/2に上書きするLoop
For i = 1 to BN_Count
A(i, i) = 1 / 2
Next i
２．Declare文を追加する
前回の結果より、境界積分のために外部関数Length，Angle，Normal_Param，Normal_Length，Normal_Foot，B_Inte_B1，B_Inte_B2，B_Inte_H1，B_Inte_H2を使用するのでした。これらのDeclare文を追加します。
追加する場所はどこか？という問題があります。個人的な趣味では、これらを使うのは「境界要素番号kで境界積分を行うLoop」なので、道具は「実行者」のすぐそばに置いておきたいのですが、個人的趣味ばかり押し出すと嫌われそうなので(^^;)、普通に[計算部]冒頭に置く事にします。
３．Gauss掃き出し法のサブルーティン名を決める
Gauss掃き出し法のサブルーティン名と引数は、Gauss_0(A, z, n, Zero_Order, Return_Code)としておきます。この名前は趣味です！（←舌の根も乾かぬ内に・・・(^^;)）。Gauss_0のDeclare文も[計算部]冒頭に追加します。
A，zの意味はいいと思うので省略します。nは未知数の数（＝2*BN_Count），Zero_Orderは正の整数でピボットの0判定条件に使います。つまりピボットの絶対値がEps＝10^(- Zero_Order)未満だったら、特異行列です。自分はたいてい、Zero_Order＝8を使います。
ピボットの意味がわからなかったら、ネコ先生の記事を読もう！(^^)。
Return_Codeは計算終了時の状態を表し、Return_Code = 0～2の値を取るものとします。Return_Code＝0は「正常終了」（ちゃんと連立方程式が解けた）、Return_Code＝1は「0の行がある」、Return_Code＝2は「特異行列である」だとしておきます。これらはユーザーには関係ないです（プログラマーでないユーザーには(^^)）。
だって正常終了しなかったら、開発中のプログラムは間違てるって事ですから。しかし「プログラマー」という、「ユーザーにとっては」まずデバック時に非常に役立ちます(^^)。だからReturn_Codeの値によってはメッセージを出し、無駄な出力も行わないで強制終了するような制御も入れておきましょう。
リリースしたらその制御は発動しない事を願うのみですが、間違って発動したら「旦那、どんなメッセージが出やした？」ってお客さんに犬のようにきけるんですよ(^^)。それは間違い修正やメンテナンスに対する、非常に得難い情報になります。そしてFortranではそのような行為を「旗を立てる」と言いました(^^)。
以上まとめると以下です。
REM [計算部]
External Function Length
External Function Angle
External Function Normal_Param
External Function Normal_Length
External Function Normal_Foot
External Function B_Inte_B1
External Function B_Inte_B2
External Function B_Inte_H1
External Function B_Inte_H2
External Sub Gauss_0
REM 境界要素番号kで境界積分を行うLoop
For i = 1 to BN_Count
REM 特異点番号iの節点iを特定
x0 = BN(i, 1) REM 節点iの(x，y)座標
y0 = BN(i, 2)
For k = 1 to B_Count
REM 要素番号kからEN_B(k，2)に従い節点を特定し、要素ごとに境界積分を実行
j1 = EN_B(k，1) REM 要素kの始端の節点番号
j2 = EN_B(k，2) REM 要素kの終端の節点番号
x1 = BN(j1, 1) REM 節点j1の(x，y)座標
y1 = BN(j1, 2)
x2 = BN(j2, 1) REM 節点j2の(x，y)座標
y2 = BN(j2, 2)
REM (x0，y0)，(x1，y1)，(x2，y2)から、要素kの境界積分値を計算
REM 積分ﾊﾟﾗﾒｰﾀ計算
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)
t = Normal_Param(x0, y0, x1，y1，x2，y2)
h = Normal_Length(x1，y1，x2，y2, t)
s = Normal_Foot(x0, y0, x1，y1，x2，y2, t)
REM ψj1に関する結果をB1，qj1に関する結果をH1で表す
REM ψj2に関する結果をB2，qj2に関する結果をH2で表す
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)
REM 上記結果に基づき、係数ﾏﾄﾘｯｸｽBとHを作成
jj1 = BN_Z(j1，1) REM ψj1の解ﾍﾞｸﾄﾙzの中での位置
jjj1 = BN_Z(j1，2) REM qj1の解ﾍﾞｸﾄﾙzの中での位置
jj2 = BN_Z(j2，1) REM ψj2の解ﾍﾞｸﾄﾙzの中での位置
jjj2 = BN_Z(j2，2) REM qj2の解ﾍﾞｸﾄﾙzの中での位置
A(i, jj1) = A(i, jj1) - B1
A(i, jjj1) = A(i, jjj1) + H1
A(i, jj2) = A(i, jj2) - B2
A(i, jjj2) = A(i, jjj2) + H2
Next k
Next i
REM Bマトリックスの滑らかな境界節点のFree Termを1/2に上書きするLoop
For i = 1 to BN_Count
A(i, i) = 1 / 2
Next i
REM 領域要素番号kで領域積分を行うLoop
For i = 1 to BN_Count
REM 特異点番号iの節点iを特定
h = 0
For k = 1 to R_Count
REM 要素番号kからEN_R(k，4)に従い節点を特定し、要素ごとに領域積分を実行
REM 特異点iと要素節点情報から積分値w0を計算
h = h + w0
Next k
REM 上記結果に基づき、既知ﾍﾞｸﾄﾙwを作成
z(i) = h
Next i
REM 境界節点番号jで境界条件位置を指定するLoop
For i = BN_Count + 1 to 2 * BN_Count
REM 行番号から節点番号を特定する
j = i - BN_Count
REM 節点番号jからBN_BP(j，2)とBN_Z(j，2) に従い、(D1 D2)ﾌﾞﾛｯｸを作成
REM 節点番号jに従いBN_BV(j，2)の値をwﾌﾞﾛｯｸの後に並べ、(d1 d2)ﾌﾞﾛｯｸを作成
If BN_BP(j，1) = 1 Then
REM ψjが既知が場合
jj = BN_Z(j，1) REM ψjの解ﾍﾞｸﾄﾙzの中での位置
A(i, jj) = 1 REM ψjの値を与えるので、Aの成分は1
z(i) = BN_BV(jj，1) REM ψjの値
ElseIf BN_BP(j，2) = 1 Then
REM qjが既知が場合
jj = BN_Z(j，2) REM qjの解ﾍﾞｸﾄﾙzの中での位置
A(i, jj) = 1 REM qjの値を与えるので、Aの成分は1
z(i) = BN_BV(jj，2) REM qjの値
End If
Next i
REM 全体係数行列が出来たので、Gaussの掃き出し法にかけ、未知ﾍﾞｸﾄﾙzを求める
Zero_Order = 8
Call Gauss_0(A, z, 2 * BN_Count, Zero_Order, Return_Code)
If Return_Code = 0 Then
Print "正常終了しました。"
ElseIf Return_Code = 1 Then
Print "0の行があります。"
End
ElseIf Return_Code = 2 Then
Print "特異行列です。"
End
End If
REM [外部手続き]
REM 境界積分ﾊﾟﾗﾒｰﾀ ****************************************
External Function Length(x1，y1，x2，y2)
Length = Sqr((x2 - x1)^2 + (y2 - y1)^2)
End Function
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
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
REM 境界積分 **********************************************
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
REM Gauss掃き出し法 **********************************************
External Sub Gauss_0(A, z, Zero_Order, Return_Code)
REM まだ何にもしない
End Sub
