はじめに

２次元FEM応力解析を行ったあと，やはりすぐに結果を見てみたいものである．
このため，昔から図化プログラムを作っており，いつか解説を作ろうと思っていたのだが，めんどうくさくてなかなかできなかった．
今回は，少し時間ができたので，それを試みたものである．

やっていること

私の場合は，２次元FEM応力解析結果に対し，Python matplotlibにより，以下の６個の出力図を作っている．

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

上記のうち，コンターといっても，要素を，その要素の平均応力の値に応じて色分けしただけである．また，最大せん断応力は，第１主応力と第２主応力から定まるモール円の半径と定義している．

解析結果は，以下の記事で紹介したプログラムによるものを対象としている．

FEM解析結果の出力フォーマットが変われば，図化プログラムの入力部分も変わってくるが，そこらへんは適宜変更ということで．

上述６個の出力図を作ったあと，関連性をチェックしやすくするため，ImageMagickにより１枚の図にまとめている．
まとめ方は以下の通り．画像を作るのが面倒なので，表として書いているがご容赦を．

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

作例はこちら．

環境

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

このプログラムの問題点，それはとても遅いことである．

GMTで作図したら，早くなる．
GMTでの作図例（https://qiita.com/damyarou/items/23607fb0a252b53ebebb）

作図プログラム解説

私は，個人的にpyplotを直接呼び込むのが好きなので，そのようにプログラムしている．
以下，プログラム全文を区切って説明する．少し長いがご勘弁を．

モジュール読み込み

sysは，コマンドライン入力を使うためにインポートしているが，ここでは使っていない．
フォント指定を行うために，FontProperties を読み込んでいる．桁数を揃えて見栄えを良くするため，一部等幅フォントを使っている．

#***************************************
# Diagrams for 2D FEM result
#***************************************
import matplotlib.pyplot as plt
import numpy as np
import sys
from matplotlib.font_manager import FontProperties

データ読み込み

この部分は，各人の持っているFEM解析プログラムに合わせる必要あり．
ここでは，関数の引数に１要素の節点数を示す変数 nod を入れており，三角形要素・四角形要素双方に対応させている．

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

応力コンター作図（要素の色塗り）

関数 SIGCOL

ここでは応力値を10段階に分けて色塗りするようにしている．
しきい値配列 sg は，解析モデルやケースにより異なるため，あとで出てくる「入力ファイル指定」の際，指定するようにしている．

関数 CON

要素に色塗りを行う関数．
塗り分ける色は10色としている．圧縮側（マイナス応力）は茶色から薄い黄色，引張側（プラス応力）は水色から濃い青として色を指定した．完全に私の趣味である．
色塗りは，plt.fill(xa,ya,facecolor=col[icol],edgecolor='#666666',lw=0.1) のように，facecolor で応力値に応じた色で塗りつぶし，edgecolor で要素境界を灰色で描くようにしている．
応力の最大・最小を示すテキストは，plt.text(0.9,0.9,sinfo,va='top',ha='right',fontproperties=fp,transform=plt.gca().transAxes) の部分で作図領域の右上に表示している．フォントは fontproperties=fp で指定し（ここでは等幅のMenlo.ttc），transform=plt.gca().transAxes でグラフ領域の相対座標指定（左下[0,0]，右上[1,1]）している．
凡例のカラーバーは自作した．グラフ領域の中に plt.axes([xcb,ycb,wcb,hcb]) で小さな領域を作り，そこを plt.axvspan で区切って塗りつぶしている．

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,node,x,y,nele,ss,fp):
    _s=[]
    for i in range(1,len(ss)-1):
        _s=_s+[float(ss[i])]
    sg=np.array(_s)
    col=['#cc0000','#ff0066','#ff9966','#ffcc66','#ffffaa','#00ffff','#00ccff','#0099ff','#0000ff','#000088']
    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))
    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))
    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))
    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]
        plt.fill(xa,ya,facecolor=col[icol],edgecolor='#666666',lw=0.1)
    print(ls0)
    sinfo=ls0+'\n'+ls1+'\n'+ls2
    plt.text(0.9,0.9,sinfo,va='top',ha='right',fontproperties=fp,transform=plt.gca().transAxes)
    # Legend (color bar)
    xcb=0.2; ycb=0.8; wcb=0.2; hcb=0.02
    plt.axes([xcb,ycb,wcb,hcb])
    plt.xlim([0,10])
    plt.ylim([0,1])
    plt.axis('off')
    for i in range(0,10):
        x1=float(i)
        x2=float(i+1)
        plt.axvspan(x1,x2, facecolor=col[i], alpha=1)
    for i in range(0,11):
        plt.text(float(i),-0.2,ss[i],fontsize=12,va='top',ha='center')
    plt.text(5,1.1,'Comp.  (MPa)  Tens.',fontsize=10,va='bottom',ha='center')

変位モード作図

変位の最大値が，図上で 0.5m （変数 scl_dis=0.5）となるように変位を変換している．
変位後の形状は plt.fill(xa,ya,facecolor='#ffffcc',edgecolor='#aaaaaa',lw=0.5) で要素を薄い黄色で，境界を薄い灰色で着色している．
変位前の形状は plt.fill(xa,ya,facecolor='none',edgecolor='#aaaaaa',lw=0.2) で要素境界のみを薄い灰色で着色している．facecolor='none' になかなか気が付かなかった．

def DIS(nod,dmax,disx,disy,x,y,node,nele,fp):
    # displacement
    scl_dis=0.5
    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))
    rx=x+disx/dmax*scl_dis
    ry=y+disy/dmax*scl_dis
    for ne in range(0,nele):
        nn=node[0:nod,ne]-1
        xa=rx[nn]
        ya=ry[nn]
        plt.fill(xa,ya,facecolor='#ffffcc',edgecolor='#aaaaaa',lw=0.5)
    for ne in range(0,nele):
        nn=node[0:nod,ne]-1
        xa=x[nn]
        ya=y[nn]
        plt.fill(xa,ya,facecolor='none',edgecolor='#aaaaaa',lw=0.2)
    print(ls0)
    sinfo=ls0+'\n'+ls1+'\n'+ls2+'\n'+ls3+'\n'+ls4
    plt.text(0.98,0.98,sinfo,va='top',ha='right',fontproperties=fp,transform=plt.gca().transAxes)

応力ベクトル作図

矢印を書くのは面倒なので，引張応力は青線，圧縮応力は赤線で表示するようにした．
応力の最大値が図上で 1.5m （変数 scl_vec=1.5）となるように線の長さを変換している．
凡例は，引張応力を示す青線と，圧縮応力を示す赤線を長さゼロのダミーで描き，凡例化し左上に表示している．

def VEC(nod,nnn,vmax,_sig,angl,x,y,node,nele,fp):
    scl_vec=1.5
    if nnn==4:
        # 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 nnn==5:
        # 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))
    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]
        plt.fill(xa,ya,facecolor='#ffffcc',edgecolor='#aaaaaa',lw=0.2)
        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))
    col=['#ff0000','#0000ff']
    for ne in range(0,nele):
        if _sig[ne]<0.0:
            icol=0
        else:
            icol=1
        plt.plot([x1[ne],x2[ne]],[y1[ne],y2[ne]],color=col[icol],lw=0.5)
    print(ls0)
    sinfo=ls0+'\n'+ls1+'\n'+ls2
    plt.text(0.98,0.98,sinfo,va='top',ha='right',fontproperties=fp,transform=plt.gca().transAxes)
    # Legend
    plt.plot([xmin],[ymin],'-',color=col[1],lw=1.5,label='Tension')
    plt.plot([xmin],[ymin],'-',color=col[0],lw=1.5,label='Compression')
    plt.legend(loc='upper left')

入力ファイル指定

３節点要素用と４節点要素用の２つの関数を準備した．作図事例は４節点要素のものである．ここで指定しているものは以下の通り．

１要素節点数
入力ファイル名（FEM解析結果ファイル）
ｘ方向作図範囲
ｙ方向作図範囲
応力区分（応力コンターを作る上での応力区分）

応力区分は，応力コンターを作る上での区分であり，等間隔でなくても構わない．
入力は，カラーバーの凡例に書き込むため，文字列入力としている．
９個の応力値に加え，最初と最後は「カラ」のリストである．
これは，凡例表示したカラーバーを見ればわかると思うが，カラーバー両端の数値は不要であることを意味する．
すなわち，def INP4() の事例でいえば，圧縮側 -0.8MPa より小さな応力は茶色，引張側 +0.8MPa より大きな応力は濃い青で表示されており，カラーバーの両端に数値は不要であるため．

もちろん，この部分は，入力ファイルや作図範囲，応力コンターの区切り値を変えたい場合には書き換える必要がある．

def INP3():                   # 関数名（３節点要素用）
    nod=3                     # １要素節点数
    fnameR='out_arch3.txt'    # 入力ファイル名（FEM解析結果ファイル）
    xmin=-7 ; xmax=7  ; dx=1  # x方向作図領域
    ymin=64 ; ymax=77 ; dy=1  # y方向作図領域
    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'    # 入力ファイル名（FEM解析結果ファイル）
    xmin=-6 ; xmax=6  ; dx=1  # x方向作図領域
    ymin=175; ymax=185; dy=1  # y方向作図領域
    ss=['','-.8','-.6','-.4','-.2','0','.2','.4','.6','.8',''] # 応力区分
    return nod,fnameR,xmin,xmax,dx,ymin,ymax,dy,ss

メインルーチン（データ読み込みと作図準備）

データを読み込んだあと，作図のためのいくつかの処理を行う．
ここでは，細かいことは言わず，1MPa = 100t/m2 である．
vmax（応力の最大値）と dmax（変位の最大値）は，作図上，応力ベクトル線長あるいは変位を制御するために用いている．

#===================
# 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
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])
# m -> mm
disx=disx*1000.0
disy=disy*1000.0
dmax=np.max([np.max(np.abs(disx)),np.max(np.abs(disy))])

メインルーチン（作図部分）

等幅フォントで表示したい部分に対し，fp でフォントとそのサイズを指定している．
出力画像ファイル名はリストで指定し，順次作図・保存していく．

fsz=14
fp = FontProperties(fname=r'/System/Library/Fonts/Menlo.ttc', size=fsz)
if nod==3: flist=['_fig_pl3_con1.png','_fig_pl3_con2.png','_fig_pl3_con3.png',
                  '_fig_pl3_disp.png','_fig_pl3_vec1.png','_fig_pl3_vec2.png']
if nod==4: flist=['_fig_pl4_con1.png','_fig_pl4_con2.png','_fig_pl4_con3.png',
                  '_fig_pl4_disp.png','_fig_pl4_vec1.png','_fig_pl4_vec2.png']
for nnn,fnameF in enumerate(flist):
    ixx=10
    fig = plt.figure(figsize=(ixx,ixx/(xmax-xmin)*(ymax-ymin)),facecolor='w')
    plt.rcParams['font.family'] = 'Arial Narrow'
    plt.rcParams["font.size"] = fsz
    plt.xlim([xmin,xmax])
    plt.ylim([ymin,ymax])
    plt.xticks(np.arange(xmin,xmax+dx,dx))
    plt.yticks(np.arange(ymin,ymax+dy,dy))
    plt.xlabel('x-direction (m)')
    plt.ylabel('y-direction (m)')
    plt.gca().spines['right'].set_visible(False)
    plt.gca().spines['top'].set_visible(False)
    plt.gca().yaxis.set_ticks_position('left')
    plt.gca().xaxis.set_ticks_position('bottom')
    plt.gca().set_aspect('equal', adjustable='box')

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

    plt.savefig(fnameF, dpi=300, bbox_inches="tight", pad_inches=0.2)
    plt.clf()

ImageMagickによる処理

Jupyter notebook を使っているので，シェルスクリプトを実行するため，先頭に %%bash を記述している．
やっていることは，以下の通り．

３個のコンター図を縦に結合
２個のベクトル図と変位モード図を縦に結合
これら２つの縦長画像を横に結合
最後に扱い良さを考慮して，横幅 1000px にリサイズしている．

%%bash

convert -append _fig_pl4_con1.png _fig_pl4_con2.png _fig_pl4_con3.png _1.png
convert -append _fig_pl4_vec1.png _fig_pl4_vec2.png _fig_pl4_disp.png _2.png
convert +append _1.png _2.png fig_pl4_cv.png
rm _1.png _2.png

convert -resize 1000x fig_pl4_cv.png fig_pl4_cvs.png

以　上