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

What you are doing

Output diagram example

tex_fig.png

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

Plane skeleton analysis with Python (3) Creation of section force diagram
Planar skeleton analysis with Python
Planar skeleton analysis with Python (4) Handling of forced displacement
3D skeleton structure analysis with Python
Example of 3D skeleton analysis by Python
Static analysis of Python code with GitLab CI
Two-dimensional elastic skeleton geometric nonlinear analysis with Python
Data analysis with python 2
Voice analysis with python
Voice analysis with python
Data analysis with Python
[OpenCV / Python] I tried image analysis of cells with OpenCV
Calculate the regression coefficient of simple regression analysis with python
Challenge principal component analysis of text data with Python
3. Natural language processing with Python 5-1. Concept of sentiment analysis [AFINN-111]
[Python] Morphological analysis with MeCab
[Co-occurrence analysis] Easy co-occurrence analysis with Python! [Python]
Perform isocurrent analysis of open channels with Python and matplotlib
Sentiment analysis with Python (word2vec)
Static analysis of Python programs
Japanese morphological analysis with Python
Basic summary of data manipulation with Python Pandas-First half: Data creation & manipulation
Muscle jerk analysis with Python
Practical exercise of data analysis with Python ~ 2016 New Coder Survey Edition ~
Image processing with Python 100 knocks # 4 Binarization of Otsu (discriminant analysis method)
From the introduction of JUMAN ++ to morphological analysis of Japanese with Python
Impedance analysis (EIS) with python [impedance.py]
Text mining with Python ① Morphological analysis
Getting Started with Python Basics of Python
Life game with Python! (Conway's Game of Life)
10 functions of "language with battery" python
Planar skeleton analysis in Python (2) Hotfix
Implementation of Dijkstra's algorithm with python
Data analysis starting with python (data visualization 1)
Logistic regression analysis Self-made with python
Coexistence of Python2 and 3 with CircleCI (1.0)
Data analysis starting with python (data visualization 2)
Basic study of OpenCV with Python
You can do it with Python! Structural analysis of two-dimensional colloidal crystals
Basics of binarized image processing with Python
[In-Database Python Analysis Tutorial with SQL Server 2017]
Marketing analysis with Python ① Customer analysis (decyl analysis, RFM analysis)
Two-dimensional saturated-unsaturated osmotic flow analysis with Python
Execute Python script with cron of TS-220
Machine learning with python (2) Simple regression analysis
Check the existence of the file with python
Algorithm learned with Python 8th: Evaluation of algorithm
Make a relation diagram of Python module
2D FEM stress analysis program with Python
Clogged with python update of GCP console ①
Easy introduction of speech recognition with Python
Sentiment analysis of tweets with deep learning
Overlay background diagram, contour diagram, vector diagram with python
Tweet analysis with Python, Mecab and CaboCha
UnicodeEncodeError struggle with standard output of python3
1. Statistics learned with Python 1-3. Calculation of various statistics (statistics)
Principal component analysis with Power BI + Python
Drawing with Matrix-Reinventor of Python Image Processing-
Recommendation of Altair! Data visualization with Python
Data analysis starting with python (data preprocessing-machine learning)
Analysis of X-ray microtomography image by Python