So-net無料ブログ作成

[境界要素法プログラムを設計する([計算部]の完成)] [ネコ騙し数学]

[境界要素法プログラムを設計する([計算部]の完成)]

 

1.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

 

2.Declare文を追加する

 前回の結果より、境界積分のために外部関数LengthAngleNormal_ParamNormal_LengthNormal_FootB_Inte_B1B_Inte_B2B_Inte_H1B_Inte_H2を使用するのでした。これらのDeclare文を追加します。

 

 追加する場所はどこか?という問題があります。個人的な趣味では、これらを使うのは「境界要素番号kで境界積分を行うLoop」なので、道具は「実行者」のすぐそばに置いておきたいのですが、個人的趣味ばかり押し出すと嫌われそうなので(^^;)、普通に[計算部]冒頭に置く事にします。

 

3.Gauss掃き出し法のサブルーティン名を決める

 Gauss掃き出し法のサブルーティン名と引数は、Gauss_0(A, z, n, Zero_Order, Return_Code)としておきます。この名前は趣味です!(←舌の根も乾かぬ内に・・・(^^;))。Gauss_0Declare文も[計算部]冒頭に追加します。

 Azの意味はいいと思うので省略します。nは未知数の数(=2*BN_Count),Zero_Orderは正の整数でピボットの0判定条件に使います。つまりピボットの絶対値がEps10^(- Zero_Order)未満だったら、特異行列です。自分はたいてい、Zero_Order8を使います。

ピボットの意味がわからなかったら、ネコ先生の記事を読もう!(^^)

 Return_Codeは計算終了時の状態を表し、Return_Code = 02の値を取るものとします。Return_Code0は「正常終了」(ちゃんと連立方程式が解けた)、Return_Code1は「0の行がある」、Return_Code2は「特異行列である」だとしておきます。これらはユーザーには関係ないです(プログラマーでないユーザーには(^^))。

 だって正常終了しなかったら、開発中のプログラムは間違てるって事ですから。しかし「プログラマー」という、「ユーザーにとっては」まずデバック時に非常に役立ちます(^^)。だから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

nice!(0)  コメント(0) 

nice! 0

コメント 0

コメントを書く

お名前:
URL:
コメント:
画像認証:
下の画像に表示されている文字を入力してください。