Pythonで平面骨組解析 (3) 断面力図作成

## 概要
Pythonで行った平面骨組み構造解析結果より，断面力図を作成するプログラムを作ってみた．
汎用性があるものではなく，必要に応じて書き換えて使うことを前提としている．書き換え箇所は，主として作図範囲と断面力の縮尺程度であろう．色なども凝ってはいないため，その時の気分で書き換えれば良いと思っている．

### 出力できる図

+ 軸力図
+ せん断力図
+ モーメント図
+ 変位モード図

### やっていること

+ Pythonプログラムにより，matplotlib を用いて断面力図・変位モード図を作成
+ 報告書への添付をイメージし，TeX で pdf を作成して，これを ImageMagick により png に変換

## 出力図事例
![tex_fig.png](https://qiita-image-store.s3.amazonaws.com/0/129300/1da30632-dff2-21d8-d3cc-5725ce1e3006.png)


## 断面力図作成プログラム
工夫した点としては，断面力や変位の最大値・最小値をグラフの中に書き込む際，書き込み位置やフォーマットの指定に対し，汎用性を持たせるため，これらの数値を「凡例（legend）」として書き出していることだろうか．

プログラムでは，matplotlibの作図における共通部分の記述を1回とするよう，断面力図毎にforループで作図するようにしている．

```python:py_force.py
import matplotlib.pyplot as plt
import numpy as np
import sys

def calc(ne,node,x,y,d1,d2):
    i=node[0,ne]-1
    j=node[1,ne]-1
    x1=x[i]
    y1=y[i]
    x2=x[j]
    y2=y[j]
    al=np.sqrt((x2-x1)**2+(y2-y1)**2)
    theta=np.arccos((x2-x1)/al)
    x4=x1-d1[ne]*np.sin(theta)
    y4=y1+d1[ne]*np.cos(theta)
    x3=x2-d2[ne]*np.sin(theta)
    y3=y2+d2[ne]*np.cos(theta)
    return x1,x2,x3,x4,y1,y2,y3,y4
    

# Main routine
args = sys.argv
fnameR=args[1] # input data file

f=open(fnameR,'r')
text=f.readline()
text=f.readline()
text=text.strip()
text=text.split()
npoin=int(text[0]) # Number of nodes
nele =int(text[1]) # Number of elements
nsec =int(text[2]) # Number of sections
npfix=int(text[3]) # Number of restricted nodes
nlod =int(text[4]) # Number of loaded nodes

x   =np.zeros(npoin,dtype=np.float64) # Coordinates of nodes
y   =np.zeros(npoin,dtype=np.float64) # Coordinates of nodes
node=np.zeros([2,nele],dtype=np.int)  # Node-element relationship
disx=np.zeros(npoin,dtype=np.float64) # Coordinates of nodes
disy=np.zeros(npoin,dtype=np.float64) # Coordinates of nodes
N1  =np.zeros(nele,dtype=np.float64)  # Section force vector
S1  =np.zeros(nele,dtype=np.float64)  # Section force vector
M1  =np.zeros(nele,dtype=np.float64)  # Section force vector
N2  =np.zeros(nele,dtype=np.float64)  # Section force vector
S2  =np.zeros(nele,dtype=np.float64)  # Section force vector
M2  =np.zeros(nele,dtype=np.float64)  # Section force vector

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

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

text=f.readline()
for i in range(0,nele):
    text=f.readline()
    text=text.strip()
    text=text.split()
    node[0,i]=int(text[1]) #node_1
    node[1,i]=int(text[2]) #node_2

text=f.readline()
for i in range(0,npoin):
    text=f.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=f.readline()
for i in range(0,nele):
    text=f.readline()
    text=text.strip()
    text=text.split()
    N1[i]=-float(text[1]) # axial force at node-1
    S1[i]= float(text[2]) # shear force at node-1
    M1[i]= float(text[3]) # moment at node-1
    N2[i]= float(text[4]) # axial force at node-2
    S2[i]=-float(text[5]) # shear force at node-2
    M2[i]=-float(text[6]) # moment at node-2
f.close()

nmax=np.max([np.max(np.abs(N1)),np.max(np.abs(N2))])
smax=np.max([np.max(np.abs(S1)),np.max(np.abs(S2))])
mmax=np.max([np.max(np.abs(M1)),np.max(np.abs(M2))])
dmax=np.max([np.max(np.abs(disx)),np.max(np.abs(disy))])

xmin=-3
xmax=13
ymin=-3
ymax=9

scl_dis=1.0
scl_axi=1.0
scl_she=1.0
scl_mom=2.0

for nnn in range(0,4):
    ax=plt.subplot(111)
    ax.set_xlim([xmin,xmax])
    ax.set_ylim([ymin,ymax])
    ax.set_xlabel('x-direction (m)')
    ax.set_ylabel('y-direction (m)')
    ax.spines['right'].set_visible(False)
    ax.spines['top'].set_visible(False)
    ax.yaxis.set_ticks_position('left')
    ax.xaxis.set_ticks_position('bottom')
    aspect = (ymax-ymin)/(xmax-xmin)*(ax.get_xlim()[1] - ax.get_xlim()[0]) / (ax.get_ylim()[1] - ax.get_ylim()[0])
    ax.set_aspect(aspect)

    if nnn==0:
        # displacement
        fnameF='fig_dis.png'
        ls1='disx_max={0:15.7e}'.format(np.max(disx))
        ls2='disx_min={0:15.7e}'.format(np.min(disx))
        ls3='disy_max={0:15.7e}'.format(np.max(disy))
        ls4='disy_min={0:15.7e}'.format(np.min(disy))
        dx=x+disx/dmax*scl_dis
        dy=y+disy/dmax*scl_dis
        for ne in range(0,nele):
            n1=node[0,ne]-1
            n2=node[1,ne]-1
            ax.plot([x[n1],x[n2]],[y[n1],y[n2]],color='gray',linewidth=0.5)
            ax.plot([dx[n1],dx[n2]],[dy[n1],dy[n2]],color='black',linewidth=1)
    if nnn==1:
        # axial force diagram
        fnameF='fig_axi.png'
        ls1='N_max={0:15.7e}'.format(np.max([np.max(N1),np.max(N2)]))
        ls2='N_min={0:15.7e}'.format(np.min([np.min(N1),np.min(N2)]))
        ls3=''
        ls4=''
        d1=N1/nmax*scl_axi
        d2=N2/nmax*scl_axi
        for ne in range(0,nele):
            x1,x2,x3,x4,y1,y2,y3,y4=calc(ne,node,x,y,d1,d2)
            if d1[ne]<=0.0: # compression
                ax.fill([x1,x2,x3,x4],[y1,y2,y3,y4],color='black',alpha=0.1)
            else: # tension
                ax.fill([x1,x2,x3,x4],[y1,y2,y3,y4],color='black',alpha=0.2)
        for ne in range(0,nele):
            n1=node[0,ne]-1
            n2=node[1,ne]-1
            plt.plot([x[n1],x[n2]],[y[n1],y[n2]],color='black',linewidth=0.5)
    if nnn==2:
        # shearing force
        fnameF='fig_she.png'
        ls1='S_max={0:15.7e}'.format(np.max([np.max(-S1),np.max(-S2)]))
        ls2='S_min={0:15.7e}'.format(np.min([np.min(-S1),np.min(-S2)]))
        ls3=''
        ls4=''
        d1=S1/smax*scl_she
        d2=S2/smax*scl_she
        for ne in range(0,nele):
            x1,x2,x3,x4,y1,y2,y3,y4=calc(ne,node,x,y,d1,d2)
            ax.fill([x1,x2,x3,x4],[y1,y2,y3,y4],color='black',alpha=0.1)
        for ne in range(0,nele):
            n1=node[0,ne]-1
            n2=node[1,ne]-1
            ax.plot([x[n1],x[n2]],[y[n1],y[n2]],color='black',linewidth=0.5)
    if nnn==3:
        # moment
        fnameF='fig_mom.png'
        ls1='M_max={0:15.7e}'.format(np.max([np.max(-M1),np.max(-M2)]))
        ls2='M_min={0:15.7e}'.format(np.min([np.min(-M1),np.min(-M2)]))
        ls3=''
        ls4=''
        d1=M1/mmax*scl_mom
        d2=M2/mmax*scl_mom
        for ne in range(0,nele):
            x1,x2,x3,x4,y1,y2,y3,y4=calc(ne,node,x,y,d1,d2)
            ax.fill([x1,x2,x3,x4],[y1,y2,y3,y4],color='black',alpha=0.1)
        for ne in range(0,nele):
            n1=node[0,ne]-1
            n2=node[1,ne]-1
            ax.plot([x[n1],x[n2]],[y[n1],y[n2]],color='black',linewidth=0.5)

    ax.plot(xmin,ymin,'.',label=ls1)
    ax.plot(xmin,ymin,'.',label=ls2)
    ax.plot(xmin,ymin,'.',label=ls3)
    ax.plot(xmin,ymin,'.',label=ls4)
    ax.legend(loc='upper right',numpoints=1,markerscale=0, frameon=False,prop={'family':'monospace','size':12})
    plt.savefig(fnameF, bbox_inches="tight", pad_inches=0.2)
    plt.clf()
```

## TeX コマンド
TeXにより，4枚のグラフをA3横1枚にアレンジしている．

```tex:tex_fig_tex
\documentclass[english]{jsarticle}
\usepackage[a3paper,landscape,top=25mm,bottom=25mm,left=25mm,right=25mm]{geometry}
\usepackage[dvipdfmx]{graphicx}
\pagestyle{empty}

\begin{document}

\begin{center}
\begin{tabular}{|c|c|}\hline

\begin{minipage}{14.0cm}\vspace{0.2zh}\includegraphics[width=14.0cm,bb={0 0 715 568}]{fig_axi.png}\end{minipage}&
\begin{minipage}{14.0cm}\vspace{0.2zh}\includegraphics[width=14.0cm,bb={0 0 715 568}]{fig_mom.png}\end{minipage}\\
\LARGE \textsf{Axial force} & \LARGE \textsf{Moment} \\ \hline
\begin{minipage}{14.0cm}\vspace{0.2zh}\includegraphics[width=14.0cm,bb={0 0 715 568}]{fig_she.png}\end{minipage}&
\begin{minipage}{14.0cm}\vspace{0.2zh}\includegraphics[width=14.0cm,bb={0 0 715 568}]{fig_dis.png}\end{minipage}\\
\LARGE \textsf{Shearing force} & \LARGE \textsf{Displacement mode} \\ \hline
\end{tabular}
\end{center}

\centerline{\LARGE \textsf{Fig Section Force Diagrams}}

\end{document}
```

## TeX コマンド実行スクリプト
convert は ImageMagick のコマンドであり，pdfから余白を取り除き，png画像に変換するもの．

```text:a_tex.txt
platex tex_fig.tex
dvipdfmx -p a3 tex_fig.dvi

convert -trim -density 400 tex_fig.pdf -bordercolor 'transparent' -border 20x20 -quality 100 tex_fig.png
```

以上