WANtaroHP (FEM: 2D Auto-mesh)

Contents

• はじめに
• 三角形要素の自動分割
• 四角形要素の自動分割
• 単純な要素分割事例
• 領域番号定義プログラム f90_fem_domain
• プログラム
• (付録) 矩形領域分割プログラム
• (付録) 節点番号最適化プログラム

はじめに

• 本プログラムは，谷口健男著，「FEMのための要素自動分割　デローニー三角分割法の利用」に掲載のFORTRANプログラムの一部(全部ではない)を，Fortran 90 に書き直したものである．
• 機能は，円孔などの内部境界も含め，任意の形状の2次元領域を三角形要素あるいは四角形要素に自動分割するものである．三角形要素と四角形要素の混在には対応していない．
• 三角形要素用，四角形要素用プログラム共，部分的に細かい要素を作成する機能は有していない．
• 三角形要素では，節点数，要素数，要素-節点関係，要素-領域関係，節点座標が出力され，物性値の異なる領域にも自動対応できる．
• 四角形要素では，入力時は領域毎にデータ入力を行うが，出力時は要素-領域関係が出力されないため，手動で対応する必要がある．
• 四角形要素では，原典と同様に，(1)「三角形要素の生成」，(2)「三角形要素をベースとした四角形要素の生成」，(3)「生成された四角形要素の形状修正」という過程を，それぞれの3つのプログラムで実現している．
• 三角形要素・四角形要素共に要素分割を行うのみで，境界条件や荷重設定を行う機能は含んでいないため，これらは別途設定する必要がある．
• ここで紹介しているプログラム・コードは以下のものである．

• このプログラムでの，外部境界・内部境界の定義は以下のとおり．
  • 外部境界：外部に要素がない，解析領域の最も外側を示す境界線．
  • 内部境界：内側に要素がない，解析領域の最も内側を示す境界線．
  • したがって，内部境界の内側は空洞になる．
• このプログラムで，領域内に物性値が違う複数の閉領域を自動的に設定するのは不可能．
• 領域内に物性値が違う複数の閉領域を設定したい場合，境界内部の節点数 NIB とその数に相当する節点座標を入力し，要素生成後に別途対応する要素にそれぞれの物性を与える作業を行う必要がある．

(参考文献) 谷口健男著：FEMのための要素自動分割　デローニー三角分割法の利用，森北出版（株），2006年8月

更に追加として，自作した，領域番号定義プログラム f90_fem_domain.f90 を紹介している． 上記プログラムの場合，三角形要素では，入力データに従い領域番号が定義されるが，四角形要素では入力時は領域を指定するが最終的には領域定義はなされない． このため，四角形要素の領域番号の定義，および三角形要素において必要な場合に領域を再定義する可能性を考慮し，領域番号定義プログラムを作成・紹介するものである．

三角形要素の自動分割

プログラムの特徴

• 指定した総節点数となるよう三角形要素群が生成される．
• 外部境界，内部境界，境界内部節点として入力した節点は，出力時も同一番号・同一座標が保証される．
• 入力時に指定した以外の節点については法則性はなく自動的に生成・番号付けされる．
• 要素番号については外部境界数で指定した順に要素-領域関係の番号が付されるが，要素番号は法則性はなく自動的に生成・番号付けされる．
• 出力は節点数，要素数，要素-節点関係，要素-領域関係，節点座標であり，物性値の異なる領域にも自動対応できる．
• 入力データ不良等により計算が途中で止まった場合は「Error終了メッセージ」を表示しプログラムを終了する．

f90_fem_mesh3 データ入力書式（半角空白区切りテキストファイル）

INPOINT                           # 総節点数(指定値)
NEX  MIN                          # 外部境界数，内部境界数
ibex(i)                           # 個々の外部境界上の節点数
   .....(i=1 to NEX).....         #
ibin(i)                           # 個々の内部境界上の節点数
   .....(i=1 to MIN).....         # 内部境界数が0なら省略
ibno(i,1) ..... ibno(i,ibex(i))   # 個々の外部境界上の節点番号を1行に入力
   .....(i=1 to NEX).....         #
NOB  NIB                          # 境界上の節点総数，境界内部の節点総数
px(i)  py(i)                      # 節点のx座標およびy座標
   .....(i=1 to NOB+NIB).....     #

注意事項

• 外部境界は時計回りに，内部境界は反時計回りに節点番号を入力する．
• 外部境界上の節点番号ibno(i,j)は循環させず1個手前の節点までの入力とする．
• 節点番号は，(1)外部境界，(2)内部境界，(3)境界内部節点の順につけ，節点座標は節点番号順に入力する．

f90_fem_mesh3 出力ファイル書式（半角空白区切りテキストファイル）

NODE  NELM
mtj(i,1)  mtj(i,2)  mtj(i,3)  idm(i)
        mtj(i,j): 要素-節点関係
        idm(i)  : 要素-領域関係
   .....(i=1 to NELM).....
px(i)  py(i)
        px(i), py(i): 節点(x,y)座標
   .....(i=1 to NODE).....

実行用バッチファイル例

./f90_fem_mesh3 inp_mesh3.txt out_mesh3.txt

ここに，f90_fem_mesh3 は実行プログラム名，inp_mesh3.txt は入力ファイル名，out_mesh3.txt は出力ファイル名である．

四角形要素の自動分割

プログラムの特徴

• 総節点数は指定はできず，自動設定される．
• 外部境界・内部境界上では指定した節点の中間に1節点追加されて四角形要素群が生成される．
• 要素の大きさは境界上の節点間隔と領域毎入力する格子間隔に支配される．
• 外部境界，内部境界，境界内部節点として入力した節点は，出力時も同一番号・同一座標が保証される．
• 入力時に指定した以外の節点については法則性はなく自動的に生成・番号付けされる．
• 出力は節点数，要素数，要素-節点関係，節点座標である．入力時は領域毎にデータ入力を行うが，出力時は要素-領域関係が出力されないため，材料種別の指定は手動で対応する必要がある．
• プログラムでは，f90_fem_mesh4a で「三角形要素の生成」，f90_fem_mesh4b で「三角形要素をベースとした四角形要素の生成」，f90_fem_mesh4c で「生成された四角形要素の形状修正」という過程を順次実行する． 各プログラムへの入出力ファイルは，それぞれのプログラム実行時にコマンドライン引数より入力する．
• 入力データ不良等により計算が途中で止まった場合は「Error終了メッセージ」を表示しプログラムを終了する．

f90_fem_mesh4a データ入力書式（空白区切りテキストファイル）

NEX  MIN                          # 外部境界数，内部境界数
ibex(i)  delx(i)                  # ibex(i): 個々の外部境界上の節点数
   .....(i=1 to NEX).....         # delx(i): 各部分領域の格子間隔
ibin(i)                           # 個々の内部境界上の節点数
   .....(i=1 to MIN).....         # 内部境界数が0なら省略
ibno(i,1) ..... ibno(i,ibex(i))   # 個々の外部境界上の節点番号を1行に入力
   .....(i=1 to NEX).....         #
NOB  NIB                          # 境界上の節点総数，境界内部の節点総数
px(i)  py(i)                      # 節点のx座標およびy座標
   .....(i=1 to NOB+NIB).....     #

注意事項

• 外部境界は時計回りに，内部境界は反時計回りに節点番号を入力する．
• 各部分領域の格子間隔delx(i)は，外部境界で指定した座標間隔と同程度が良いようである．
• 外部境界上の節点番号ibno(i,j)は循環させず1個手前の節点までの入力とする．
• データ入力における三角形要素生成プログラムとの違いは，総節点番号(INPNODT)の入力がないこと，各部分領域の格子間隔(delx)の入力を必要とすることの2点のみである．

f90_fem_mesh4 出力ファイル書式（空白区切りテキストファイル）

NODE  NELM    # 総節点数，総要素数
mtj(i,1)  mtj(i,2)  mtj(i,3)  jac(i,1)  jac(i,2) jac(i,3) idm(i)
        mtj(i,j): 要素-節点関係
        jac(i,j): 要素-隣接要素関係
        idm(i)  : 要素-領域関係
   .....(i=1 to NELM).....
px(i)  py(i)  ifix(i)
        px(i), py(i): 節点(x,y)座標
        ifix(i): 固定節点(=1)・可動節点(=0)の判別
   .....(i=1 to NODE).....

f90_fem_mesh4b 出力ファイル書式（空白区切りテキストファイル）

NODE  NELM    # 総節点数，総要素数
mtj(i,1)  mtj(i,2)  mtj(i,3)  mtj(i,4)
        mtj(i,j): 要素-節点関係
   .....(i=1 to NELM).....
px(i)  py(i)  ifix(i)
        px(i), py(i): 節点(x,y)座標
        ifix(i): 固定節点(=1)・可動節点(=0)の判別
   .....(i=1 to NODE).....

f90_fem_mesh4c 出力ファイル書式（空白区切りテキストファイル）

NODE  NELM    # 総節点数，総要素数
mtj(i,1)  mtj(i,2)  mtj(i,3)  mtj(i,4)
        mtj(i,j): 要素-節点関係
   .....(i=1 to NELM).....
px(i)  py(i)
        px(i), py(i): 節点(x,y)座標
   .....(i=1 to NODE).....

実行用バッチファイル例

./f90_fem_mesh4a inp_mesh4a.txt out_mesh4a.txt
./f90_fem_mesh4b out_mesh4a.txt out_mesh4b.txt
./f90_fem_mesh4c out_mesh4b.txt out_mesh4c.txt

手動でファイル作成が必要なのは，最初の入力ファイル inp_mesh4a.txt のみであり，最終的な成果は出力ファイルは out_mesh4c.txt である． その他のファイルは，プログラム間でデータのやりとりに用いているのみであるため，常に同一名を使用しても中間状態の確認を行う必要性がない限り問題はない．

単純な要素分割事例

下図に示す3個の領域（I，II，III）と領域 II の内部に 1 個の内部境界を有するモデルを要素分割する． 境界上の総節点数は 10 とし，内部節点 1 個を領域 I 内に配置する．

三角形要素分割入力データ事例 (Case-T1)

三角形要素においては，上記座標を基本として，Case-T2 では指定節点数を INPNODT=100，Case-T3 では指定節点数を INPNODT=1000 としている．

三角形要素分割の場合，境界上の入力点ピッチが粗くても，節点数を大きくすることにより，概ね均等な三角形要素が生成される．

四角形要素分割入力データ事例 (Case-Q1)

四角形要素のケース Case-Q1 では三角形要素のケース Case-T1 と同一の座標を指定している． このケースでは，節点番号28の座標位置が不適切で，要素番号 15 で形状不正となっている．

Case-Q2 では，領域IIの上辺と下辺に1点づつ節点を追加した． これにより，要素の形状不正はなくなっている．

Case-Q3 では，境界上の節点の指定を座標値の増分で5以下となるように指定した． これにより，要素の扁平さは緩和されている．

また，要素の生成過程より，初期に生成された三角形要素に対し，各三角形の重心と辺の中点に節点を増設し，1つの三角形を3つの四角形要素に分割することにより要素が生成されていく過程が理解できる． すなわち，初期に与えた外部境界点の隣り合う2点の中点に必ず1節点が追加されることになる． よって外部境界点座標を規則的に与えた場合，その点間隔と同じもしくは半分の寸法を格子間隔 delx として指定すると比較的イメージに近い分割が行える．

このように，四角形要素の場合は，要素形状は入力節点に大きく左右されるため，最終的に得たい要素の寸法や形状を予め考慮して入力節点数や座標を設定する必要がある．

Case-T1 節点番号図 節点数指定：20 要素数：27
Case-T1 要素番号図 節点数指定：20 要素数：27
Case-T2 節点番号図 節点数指定：100 要素数：162
Case-T2 要素番号図 節点数指定：100 要素数：162
Case-T 要素図 節点数指定：1000 要素数：1870
三角形要素分割図

三角形要素生成
四角形要素生成
四角形要素整形結果 52節点，42要素
四角形要素整形結果 52節点，42要素
四角形要素分割図 (Case-Q1) 図中黒丸は，初期入力節点

三角形要素生成
四角形要素生成
四角形要素整形結果 60節点，48要素
四角形要素整形結果 60節点，48要素
四角形要素分割図 (Case-Q2) 図中黒丸は，初期入力節点

三角形要素生成
四角形要素生成
四角形要素整形結果 202節点，180要素
四角形要素整形結果 202節点，180要素
四角形要素分割図 (Case-Q3) 図中黒丸は，初期入力節点

領域番号定義プログラム f90_fem_domain

三角形要素では，入力データに従い領域番号が定義されるが，四角形要素では入力時は領域を指定するが最終的には領域定義はなされない．

このため，四角形要素の領域番号の定義，および三角形要素において必要な場合に領域を再定義する可能性を考慮し，領域番号定義プログラムを作成した．

このプログラムの特徴は以下のとおり．

• 三角形要素・四角形要素ともに対応している．ただしコマンドライン引数で1要素当たりの節点数を入力する必要がある．
• 指定した領域にある要素が含まれるか否かは，以下の理論により判定している．
  • 各要素の重心座標を計算する．
  • 要素重心座標が多角形で構成される領域もしくは境界線上に含まれる場合，要素重心と多角形の各点を結んだ線分で構成される三角形の要素重心側の角度の合計は360度もしくは-360度となる．
  • 要素重心座標が多角形で構成される領域に含まれない場合，要素重心と多角形の各点を結んだ線分で構成される三角形の要素重心側の角度の合計は0度になる．
  • 要素重心座標が多角形を構成する点と一致する場合，その点との距離が0となり，なす角を計算できない．この場合角度を0度とみなすことにより，角度の合計は0を超過し360度未満となる．
• 領域境界点は時計回り・反時計回りいずれも可．境界が閉じる1点前までの座標を入力する．
• 後から領域定義したものが優先される．
• 領域境界線上の点はその領域に含まれる．ただし後からの定義が優先される．(プログラム上，合計角度が0とみなせない場合は，境界上も含め，着目要素の重心はその多角形領域に含まれると判断する)

実行用バッチコマンド書式

./f90_fem_domain nod fnameR.txt fnameI.txt fnameW.txt

f90_fem_domain 領域指定データファイル入力書式（半角空白区切りテキストファイル）

MATEL                           # 指定する領域数
nt(1)                           # 領域 1 で指定する座標の数
x_dom(1,j)  y_dom(1,j)          # 領域 1 での(x,y)座標
   .....(j=1 to nt(1).....      #
   .....             .....      #
nt(k)                           # (必要数繰り返し)
x_dom(k,j)  y_dom(k,j)          #
   .....(j=1 to nt(k).....      #
   .....             .....      #
nt(MATEL)                       # 領域 MATEL で指定する座標の数
x_dom(MATEL,j)  y_dom(MATEL,j)  # 領域 MATEL での(x,y)座標
   .....(j=1 to nt(MATEL))..... #

f90_fem_domain 出力ファイル書式（半角空白区切りテキストファイル）

nod  NODT  NELT  MATEL
kakom(ne,1)  kakom(ne,2)  kakom(ne,3)  matno(ne)
     kakom(ne,j): 要素−節点関係 (四角形要素では j=1 to 4)
     matno(ne)  : 要素−領域関係 (材料種別番号)
   .....(ne=1 to NELT).....
x(i)  y(i)
     x(i),y(i): 節点の(x,y)座標
   .....(i=1 to NODT).....

出力事例

f90_fem_domain 領域指定データファイル入力事例

Case-T3 1000節点 1870要素
Case-Q3 202節点 180要素
要素領域設定・着色事例

プログラム

Fortran programs

Shell scripts

Input data samples

(付録) 矩形領域分割プログラム

テスト用データ作成などに便利な，矩形領域を四角形要素に分割するプログラムを示します．

./f90_rectmesh aa bb nn mm x0 y0 > out_mesh.txt

(付録) 節点番号最適化プログラム

自動メッシュ出力の段階で節点番号最適化を行うプログラムです．

./f90_num_mesh nod fnameR fnameW