PythonとGMTによる２次元FEM応力解析結果の図化

はじめに

こちら（https://qiita.com/damyarou/items/0e0e68a5f1b96381b490） の記事では，Pythonとmatplotlibにより２次元FEM解析結果の図化プログラムを紹介した．
しかし，このプログラムはとても遅い！それが気になって，Pythonで図化用にデータ処理した後，GMT（Generic Mapping Tools）で図化することを試みた．GMTを使うには，それなりの前処理が必要であるが，GMTのスクリプトそのものは比較的シンプルな記載で思ったことを実行でき，私は個人的には好きである．

作成する図は，matplotlib を用いたものと同じで以下の通り．

• 第１主応力コンター
• 第２主応力コンター
• 最大せん断応力コンター
• 第１主応力ベクトル図
• 第２主応力ベクトル図
• 変位モード図

作例はこちら．

環境

• MacBook Pro (Retina, 13-inch, Mid 2014)
• macOS High Sierra
• Python 3.6.4
• GMT 5.4.3
• ImageMagick 7.0.7

それでスピードは？

スピードの比較結果は以下の通り．GMT の場合，EPS 出力までは早く，ImageMagick による eps から png への変換で少し時間を食うが，トータルでみると GMT を用いた方法のほうが圧倒的に早い．なお時間計測は粗っぽくてもかまわないので，bash の SECONDSで行いました．

項目	Python作図（前回）	GMT作図（今回）
Python作図（png）	242 sec	---
Pythonデータ処理	---	2 sec
GMT EPS出力	---	7 sec
Imagemagick処理	29 sec	47 sec
合計	271 sec	56 sec

全体構成

ここで紹介する GMT による作図では，以下の３つのステップを踏んでいる．

• Python プログラムによる作図のためのデータ作成
• GMT による作図（epsファイル作成）
• ImageMagick による処理（eps の png 化および画像の結合・リサイズ）

上記を実行するスクリプトは以下の通り．SECONDS による時間計測も含めてある．

%%bash

SECONDS=0
python3 py_2dfem_gmt_data.py # Pythonプログラム（作図用データ作成）
time=$SECONDS
echo $time

source a_gmt.txt             # GMT スクリプト（作図）
time=$SECONDS
echo $time

source a_im.txt              # Imagemagickコマンド（画像処理）
time=$SECONDS
echo $time

GMT 作図のためのデータ作成

上記の図を作成するために，以下の作業用ファイルを準備する．たくさんあるが計算機がやってくれるためスペルミスしたり順番を間違えない限りなんの問題もない．これらファイルを作成するプログラムは，Pythonでの作図プログラムから作図部分を取り除き，数値をファイル出力するようにしただけのものである．長いので最後に掲載する．

応力コンター図作成のための作業ファイル

下記着色用ファイル群は，機械的にファイルをオープン・クローズしており，中にはカラのものも存在する．しかしたとえファイルの中身がカラであっても GMT は自動的に読み飛ばしてくれるので，心配はない．

第１主応力コンター	第２主応力コンター	最大せん断応力コンター	備考
_c1text.txt	_c2text.txt	_c3text.txt	情報表示
_c1col00.txt	_c2col00.txt	_c3col00.txt	着色用 ファイル群
_c1col01.txt	_c2col01.txt	_c3col01.txt
_c1col02.txt	_c2col02.txt	_c3col02.txt
_c1col03.txt	_c3col03.txt	_c3col03.txt
_c1col04.txt	_c2col04.txt	_c3col04.txt
_c1col05.txt	_c2col05.txt	_c3col05.txt
_c1col06.txt	_c2col06.txt	_c3col06.txt
_c1col07.txt	_c2col07.txt	_c3col07.txt
_c1col08.txt	_c2col08.txt	_c3col08.txt
_c1col09.txt	_c2col09.txt	_c3col09.txt

変位モード図作成のための作業ファイル

変位モード用	備考
_v0text.txt	情報表示
_v0dis.txt	変位後形状格納
_v0mes.txt	変位前形状格納

応力ベクトル図作成のための作業ファイル

第１主応力ベクトル	第２主応力ベクトル	備考
_v1text.txt	_v2text.txt	情報表示
_v1ten.txt	_v2ten.txt	引張応力格納
_v1com.txt	_v2com.txt	圧縮応力格納

作成される画像ファイル

画像ファイル	説明	備考
fig_gmt_cont1.png	第１主応力コンター	GMTが作成したepsを ImageMagickで pngに変換
fig_gmt_cont2.png	第２主応力コンター
fig_gmt_cont3.png	最大せん断応力コンター
fig_gmt_vect1.png	第１主応力ベクトル
fig_gmt_vect2.png	第２主応力ベクトル
fig_gmt_disp0.png	変位モード
fig_gmt_cv.png	上記６枚の結合画像	ImageMagick による処理
fig_gmt_cvs.png	横幅 1000px のりサイズ画像

GMT作図スクリプト

最初の部分

最初の部分で，グラフ範囲の定義，応力コンター図塗り分け色の定義，GMT psbasemap のパラメータ設定を行っている．

# Range of graph
xmin=-6
xmax=6
dx=1
ymin=175
ymax=185
dy=1

# Color for stress contour
# 省略

# parameters for GMT psbasemap
rangeSW=$xmin/$xmax/$ymin/$ymax
scale=$(($xmax-$xmin))/$(($ymax-$ymin))
afgS=a${dx}f${dx}
afgW=a${dy}f${dy}
labelS='x-direction (m)'
labelW='y-direction (m)'

# confirmation
echo $rangeSW
echo $scale
echo $afgS
echo $afgW
echo $labelS
echo $labelW

psbasemap パラメータを echo で出力させているが，結果は以下の通り．期待通りの値が入っている．
scale=\$((\$xmax-\$xmin))/\$((\$ymax-\$ymin)) は，シェルスクリプトの演算機能（引き算）を使ってみたもの．

-6/6/175/185
12/10
a1f1
a1f1
x-direction (m)
y-direction (m)

共通の psbasemap

gmt psbasemap -R$rangeSW -JX$scale -Bx$afgS+l"$labelS" -By$afgW+l"$labelW" -BWS -P -K > $fig
# -R-6/6/175/185               軸範囲を x: -6 ~ 6，y: 175 ~ 185 に指定
# -JX12/10                     作図領域（グラフ領域）を横12cm，縦10cmに指定
# -Bxa1f1+l"x-direction (m)"   ｘ軸について１ピッチで数値を１ピッチで目盛り線を入れる．
#                              軸ラベルは 'x-direction (m)'
# -Bya1f1+l"y-direction (m)"   ｙ軸について１ピッチで数値を１ピッチで目盛り線を入れる．
#                              軸ラベルは 'y-direction (m)'
# -BWS                         ｙ軸をW（西側＝左側），ｘ軸をS（南側＝下側）に描画

応力コンター

プロット本体

応力コンター（stress distribution contour）では，予め作られた塗り分け色別の10個のファイルからデータを読み取り，色塗りを行っている．

凡例

応力コンターの凡例（legend）では，以下のように，psbasemapで現在のグラフの原点から，原点をｘ方向に0.5cm，ｙ方向に8.5cm移動させた位置に小さな領域（横5cm，縦0.25cm）を定義し，そこにグラフを書くように設定している．
カラーバーのような凡例は，棒グラフを活用している（psxy ... -Sb1u）．高さ１の棒グラフを横並びに10個立て，それを応力コンターで使っている色で塗りつぶしている．

# 凡例用 psbasemap
gmt psbasemap -R0/10/0/1 -JX5/0.25 -Bx -By -X0.5 -Y8.5 -K -O >> $fig
echo '9.5 1' | gmt psxy -R -J -Sb1u -G$col09 -K -O >> $fig
# -R0/10/0/1    グラフ軸範囲指定（x: 0~10, y: 0~1）
# -JX5/0.25     グラフのサイズ（横5cm，縦0.25cm）
# -Bx -By       軸設定は何もしない
# -X0.5 -Y8.5   現在の座標原点（左下隅）よりｘ方向に0.5cm，ｙ方向に8.5cm移動した位置を新しい原点とする
# -Sb1u         ｘ座標 9.5 の位置に高さ 1，横幅 1 の棒グラフを立てる
# -G$col09      指定色（$col09）で棒グラフ内部を塗りつぶす

変位モード図

変位モード図（displacement mode）では，データファイル _v0dis.txt により変位後の形状を薄い黄色で塗りつぶした図を作り，その上にデータファイル _v0mes.txt で変位前のオリジナルのメッシュ図を重ね書きしている．psxy のオプション -L は，閉領域をプロットするときに使うオプションである．この -L オプションには今回気がついた．

gmt psxy _v0dis.txt -R -J -G#ffffcc -W0.5,#aaaaaa -K -O >> $fig
gmt psxy _v0mes.txt -R -J -L -W0.2,#aaaaaa -K -O >> $fig

応力ベクトル図

プロット本体

応力ベクトル図（stress vector）では，まず全要素を黄色で塗りつぶした領域図を作り，その上に青線で示された引張応力ベクトルと赤線で示された圧縮応力ベクトルを上書きしている．

凡例

応力ベクトルの凡例（legend）では，以下のように，psbasemapで現在のグラフの原点から，原点をｘ方向に1.0cm，ｙ方向に8.5cm移動させた位置に小さな領域（横3cm，縦1cm）を定義し，そこにグラフを書くように設定している．グラフといっても，青線と赤線を描画し，その横にテキスト（Tension, Compression）を配置しただけのものである．

# 凡例用 psbasemap
gmt psbasemap -R0/5/0/2 -JX3/1 -Bx -By -X1.0 -Y8.5 -K -O >> $fig
# -R0/5/0/2     グラフ軸範囲指定（x: 0~5, y: 0~2）
# -JX3/1        グラフのサイズ（横3cm，縦1cm）
# -Bx -By       軸設定は何もしない
# -X1.0 -Y8.5   現在の座標原点（左下隅）よりｘ方向に1.0cm，ｙ方向に8.5cm移動した位置を新しい原点とする

GMT スクリプト全文（a_gmt.txt）

a_gmt.txt

# Range of graph
xmin=-6
xmax=6
dx=1
ymin=175
ymax=185
dy=1

# Color for stress contour
col00=#cc0000
col01=#ff0066
col02=#ff9966
col03=#ffcc66
col04=#ffffaa
col05=#00ffff
col06=#00ccff
col07=#0099ff
col08=#0000ff
col09=#000088

# parameters for GMT psbasemap
rangeSW=$xmin/$xmax/$ymin/$ymax
scale=$(($xmax-$xmin))/$(($ymax-$ymin))
afgS=a${dx}f${dx}
afgW=a${dy}f${dy}
labelS='x-direction (m)'
labelW='y-direction (m)'

# confirmation
echo $rangeSW
echo $scale
echo $afgS
echo $afgW
echo $labelS
echo $labelW

gmt set MAP_ANNOT_OFFSET_PRIMARY 0.3c
gmt set MAP_LABEL_OFFSET 0.3c
gmt set MAP_TICK_LENGTH_PRIMARY -0.2c
gmt set MAP_FRAME_PEN thinner,black
gmt set FONT_ANNOT_PRIMARY 10p,Helvetica,black
gmt set FONT_ANNOT_SECONDARY 10p,Helvetica,black
gmt set FONT_LABEL 10p,Helvetica,black

#=================================
# Stress distribution contour
#=================================
for i in 1 2 3
do
fig=fig_gmt_cont${i}.eps
# data plot
gmt psbasemap -R$rangeSW -JX$scale -Bx$afgS+l"$labelS" -By$afgW+l"$labelW" -BWS -P -K > $fig
gmt psxy _c${i}col00.txt -R -J -G$col00 -W0.1,#aaaaaa -K -O >> $fig
gmt psxy _c${i}col01.txt -R -J -G$col01 -W0.1,#aaaaaa -K -O >> $fig
gmt psxy _c${i}col02.txt -R -J -G$col02 -W0.1,#aaaaaa -K -O >> $fig
gmt psxy _c${i}col03.txt -R -J -G$col03 -W0.1,#aaaaaa -K -O >> $fig
gmt psxy _c${i}col04.txt -R -J -G$col04 -W0.1,#aaaaaa -K -O >> $fig
gmt psxy _c${i}col05.txt -R -J -G$col05 -W0.1,#aaaaaa -K -O >> $fig
gmt psxy _c${i}col06.txt -R -J -G$col06 -W0.1,#aaaaaa -K -O >> $fig
gmt psxy _c${i}col07.txt -R -J -G$col07 -W0.1,#aaaaaa -K -O >> $fig
gmt psxy _c${i}col08.txt -R -J -G$col08 -W0.1,#aaaaaa -K -O >> $fig
gmt psxy _c${i}col09.txt -R -J -G$col09 -W0.1,#aaaaaa -K -O >> $fig
gmt pstext _c${i}text.txt -R0/1/0/1 -J -F+f+a+j -N -K -O >> $fig
# legend
gmt psbasemap -R0/10/0/1 -JX5/0.25 -Bx -By -X0.5 -Y8.5 -K -O >> $fig
echo '9.5 1' | gmt psxy -R -J -Sb1u -G$col09 -K -O >> $fig
echo '8.5 1' | gmt psxy -R -J -Sb1u -G$col08 -K -O >> $fig
echo '7.5 1' | gmt psxy -R -J -Sb1u -G$col07 -K -O >> $fig
echo '6.5 1' | gmt psxy -R -J -Sb1u -G$col06 -K -O >> $fig
echo '5.5 1' | gmt psxy -R -J -Sb1u -G$col05 -K -O >> $fig
echo '4.5 1' | gmt psxy -R -J -Sb1u -G$col04 -K -O >> $fig
echo '3.5 1' | gmt psxy -R -J -Sb1u -G$col03 -K -O >> $fig
echo '2.5 1' | gmt psxy -R -J -Sb1u -G$col02 -K -O >> $fig
echo '1.5 1' | gmt psxy -R -J -Sb1u -G$col01 -K -O >> $fig
echo '0.5 1' | gmt psxy -R -J -Sb1u -G$col00 -K -O >> $fig
gmt pstext -R -J -F+f+a+j -N -K -O << EOT >> $fig
1 -0.2  9p,Helvetica-Narrow 0 TC -.8
3 -0.2  9p,Helvetica-Narrow 0 TC -.4
5 -0.2  9p,Helvetica-Narrow 0 TC 0
7 -0.2  9p,Helvetica-Narrow 0 TC .4
9 -0.2  9p,Helvetica-Narrow 0 TC .8
5  1.2  9p,Helvetica 0 BC Comp.  (MPa)  Tens.
EOT
echo '0 0' | gmt psxy -R -J -Sp -O  >> $fig
echo $fig
done

#=================================
# Displacement mode
#=================================
fig=fig_gmt_disp0.eps
gmt psbasemap -R$rangeSW -JX$scale -Bx$afgS+l"$labelS" -By$afgW+l"$labelW" -BWS -P -K > $fig
gmt psxy _v0dis.txt -R -J -G#ffffcc -W0.5,#aaaaaa -K -O >> $fig
gmt psxy _v0mes.txt -R -J -L -W0.2,#aaaaaa -K -O >> $fig
gmt pstext _v0text.txt -R0/1/0/1 -J -F+f+a+j -N -K -O >> $fig
gmt psxy -R -J -Sp -O << EOT  >> $fig
$xmin $ymin
EOT
echo $fig

#=================================
# Stress vector
#=================================
for i in 1 2
do
fig=fig_gmt_vect${i}.eps
# data plot
gmt psbasemap -R$rangeSW -JX$scale -Bx$afgS+l"$labelS" -By$afgW+l"$labelW" -BWS -P -K > $fig
gmt psxy _v0mes.txt -R -J -G#ffffcc -W0.2,#aaaaaa -K -O >> $fig
gmt psxy _v${i}ten.txt -R -JX -Sv0.2+s -W0.3,#0000ff -K -O >> $fig
gmt psxy _v${i}com.txt -R -JX -Sv0.2+s -W0.3,#ff0000 -K -O >> $fig
gmt pstext _v${i}text.txt -R0/1/0/1 -J -F+f+a+j -N -K -O >> $fig
# legend
gmt psbasemap -R0/5/0/2 -JX3/1 -Bx -By -X1.0 -Y8.5 -K -O >> $fig
echo '0 1.5 1.5 1.5' | gmt psxy -R -JX -Sv0.2+s -W1,#0000ff -K -O >> $fig
echo '0 0.5 1.5 0.5' | gmt psxy -R -JX -Sv0.2+s -W1,#ff0000 -K -O >> $fig
gmt pstext -R -J -F+f+a+j -N -K -O << EOT >> $fig
2.0 1.5  9p,Helvetica 0 ML Tension
2.0 0.5  9p,Helvetica 0 ML Compression
EOT
echo '0 0' | gmt psxy -R -J -Sp -O  >> $fig
echo $fig
done

ImageMagickコマンド（a_im.txt）

最後に，Imagemagickコマンドで，６個のグラフの結合と，リサイズを行っている．

a_im.txt

#=================================
# ImageMagick
#=================================
mogrify -trim -density 300 -bordercolor 'transparent' -border 10x10 -format png *.eps
convert -append fig_gmt_cont1.png fig_gmt_cont2.png fig_gmt_cont3.png _1.png
convert -append fig_gmt_vect1.png fig_gmt_vect2.png fig_gmt_disp0.png _2.png
convert +append _1.png _2.png fig_gmt_cv.png
rm _1.png _2.png
rm *.eps

convert -resize 1000x fig_gmt_cv.png fig_gmt_cvs.png

データファイル作成のためのPythonプログラム全文

プログラム名：python:py_2dfem_gmt_data.py

py_2dfem_gmt_data.py

#***************************************
# Data creation for GMT diagrams
#***************************************
import numpy as np
import sys

def DATAINP(nod,fnameR):
    # Data input
    fin=open(fnameR,'r')
    text=fin.readline()
    text=fin.readline()
    text=text.strip()
    text=text.split()
    npoin=int(text[0])
    nele =int(text[1])
    nsec =int(text[2])
    npfix=int(text[3])
    nlod =int(text[4])
    NSTR =int(text[5])

    node=np.zeros([nod+1,nele],dtype=np.int)
    x   =np.zeros(npoin,dtype=np.float64)
    y   =np.zeros(npoin,dtype=np.float64)
    disx=np.zeros(npoin,dtype=np.float64)
    disy=np.zeros(npoin,dtype=np.float64)
    sig1=np.zeros(nele,dtype=np.float64)
    sig2=np.zeros(nele,dtype=np.float64)
    angl=np.zeros(nele,dtype=np.float64)
    taum=np.zeros(nele,dtype=np.float64)

    text=fin.readline()
    for i in range(0,nsec):
        text=fin.readline()

    text=fin.readline()
    for i in range(0,npoin):
        text=fin.readline()
        text=text.strip()
        text=text.split()
        x[i]=float(text[1])
        y[i]=float(text[2])

    text=fin.readline()
    for i in range(0,npfix):
        text=fin.readline()

    text=fin.readline()
    for i in range(0,nele):
        text=fin.readline()
        text=text.strip()
        text=text.split()
        for j in range(0,nod):
            node[j,i]=int(text[j+1]) #node,sectuon number

    text=fin.readline()
    for i in range(0,npoin):
        text=fin.readline()
        text=text.strip()
        text=text.split()
        disx[i]=float(text[1]) # displacement in x-direction
        disy[i]=float(text[2]) # displacement in y-direction

    text=fin.readline()
    for i in range(0,nele):
        text=fin.readline()
        text=text.strip()
        text=text.split()
        sig1[i]=float(text[4])
        sig2[i]=float(text[5])
        angl[i]=float(text[6])
        taum[i]=0.5*abs(sig1[i]-sig2[i])
    fin.close()
    return npoin,nele,node,x,y,disx,disy,sig1,sig2,angl,taum


def SIGCOL(sig,sg):
    if sg[8]<sig: icol=9
    if sg[7]<sig  and  sg[8]>=sig: icol=8
    if sg[6]<sig  and  sg[7]>=sig: icol=7
    if sg[5]<sig  and  sg[6]>=sig: icol=6
    if sg[4]<sig  and  sg[5]>=sig: icol=5
    if sg[4]>=sig and  sg[3]<=sig: icol=4
    if sg[3]>sig  and  sg[2]<=sig: icol=3
    if sg[2]>sig  and  sg[1]<=sig: icol=2
    if sg[1]>sig  and  sg[0]<=sig: icol=1
    if sg[0]>sig: icol=0
    return icol

def CON(nod,nnn,_sig,x,y,node,nele,ss):
    _s=[]
    for i in range(1,len(ss)-1):
        _s=_s+[float(ss[i])]
    sg=np.array(_s)
    if nnn==0:
        # first principal stress
        ls0='* First principal stress'
        ls1='sig1_max={0:8.3f} MPa'.format(np.max(_sig))
        ls2='sig1_min={0:8.3f} MPa'.format(np.min(_sig))
        ft=open('_c1text.txt','w')
        print('0.95 0.95  10p,Courier 0 TR '+ls0,file=ft)
        print('0.95 0.90  10p,Courier 0 TR '+ls1,file=ft)
        print('0.95 0.85  10p,Courier 0 TR '+ls2,file=ft)
        ft.close()
    if nnn==1:
        # second principal stress
        ls0='* Second principal stress'
        ls1='sig2_max={0:8.3f} MPa'.format(np.max(_sig))
        ls2='sig2_min={0:8.3f} MPa'.format(np.min(_sig))
        ft=open('_c2text.txt','w')
        print('0.95 0.95  10p,Courier 0 TR '+ls0,file=ft)
        print('0.95 0.90  10p,Courier 0 TR '+ls1,file=ft)
        print('0.95 0.85  10p,Courier 0 TR '+ls2,file=ft)
        ft.close()
    if nnn==2:
        # maximum shearing stress
        ls0='* Maximum shearing stress'
        ls1='tau_max={0:8.3f} MPa'.format(np.max(_sig))
        ls2='tau_min={0:8.3f} MPa'.format(np.min(_sig))
        ft=open('_c3text.txt','w')
        print('0.95 0.95  10p,Courier 0 TR '+ls0,file=ft)
        print('0.95 0.90  10p,Courier 0 TR '+ls1,file=ft)
        print('0.95 0.85  10p,Courier 0 TR '+ls2,file=ft)
        ft.close()
    f0=open('_c'+str(nnn+1)+'col00.txt','w')
    f1=open('_c'+str(nnn+1)+'col01.txt','w')
    f2=open('_c'+str(nnn+1)+'col02.txt','w')
    f3=open('_c'+str(nnn+1)+'col03.txt','w')
    f4=open('_c'+str(nnn+1)+'col04.txt','w')
    f5=open('_c'+str(nnn+1)+'col05.txt','w')
    f6=open('_c'+str(nnn+1)+'col06.txt','w')
    f7=open('_c'+str(nnn+1)+'col07.txt','w')
    f8=open('_c'+str(nnn+1)+'col08.txt','w')
    f9=open('_c'+str(nnn+1)+'col09.txt','w')
    for ne in range(0,nele):
        sig=_sig[ne]
        icol=SIGCOL(sig,sg)
        nn=node[0:nod,ne]-1
        xa=x[nn]
        ya=y[nn]
        if icol==0: WR(nod,xa,ya,f0)
        if icol==1: WR(nod,xa,ya,f1)
        if icol==2: WR(nod,xa,ya,f2)
        if icol==3: WR(nod,xa,ya,f3)
        if icol==4: WR(nod,xa,ya,f4)
        if icol==5: WR(nod,xa,ya,f5)
        if icol==6: WR(nod,xa,ya,f6)
        if icol==7: WR(nod,xa,ya,f7)
        if icol==8: WR(nod,xa,ya,f8)
        if icol==9: WR(nod,xa,ya,f9)
    f0.close()
    f1.close()
    f2.close()
    f3.close()
    f4.close()
    f5.close()
    f6.close()
    f7.close()
    f8.close()
    f9.close()
    print(ls0)


def DIS(nod,dmax,disx,disy,x,y,node,nele):
    scl_dis=0.5
    # displacement
    ls0='* Displacement mode'
    ls1='disx_max={0:8.3f} mm'.format(np.max(disx))
    ls2='disx_min={0:8.3f} mm'.format(np.min(disx))
    ls3='disy_max={0:8.3f} mm'.format(np.max(disy))
    ls4='disy_min={0:8.3f} mm'.format(np.min(disy))
    ft=open('_v0text.txt','w')
    print('0.95 0.95  10p,Courier 0 TR '+ls0,file=ft)
    print('0.95 0.90  10p,Courier 0 TR '+ls1,file=ft)
    print('0.95 0.85  10p,Courier 0 TR '+ls2,file=ft)
    print('0.95 0.80  10p,Courier 0 TR '+ls3,file=ft)
    print('0.95 0.75  10p,Courier 0 TR '+ls4,file=ft)
    ft.close()
    rx=x+disx/dmax*scl_dis
    ry=y+disy/dmax*scl_dis
    ff=open('_v0dis.txt','w')
    for ne in range(0,nele):
        nn=node[0:nod,ne]-1
        xa=rx[nn]
        ya=ry[nn]
        WR(nod,xa,ya,ff)
    ff.close()
    ff=open('_v0mes.txt','w')
    for ne in range(0,nele):
        nn=node[0:nod,ne]-1
        xa=x[nn]
        ya=y[nn]
        WR(nod,xa,ya,ff)
    ff.close()
    print(ls0)

def VEC(nod,nnn,vmax,_sig,angl,x,y,node,nele):
    mmm=nnn-3
    scl_vec=1.5
    if mmm==1:
        # first principal stress
        _ang=angl
        ls0='* First principal stress'
        ls1='sig1_max={0:8.3f} MPa'.format(np.max(_sig))
        ls2='sig1_min={0:8.3f} MPa'.format(np.min(_sig))
    if mmm==2:
        # second principal stress
        _ang=angl+90.0
        ls0='* Second principal stress'
        ls1='sig2_max={0:8.3f} MPa'.format(np.max(_sig))
        ls2='sig2_min={0:8.3f} MPa'.format(np.min(_sig))
    ft=open('_v'+str(mmm)+'text.txt','w')
    print('0.95 0.95  10p,Courier 0 TR '+ls0,file=ft)
    print('0.95 0.90  10p,Courier 0 TR '+ls1,file=ft)
    print('0.95 0.85  10p,Courier 0 TR '+ls2,file=ft)
    ft.close()
    xg=np.zeros(nele)
    yg=np.zeros(nele)
    for ne in range(0,nele):
        nn=node[0:nod,ne]-1
        xa=x[nn]
        ya=y[nn]
        xg[ne]=np.average(xa)
        yg[ne]=np.average(ya)
    al=abs(_sig)/vmax*scl_vec
    x1=xg-0.5*al*np.cos(np.radians(_ang))
    y1=yg-0.5*al*np.sin(np.radians(_ang))
    x2=xg+0.5*al*np.cos(np.radians(_ang))
    y2=yg+0.5*al*np.sin(np.radians(_ang))
    f0=open('_v'+str(mmm)+'ten.txt','w')
    f1=open('_v'+str(mmm)+'com.txt','w')
    for ne in range(0,nele):
        if _sig[ne]<0.0: # compression
            print('{0:10.3f} {1:10.3f} {2:10.3f} {3:10.3f}'.format(x1[ne],y1[ne],x2[ne],y2[ne]),file=f1)
        else: # tension
            print('{0:10.3f} {1:10.3f} {2:10.3f} {3:10.3f}'.format(x1[ne],y1[ne],x2[ne],y2[ne]),file=f0)
    f0.close()
    f1.close()
    print(ls0)


def WR(nod,xa,ya,ff):
    print('>',file=ff)
    for j in range(0,nod):
        print('{0:10.3f} {1:10.3f}'.format(xa[j],ya[j]),file=ff)


def INP3():
    nod=3
    fnameR='out_arch3.txt'
    xmin=-7 ; xmax=7  ; dx=1
    ymin=64 ; ymax=77 ; dy=1
    ss=['','-10','-6',-4,'-2','0','1','2','3','5','']
    return nod,fnameR,xmin,xmax,dx,ymin,ymax,dy,ss


def INP4():
    nod=4
    fnameR='out_arch4.txt'
    # range of drawing
    xmin=-6 ; xmax=6  ; dx=1
    ymin=175; ymax=185; dy=1
    ss=['','-.8','-.6','-.4','-.2','0','.2','.4','.6','.8','']
    return nod,fnameR,xmin,xmax,dx,ymin,ymax,dy,ss


#===================
# Main routine
#===================
#nod,fnameR,xmin,xmax,dx,ymin,ymax,dy,ss=INP3()
nod,fnameR,xmin,xmax,dx,ymin,ymax,dy,ss=INP4()
npoin,nele,node,x,y,disx,disy,sig1,sig2,angl,taum=DATAINP(nod,fnameR)

print('xmin, xmax ={0:10.3f} {1:10.3f}'.format(np.min(x),np.max(x)))
print('ymin, ymax ={0:10.3f} {1:10.3f}'.format(np.min(y),np.max(y)))
#stry=input('Continue? y/n ')
#if stry=='n': sys.exit()

# t/m2 -> MPa
sig1=sig1*0.01
sig2=sig2*0.01
taum=taum*0.01
# m -> mm
disx=disx*1000.0
disy=disy*1000.0

dmax=np.max([np.max(np.abs(disx)),np.max(np.abs(disy))])
vec1_max=max(abs(np.max(sig1)),abs(np.min(sig1)))
vec2_max=max(abs(np.max(sig2)),abs(np.min(sig2)))
vmax=np.max([vec1_max,vec2_max])

for nnn in range(0,6):
    if nnn==0: CON(nod,nnn,sig1,x,y,node,nele,ss)
    if nnn==1: CON(nod,nnn,sig2,x,y,node,nele,ss)
    if nnn==2: CON(nod,nnn,taum,x,y,node,nele,ss)
    if nnn==3: DIS(nod,dmax,disx,disy,x,y,node,nele)
    if nnn==4: VEC(nod,nnn,vmax,sig1,angl,x,y,node,nele)
    if nnn==5: VEC(nod,nnn,vmax,sig2,angl,x,y,node,nele)

以　上