Plane skeleton analysis with Python (3) Creation of section force diagram
Overview
From the results of the plane frame structure analysis performed with Python, I made a program to create a section force diagram. It is not versatile, and it is assumed that it will be rewritten and used as needed. The rewritten part will be mainly the scale of the drawing range and the cross-sectional force. Since the colors are not elaborate, I think that it should be rewritten according to the mood at that time.
Diagram that can be output
- Axial force diagram
- Shear force diagram
- Moment diagram
- Displacement mode diagram
What you are doing
- Create a section force diagram / displacement mode diagram using matplotlib with a Python program
- Imagine an attachment to a report, create a pdf with TeX, and convert it to png with ImageMagick
Output diagram example
Sectional force diagram creation program
As a point to devise, when writing the maximum and minimum values of the sectional force and displacement in the graph, these numerical values are referred to as "legend" in order to give versatility to the specification of the writing position and format. Is it written as?
In the program, the common part in the drawing of matplotlib is described once, and each section force diagram is drawn in a for loop.
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 command
By TeX, four graphs are arranged into one A3 horizontal sheet.
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 command execution script
convert is an ImageMagick command that removes margins from pdf and converts it to a png image.
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
that's all
Recommended Posts