概要
Pythonで行った平面骨組み構造解析結果より，断面力図を作成するプログラムを作ってみた．
汎用性があるものではなく，必要に応じて書き換えて使うことを前提としている．書き換え箇所は，主として作図範囲と断面力の縮尺程度であろう．色なども凝ってはいないため，その時の気分で書き換えれば良いと思っている．
出力できる図
- 軸力図
- せん断力図
- モーメント図
- 変位モード図
やっていること
- Pythonプログラムにより，matplotlib を用いて断面力図・変位モード図を作成
- 報告書への添付をイメージし，TeX で pdf を作成して，これを ImageMagick により png に変換
出力図事例
断面力図作成プログラム
工夫した点としては，断面力や変位の最大値・最小値をグラフの中に書き込む際，書き込み位置やフォーマットの指定に対し，汎用性を持たせるため，これらの数値を「凡例（legend）」として書き出していることだろうか．
プログラムでは，matplotlibの作図における共通部分の記述を1回とするよう，断面力図毎にforループで作図するようにしている．
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_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画像に変換するもの．
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
以上