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

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

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

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

作例はこちら．
<img src="https://qiita-image-store.s3.amazonaws.com/0/129300/c26abcf6-b225-488c-145c-9336a782f580.png" width=90%>


## 環境

+ 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で行いました．

<table border="1">
<tr><th>項目</th><th>Python作図（前回）</th><th>GMT作図（今回）</th></tr>
<tr><td>Python作図（png）</td><td>242 sec</td><td>  ---  </td></tr>
<tr><td>Pythonデータ処理 </td><td>  ---  </td><td>  2 sec</td></tr>
<tr><td>GMT EPS出力      </td><td>  ---  </td><td>  7 sec</td></tr>
<tr><td>Imagemagick処理  </td><td> 29 sec</td><td> 47 sec</td></tr>
<tr><td>合計             </td><td>271 sec</td><td> 56 sec</td></tr>
</table>

# 全体構成
ここで紹介する 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 は自動的に読み飛ばしてくれるので，心配はない．

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

## 変位モード図作成のための作業ファイル
<table border="1">
<tr><th>変位モード用</th><th>備考</th></tr>
<tr><td>_v0text.txt</td><td>情報表示      </td></tr>
<tr><td>_v0dis.txt </td><td>変位後形状格納</td></tr>
<tr><td>_v0mes.txt </td><td>変位前形状格納</td></tr>
</table>

## 応力ベクトル図作成のための作業ファイル
<table border="1">
<tr><th>第１主応力ベクトル</th><th>第２主応力ベクトル</th><th>備考</th></tr>
<tr><td>_v1text.txt</td><td>_v2text.txt</td><td>情報表示    </td></tr>
<tr><td>_v1ten.txt </td><td>_v2ten.txt </td><td>引張応力格納</td></tr>
<tr><td>_v1com.txt </td><td>_v2com.txt </td><td>圧縮応力格納</td></tr>
</table>

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


# 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）

```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コマンドで，６個のグラフの結合と，リサイズを行っている．

```txt: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

```python: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)
```

以　上