WANtaroHP (F90 FEM contour by GMT)

Contents

• 熱伝導解析結果の表示 (メッシュ図および3次元モノクロコンター)
• 浸透流解析結果の表示 (メッシュ図および2次元カラーコンター)
• 2次元応力解析結果の表示 (メッシュ図および2次元カラーコンター)
• 凡例を含む2次元プロット (2次元プロット)

このページでは，GMTを使ってFEMによる計算結果を表現する方法を紹介しています．

GMTで作成されたepsファイルは，ImageMagickを用いてpng形式に変換されています．

GMTによるFEMメッシュとコンターマップ作成

• Fortran プログラムは以下の事項を実行します．
  • FEMプログラム出力の読み込み
  • データ処理とGMTのためのデータファイル出力
  • GMTを実行するためのバッチファイル出力 (filename: bat_gmt_xxxx.bat)
  • Fortran90の組み込みサブルーチン ''system('bat_gmt_xxxx.bat')'' を用いたGMTバッチファイルの実行

GMTによる凡例を含む2次元プロット

• gnuplotでは凡例の挿入は簡単に行えますが，GMTでは凡例の挿入は簡単ではありません． そこでGMTにおいて凡例付きの2次元プロットを行うためのバッチファイルを作成するプログラムを作ってみました．
• このプログラムは以下を考慮して作られています．
  • 基になるプロットに凡例を重ね合わせる．
  • プロット記号と文字数を用いて凡例領域の寸法を定義する．
  • 凡例領域を白で塗りつぶし，記号と文字列を描画する．
  • gnuplotのように凡例領域の位置を変えることができる．

熱伝導解析結果の表示

概要

このプログラムは，GMTによりFEMメッシュと温度分布の3次元表示を行うプログラムです．

Fortran プログラム (f90_GMTheat.f90)

若干の説明を下に示します．

• GMTコマンド'psxyz'を用いて，温度分布の3次元表示を行います．
• xおよびy座標は，要素の平面座標を示します．z座標は指定された時刻における温度を示します．
• 出力ファイル名は，プログラム内で自動設定されます．

以下のファイルはプログラム実行前に準備されていなければなりません．

FEM解析の出力ファイルは，管理人作：'f90_FEM_HEAT.exe'による出力ファイルです．

処理を実行するバッチファイルを下に示します． Fortranプログラムで読み込まれるいくつかのファイルがバッチファイルの中で作成されています．

001  gfortran -o f90_GMTheat.exe f90_GMTheat.f90
002
003  rem definition of mesh color (1:MATEL)
004  echo 255,192,203 > _col_mesh.txt
005  echo 255,250,205 >> _col_mesh.txt
006
007  rem definition of basic data
008  rem integer::ind,ibc
009  rem real(4)::baselength,baseheight
010  rem real(4)::bavx,bfvx
011  rem real(4)::bavz,bfvz
012  rem character::strx*100,stry*100
013  echo 0 0 > _inp_base.txt
014  echo 10.0 5.0 >> _inp_base.txt
015  echo 10.0 5.0 >> _inp_base.txt
016  echo 5.0 1.0 >> _inp_base.txt
017  echo "x-distance (m)" >> _inp_base.txt
018  echo "y-distance (m)" >> _inp_base.txt
019
020  f90_GMTheat out_CPT_ondom.csv
021
022
023  del _*

001

プログラム 'f90_GMTheat.f90'のコンパイル

004-005

ファイル'_col_mesh.txt.' の作成． メッシュ色はこのファイルで定義され，3種類の数値は 'R G B' を示す． この事例では，'echo'コマンドを用いて新しいファイルを作成しています．また2材料が使われている(パラメータMATEL=2)ので，2種類の色が定義されています．

013-018

ファイル '_inp_base.txt'の作成． このファイルはGMTのためのパラメータを定義します．

ind	FEMメッシュ図に節点番号と要素番号を描画する (0: 描画しない, 1: 描画する)
ibc	ＦＥＭメッシュ図に境界条件を描画する (0: 描画しない, 1: 描画する)
baselength	x軸とy軸の寸法(単位はcm)
baseheight	z軸の寸法(単位はcm)
bavx	GMTのx軸とy軸に対する-Bオプションのaに続く値
bfvx	GMTのx軸とy軸に対する-Bオプションのfに続く値
bavz	GMTのz軸に対する-Bオプションのaに続く値
bfvz	GMTのz軸に対する-Bオプションのfに続く値
strx	x軸のラベル．文字列のためダブルプライム(")が必要
stry	y軸のラベル．文字列のためダブルプライム(")が必要

注意) z軸のラベルは 'Temperature'に固定されています．

020

プログラム'f90_GMTheat.exe'の実行．このexeファイルの機能は，バッチファイル'bat_gmt_xyz.bat'を作成し，組み込みサブルーチン'system(***)'によりこれを実行するものです．

FEM解析結果が入った入力ファイル'out_CPT_ondom.csv'は，Windowsのコマンドプロンプトのコマンドラインより読み込まれます．

出力図は以下に示すとおり．

Filename	Description	Remarks
fig_gmt_mesh.eps	FEMメッシュ図(eps)	出力ファイル名は固定
fig_gmt_xyz_xx.eps	3次元温度分布図(eps)	出力ファイル名は固定．'xx'は時刻を意味します.

FEMメッシュにおける境界条件の記号
Symbol	Description
○	温度指定境界
△	熱伝達境界

023

作業ファイルの削除

プログラムと出力図

ソースプログラムと実行用バッチファイル

トンネル周辺の温度の出力サンプル

以下のサンプルは岩盤内トンネル周辺の温度変化を示します．年平均温度は 10°C, 年間の温度変化振幅は ±10°Cです．

浸透流解析結果の表示

概要

このプログラムは，GMTにより，FEMメッシュ・圧力水頭コンター・全水頭コンターおよび流速ベクトルを示すプログラムです．

Fortran プログラム (f90_GMTseep.f90)

簡単な説明を以下に示します．

• FEMメッシュは，以下の繰り返しにより描画されます．
  • 各要素の座標指定
  • 各要素の塗りつぶし
  • 要素境界の描画
• コンターマップは以下の手順で作成されます．
  • 'grdmask'コマンドを用い，要素が存在する領域を定義する mask ファイルを作成
  • 'surface'コマンドによりコンター作成のための grid ファイルを作成
  • 'grdmath'コマンドにより mask ファイルと grid ファイルを合成し，プロット用 grid ファイルを作成
  • 'grdcontour'コマンドによりコンターを作成し，'grdimage'コマンドにより着色する
• コンター作成により得られた圧力水頭あるいは全水頭の着色図に流速ベクトルの矢印を重ね合わせることにより作成されます．
• さらに，流速ベクトルが全ての要素に描かれると図が醜くなるので，このサンプルでは流速ベクトルは3要素に1回描かれています．

以下のファイルはプウログラム実行前に準備されていなければなりません．

FEM解析の出力ファイルは，管理人作 'f90_FEM_SEEP.exe' による出力ファイルです．

処理を実行するためのバッチファイルを以下に示します． Fortranプログラムで読み込まれるいくつかのファイルは，このバッチファイルの中で作成されます．

001  gfortran -o f90_GMTseep.exe f90_GMTseep.f90
002
003  copy 0_cpt_dam.cpt 0_cpt.cpt
004
005  rem define of mesh color (1:MATEL)
006  echo 255 255 204 > _col_mesh.txt
007  echo   0 255 255 >> _col_mesh.txt
008  echo 255   0 255 >> _col_mesh.txt
009  echo   0 255 255 >> _col_mesh.txt
010  echo 255 255 204 >> _col_mesh.txt
011
012  rem integer::ind,ibc,iph,ivc,isc
013  rem real(4)::baselength,grdspc
014  rem real(4)::bav,bfv
015  rem real(4)::conint,annint
016  rem character::strx*100,strz*100
017  echo 0 1 1 1 1 > _inp_base.txt
018  echo 12.0 0.1 >> _inp_base.txt
019  echo 10.0 5.0 >> _inp_base.txt
020  echo  1.0 5.0 >> _inp_base.txt
021  echo "Distance (m)" >> _inp_base.txt
022  echo "Altitude (m)" >> _inp_base.txt
023
024  f90_GMTseep out_fem_zone4m.csv
025  copy fig_gmt_mesh.eps fig_gmt_mesh_dam.eps
026  copy fig_gmt_conp.eps fig_gmt_conp_dam.eps
027  copy fig_gmt_vect.eps fig_gmt_vect_dam.eps
028
029  del fig_gmt_mesh.eps
030  del fig_gmt_conp.eps
031  del fig_gmt_conh.eps
032  del fig_gmt_vect.eps
033  del _*

001

プログラム 'f90_GMTseep.f90'のコンパイル

003

事前に準備された '0_cpt_dam.cpt' を '0_cpt.cpt'にコピー

006-010

ファイル '_col_mesh.txt' の作成． メッシュ色はこのファイルで定義されます．このファイル中の3つの数値は'R G B'を意味します． このサンプルでは，'echo'コマンドを用いて，新しいファイルを作成しています．また5種類の材料が用いられている(MATEL=5)ため，5種類の色が定義されています．

017-022

ファイル '_inp_base.txt'の作成．このファイルはGMTのためのパラメータを定義します．

ind	メッシュ図への節点番号・要素番号の描画 (0: 描画しない, 1: 描画する)
ibc	メッシュ図への境界条件の描画 (0: 描画しない, 1: 描画する)
iph	圧力水頭コンターへの全水頭線の重ねがき，あるいはその逆 (0: 行わない, 1: 行う)
ivc	速度ベクトル図に重ねがきするコンターの種類 (1: 圧力水頭コンター, 2: 全水頭コンター)
isc	コンター用スケールマップの描画 (0: 描画しない, 1: 描画する)
baselength	図の寸法(単位：cm)
grdspc	'grd'ファイルにおけるメッシュ間隔
bav	GMTのx軸およびy軸に対する-Bオプションのaに続く数値
bfv	GMTのx軸およびy軸に対する-Bオプションのfに続く数値
conint	'grdcontour'コマンドの-Cオプションにおけるコンター間隔の値
annint	'grdcontour'コマンドの-Aオプションにおける目盛間隔の値
strx	x(水平)軸ラベル．文字列のためダブルプライム(")が必要
strz	z(垂直)軸ラベル．文字列のためダブルプライム(")が必要

024

プログラム 'f90_GMTseep.exe' の実行．このexeファイルの機能は，バッチファイル 'bat_gmt_seep.bat' を作成し，組み込みサブルーチン'system(***)'によりこれを実行することです．

FEM解析結果が収められた入力ファイル 'out_fem_zone4m.csv' は，Windowsコマンドプロンプトのコマンドラインから読み込まれます．

出力画像ファイルは以下のとおり．

Filename	Description	Remarks
fig_gmt_mesh.eps	FEMメッシュ図(eps)	ファイル名は固定
fig_gmt_conp.eps	圧力水頭コンター図(eps)	ファイル名は固定
fig_gmt_conh.eps	全水頭コンター図(eps)	ファイル名は固定
fig_gmt_vect.eps	流速ベクトル図(eps)	ファイル名は固定

FEMメッシュ図における境界条件の記号
Symbol	Description
□	全水頭指定境界
○	流量指定境界
△	浸出点境界

025-027

出力ファイル名は固定されているため，copyコマンドにより出力ファイル名を変更

029-033

作業ファイルの削除

プログラムと図面

ソースプログラム

ゾーン型フィルダム

• これは高さ20mのゾーン型フィルダムの解析事例です．
• このケースでは，上流側水深は18m，下流側水深は4mです．
• 四角印の節点では全水頭が与えられ，三角印の節点では浸潤点条件が与えられています．
• なお，ゾーン型フィルダムでは，実際にはコア背面に沿って自由水面が低下し，この解析結果が示すようにコア内での自由水面低下は発生しないことは既知の事項です．

均一型フィルダム

• これは高さ20mの均一型フィルダムの解析事例です．
• このケースでは，上流側水深は18m，下流側水深は4mです．
• 四角印の節点では全水頭が与えられ，三角印の節点では浸潤点条件が与えられています．

トンネル周辺の圧力水頭分布

• これは，全水頭300mを有する上部半円下部矩形断面トンネル周辺の浸透流解析事例です．
• トンネル掘削面での圧力水頭はゼロとしています．

2次元応力解析結果の表示

概要

このプログラムは，GMTによりFEMメッシュ図と主応力コンター図を示すプログラムです．

Fortran プログラム (f90_GMTplnt.f90)

• プログラミングについては，浸透流解析に対するものとほとんど同じです．
• このプログラムでは，コンター作成に用いない要素(材料)を指定することができます．
• 浸透流解析で得られた圧力水頭あるいは全水頭は，連続的に変化しますが，応力解析で得られた応力は，材料特性に依存し連続的には変化しないかもしれません．
• それで全ての応力値がコンター作成に用いられた場合，予期できない結果が得られるかもしれません．
• このため，上記の特別な取り扱いを可能にしました．

下記ファイルはプログラム実行前に準備されていなければなりません．

FEM解析による出力ファイルは，管理人作 'f90_FEM_PLNT.exe' の出力ファイルです．

処理実行用バッチファイルを以下に示します． Fortranプログラムから読み込まれるいくつかのファイルが，このバッチファイルの中で作られています．

001  gfortran -o f90_GMTplnt.exe f90_GMTplnt.f90
002
003  makecpt -Cno_green -T-1/1/0.1 -I > 1_cpt.cpt
004  gawk "{a[NR]=$0}END{for(i=1;i<=NR-1;i++){print a[i]};gsub(/(128)/,\"255\",a[NR]);print a[NR]}" 1_cpt.cpt > 0_cpt.cpt
005
006  rem define of drawing data
007  rem integer::ind,ibc,ifc,isc
008  rem real(4)::baselength,grdspc
009  rem real(4)::bav,bfv
010  rem real(4)::conint,annint
011  rem character::strx*100,stry*100
012  rem integer::ncut,numcut(10)
013  echo 0 1 0 1 > _inp_base.txt
014  echo 12.0 50 >> _inp_base.txt
015  echo 5000.0 5000.0 >> _inp_base.txt
016  echo  0.5 0.5 >> _inp_base.txt
017  echo "x-distance (mm)" >> _inp_base.txt
018  echo "y-distance (mm)" >> _inp_base.txt
019  echo  1 1 >> _inp_base.txt
020
021  rem define of mesh colour (1:MATEL)
022  echo 0,0,0 > _col_mesh.txt
023  echo 255,192,203 >> _col_mesh.txt
024  echo 255,250,205 >> _col_mesh.txt
025
026  f90_GMTplnt out4nodPS0_L.csv
027  copy fig_gmt_mesh.eps fig_gmt_mesh_ps.eps
028  copy fig_gmt_ps1.eps fig_gmt_ps1_ps.eps
029  copy fig_gmt_ps2.eps fig_gmt_ps2_ps.eps
030
031  del fig_gmt_mesh.eps
032  del fig_gmt_ps1.eps
033  del fig_gmt_ps2.eps
034  del _*

001

プログラム 'f90_GMTplnt.f90'のコンパイル

003

GMTコマンド 'makecpt' による，ファイル '1_cpt.cpt' の作成． このケースでは， 'no_green' が基になるファイルとして使われており，レンジは-1から1の範囲，符号はデフォルト値の逆にセットされています．

004

'gawk'により，ファイル '0_cpt.cpt' を作成． ここでの扱いは，ファイル '1_cpt.cpt' の最終行で，色'gray'を'white'に変えています．

013-019

ファイル '_inp_base.txt'の作成． このファイルは，GMT用のパラメータを定義します．

ind	メッシュ図における節点番号・要素番号の表示 (0: 表示しない 1: 表示する)
ibc	メッシュ図における境界条件の表示 (0: 表示しない, 1: 表示する)
ifc	メッシュ図における荷重条件の表示 (0: 表示しない, 1: 表示する)
isc	コンター用スケールマップの表示 (0: 表示しない, 1: 表示する)
baselength	図の寸法(単位：cm)
grdspc	'grd'ファイルにおけるメッシュ間隔
bav	GMTのx軸・y軸における-Bオプションのaに続く値
bfv	GMTのx軸・y軸における-Bオプションのfに続く値Value of 'f' in '-B' option for x and y-axis on GMT
conint	'grdcontour'コマンドのオプション-Cにおけるコンター間隔の値
annint	'grdcontour'コマンドのオプション-Aの目盛間隔の値
strx	x(水平)軸ラベル．文字列のためダブルプライム(")が必要．
stry	y(垂直)軸ラベル．文字列のためダブルプライム(")が必要．
ncut	コンター作成に用いない材料の数．このケースでは鉄管に対する値'1'
numcut(i)	コンター作成に用いない材料種別番号．このケースでは鉄管の残量種別番号に対する'1'

022-024

ファイル '_col_mesh.txt' の作成．このファイルでメッシュ色が定義される．ファイル内の3つの数値は'R G B'を意味する． この事例では，バッチファイル内で'echo'コマンドを使ってファイルが生成され，3種類の材料(MATEL=3)が用いられているため3種類の色が定義される．

026

プログラム 'f90_GMTplnt.exe' の実行．このexeファイルの機能は，バッチファイル 'bat_gmt_plnt.bat' を生成し，組み込みサブルーチン 'system(***)'を用いてこれを実行することです．

FEM解析結果を格納した入力ファイル 'out_4nodPS0_L.csv' はWindowsコマンドプロンプトのコマンドラインから読み込まれます．

出力図は下に示すとおりです．

Filename	Description	Remarks
fig_gmt_mesh.eps	FEMメッシュ図(eps)	ファイル名は固定
fig_gmt_ps1.eps	第一主応力コンター図(eps)	ファイル名は固定
fig_gmt_ps2.eps	第二主応力コンター図(eps)	ファイル名は固定

FEMメッシュ図における境界条件の記号
Symbol	Description
◇	x方向変位指定境界
△	y方向変位指定境界
○	xおよびy方向変位指定境界

027-029

出力ファイル名は固定されているので，copyコマンドを用いてファイル名を変更することが必要．

031-034

作業用ファイルの削除

プログラムと図面

ソースコード

一様引張りを受ける孔を有する板の主応力分布事例

ここでは，一様引張りを受ける孔を有する板の主応力分布の事例が示されています．

埋設水圧鉄管周辺の主応力分布事例

• 内水圧を受ける埋設水圧鉄管の事例が示されています．
• 1mの厚さの充填コンクリートはno-tension材料として扱われています．
• コンクリート要素の第一主応力は小さく，充填コンクリートのno-tension状態が実現されています．
• 水圧鉄管の引張り応力は400MPa程度です．一方，充填コンクリートや岩盤の応力は鉄管応力の1%程度以下です．
• この理由により，全ての材料の応力を用いGMTの'surface'コマンドを用いてコンター図を作成した場合，予期できない結果が得られます．
• それ故に，鉄管部の応力は，コンター図作成のためには参照されていません．

凡例を含む2次元プロット

概要

• このプログラムは，凡例を含む2次元プロット用のバッチファイルを作成し実行するものです．
• 凡例領域の場所は，gnuplotのように変えられます．
• 4種類の線種と4種類の記号を凡例の中で用いることができます．

f90プログラム実行用バッチファイルを下に示します．

gfortran -o f90_GMTgra.exe f90_GMTgra.f90
f90_GMTgra inp_form.txt bat_test.bat

データ読み込みのためのf90プログラムコードを以下に示します．

call getarg(1,fnameR) !input filename for control
call getarg(2,fnameW) !output filename (batch filename for GMT execution)

open(11,file=fnameR,status='old')
read(11,'(a)') fnameF  !output image filename (eps)
read(11,'(a)') xlabel  !label of x-axis
read(11,'(a)') ylabel  !label of y-axis
read(11,'(a)') range   !range of x-axis/range of y-axis
read(11,'(a)') scale   !length of x-axis/length of y-axis
read(11,'(a)') xga     !grid interval and tickmark for x-axis
read(11,'(a)') yga     !grid interval and tickmark for y-axis
read(11,*) nplot,dplot !number of plot data, size of symbol
allocate(fn(1:nplot),sxp(1:nplot),syp(1:nplot),sleg(1:nplot))
allocate(nline(1:nplot),nsymb(1:nplot))
do i=1,nplot
!data filename, column for x, column for y, line type, symbol, title
read(11,*) fn(i),sxp(i),syp(i),nline(i),nsymb(i),sleg(i)
end do
read(11,*) strP,dx,dy   !location of legend, distance from x-axis, distance from y-axis
read(11,*) ipt,scl,llen !font size, font w/h, line length (number of characters)
close(11)

dx, dy

: 凡例領域外枠のx軸・y軸からの距離(最短距離)

ipt

: フォントサイズ (pt)

scl

: 1フォントの幅と高さの比 (width/height)

llen

: 文字数で数えた線の長さ

プログラムと図面

ソースコード

Samples