Planar skeleton analysis with Python (4) Handling of forced displacement

Overview

In this program, boundary conditions are introduced without changing the size of the stiffness matrix, and simultaneous equations are solved. While using this method, consider a solution when a non-zero forced displacement is given. This is the replacement of unknowns and known numbers in simultaneous linear equations, and is not so difficult as a way of thinking.

As a simple example, consider the following simultaneous linear equations. $ \begin{align} & k_{11} x_1 + k_{12} x_2 + k_{13} x_3 = f_1 \\\ & k_{21} x_2 + k_{22} x_2 + k_{23} x_3 = f_2 \\\ & k_{31} x_3 + k_{32} x_2 + k_{33} x_3 = f_3 \\\ \end{align} $

Here, if the unknown number is $ x_1, x_3, f_2 $ and the known number is $ f_1, f_3, x_2 , it is necessary to rewrite as follows in order to perform the matrix operation. $ \begin{align} & k_{11} x_1 + 0 + k_{13} x_3 = f_1 - k_{12} x_2 \
& k_{21} x_2 - f_2 + k_{23} x_3 = 0 - k_{22} x_2 \
& k_{31} x_3 + 0 + k_{33} x_3 = f_3 - k_{32} x_2 \
\end{align} $$

The above relationship can be expanded and considered as follows.

Original stiffness equation

\begin{equation} \begin{bmatrix} k_{11} & \ldots & k_{1i} & \ldots & k_{1j} & \ldots & k_{1n} \\\ \vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \vdots \\\ k_{i1} & \ldots & k_{ii} & \ldots & k_{ij} & \ldots & k_{in} \\\ \vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \vdots \\\ k_{j1} & \ldots & k_{ji} & \ldots & k_{jj} & \ldots & k_{jn} \\\ \vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \vdots \\\ k_{n1} & \ldots & k_{ni} & \ldots & k_{nj} & \ldots & k_{nn} \end{bmatrix} \begin{Bmatrix} \delta_1 \\\ \vdots \\\ \delta_i \\\ \vdots \\\ \delta_j \\\ \vdots \\\ \delta_n \end{Bmatrix} =\begin{Bmatrix} f_1 \\\ \vdots \\\ f_i \\\ \vdots \\\ f_j \\\ \vdots \\\ f_n \end{Bmatrix} \end{equation}

Stiff equation after the introduction of boundary conditions

If $ \ delta_i $ and $ \ delta_j $ are known, let $ k_ {ii} $ and $ k_ {jj} $ in the stiffness matrix be $ 1 $, and the other elements of the $ i $ and $ j $ columns. Is set to zero. It also transfers the effects of the $ i $ and $ j $ columns in the stiffness matrix to the right side. $ \begin{equation} \begin{bmatrix} k_{11} & \ldots & 0 & \ldots & 0 & \ldots & k_{1n} \\\ \vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \vdots \\\ k_{i1} & \ldots & 1 & \ldots & 0 & \ldots & k_{in} \\\ \vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \vdots \\\ k_{j1} & \ldots & 0 & \ldots & 1 & \ldots & k_{jn} \\\ \vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \vdots \\\ k_{n1} & \ldots & 0 & \ldots & 0 & \ldots & k_{nn} \end{bmatrix} \begin{Bmatrix} \delta_1 \\\ \vdots \\\ -f_i \\\ \vdots \\\ -f_j \\\ \vdots \\\ \delta_n \end{Bmatrix} =\begin{Bmatrix} f_1 \\\ \vdots \\\ 0 \\\ \vdots \\\ 0 \\\ \vdots \\\ f_n \end{Bmatrix} -\delta_i \begin{Bmatrix} k_{1i} \\\ \vdots \\\ k_{ii} \\\ \vdots \\\ k_{ji} \\\ \vdots \\\ k_{ni} \end{Bmatrix} -\delta_j \begin{Bmatrix} k_{1j} \\\ \vdots \\\ k_{ij} \\\ \vdots \\\ k_{jj} \\\ \vdots \\\ k_{nj} \end{Bmatrix} \end{equation} $

Hotfix

`py_frame_r.py`


import numpy as np
import sys
import time

def print_inp(fnameW,npoin,nele,nsec,npfix,nlod,ae,x,fp,mpfix,rdis,node):
    fout=open(fnameW,'w')
    # print out of input data
    print('{0:>5s} {1:>5s} {2:>5s} {3:>5s} {4:>5s}'.format('npoin','nele','nsec','npfix','nlod'),file=fout)
    print('{0:5d} {1:5d} {2:5d} {3:5d} {4:5d}'.format(npoin,nele,nsec,npfix,nlod),file=fout)
    print('{0:>5s} {1:>15s} {2:>15s} {3:>15s}'.format('sec','A','I','E'),file=fout)
    for i in range(0,nsec):
        print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e}'.format(i+1,ae[0,i],ae[1,i],ae[2,i]),file=fout)
    print('{0:>5s} {1:>15s} {2:>15s} {3:>15s} {4:>15s} {5:>15s} {6:>5s} {7:>5s} {8:>5s}'
    .format('node','x','y','fx','fy','fr','kox','koy','kor'),file=fout)
    for i in range(0,npoin):
        lp=i+1
        print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e} {4:15.7e} {5:15.7e} {6:5d} {7:5d} {8:5d}'
        .format(lp,x[0,i],x[1,i],fp[3*lp-3],fp[3*lp-2],fp[3*lp-1],mpfix[0,i],mpfix[1,i],mpfix[2,i]),file=fout)
    print('{0:>5s} {1:>5s} {2:>5s} {3:>5s} {4:>15s} {5:>15s} {6:>15s}'
    .format('node','kox','koy','kor','rdis_x','rdis_y','rdis_r'),file=fout)
    for i in range(0,npoin):
        if 0<mpfix[0,i]+mpfix[1,i]+mpfix[2,i]:
            lp=i+1
            print('{0:5d} {1:5d} {2:5d} {3:5d} {4:15.7e} {5:15.7e} {6:15.7e}'
            .format(lp,mpfix[0,i],mpfix[1,i],mpfix[2,i],rdis[0,i],rdis[1,i],rdis[2,i]),file=fout)
    print('{0:>5s} {1:>5s} {2:>5s} {3:>5s}'.format('elem','i','j','sec'),file=fout)
    for ne in range(0,nele):
        print('{0:5d} {1:5d} {2:5d} {3:5d}'.format(ne+1,node[0,ne],node[1,ne],node[2,ne]),file=fout)
    fout.close()

def print_out(fnameW,npoin,nele,disg,fsec):
    fout=open(fnameW,'a')
    # displacement
    print('{0:>5s} {1:>15s} {2:>15s} {3:>15s}'.format('node','dis-x','dis-y','dis-r'),file=fout)
    for i in range(0,npoin):
        lp=i+1
        print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e}'.format(lp,disg[3*lp-3],disg[3*lp-2],disg[3*lp-1]),file=fout)
    # section force
    print('{0:>5s} {1:>15s} {2:>15s} {3:>15s} {4:>15s} {5:>15s} {6:>15s}'
    .format('elem','N_i','S_i','M_i','N_j','S_j','M_j'),file=fout)
    for ne in range(0,nele):
        print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e} {4:15.7e} {5:15.7e} {6:15.7e}'
        .format(ne+1,fsec[0,ne],fsec[1,ne],fsec[2,ne],fsec[3,ne],fsec[4,ne],fsec[5,ne]),file=fout)
    fout.close()


def STIFF(i,node,x,ae):
    ek=np.zeros([6,6],dtype=np.float64) # local stiffness matrix
    tt=np.zeros([6,6],dtype=np.float64) # transformation matrix
    mpf=node[0,i]
    mps=node[1,i]
    x1=x[0,mpf-1]
    y1=x[1,mpf-1]
    x2=x[0,mps-1]
    y2=x[1,mps-1]
    xx=x2-x1
    yy=y2-y1
    el=np.sqrt(xx**2+yy**2)
    mb=node[2,i]
    aa=ae[0,mb-1]
    am=ae[1,mb-1]
    ee=ae[2,mb-1]
    EA=ee*aa
    EI=ee*am
    ek[0,0]= EA/el; ek[0,1]=         0.0; ek[0,2]=        0.0; ek[0,3]=-EA/el; ek[0,4]=         0.0; ek[0,5]=        0.0
    ek[1,0]=   0.0; ek[1,1]= 12*EI/el**3; ek[1,2]= 6*EI/el**2; ek[1,3]=   0.0; ek[1,4]=-12*EI/el**3; ek[1,5]= 6*EI/el**2
    ek[2,0]=   0.0; ek[2,1]=  6*EI/el**2; ek[2,2]= 4*EI/el   ; ek[2,3]=   0.0; ek[2,4]= -6*EI/el**2; ek[2,5]= 2*EI/el
    ek[3,0]=-EA/el; ek[3,1]=         0.0; ek[3,2]=        0.0; ek[3,3]= EA/el; ek[3,4]=         0.0; ek[3,5]=        0.0
    ek[4,0]=   0.0; ek[4,1]=-12*EI/el**3; ek[4,2]=-6*EI/el**2; ek[4,3]=   0.0; ek[4,4]= 12*EI/el**3; ek[4,5]=-6*EI/el**2
    ek[5,0]=   0.0; ek[5,1]=  6*EI/el**2; ek[5,2]= 2*EI/el   ; ek[5,3]=   0.0; ek[5,4]= -6*EI/el**2; ek[5,5]= 4*EI/el
    s=yy/el
    c=xx/el
    tt[0,0]=  c; tt[0,1]=  s; tt[0,2]=0.0; tt[0,3]=0.0; tt[0,4]=0.0; tt[0,5]=0.0
    tt[1,0]= -s; tt[1,1]=  c; tt[1,2]=0.0; tt[1,3]=0.0; tt[1,4]=0.0; tt[1,5]=0.0
    tt[2,0]=0.0; tt[2,1]=0.0; tt[2,2]=1.0; tt[2,3]=0.0; tt[2,4]=0.0; tt[2,5]=0.0
    tt[3,0]=0.0; tt[3,1]=0.0; tt[3,2]=0.0; tt[3,3]=  c; tt[3,4]=  s; tt[3,5]=0.0
    tt[4,0]=0.0; tt[4,1]=0.0; tt[4,2]=0.0; tt[4,3]= -s; tt[4,4]=  c; tt[4,5]=0.0
    tt[5,0]=0.0; tt[5,1]=0.0; tt[5,2]=0.0; tt[5,3]=0.0; tt[5,4]=0.0; tt[5,5]=1.0
    return ek,tt,mpf,mps

# Main routine
start=time.time()
args = sys.argv
fnameR=args[1] # input data file
fnameW=args[2] # output data file

f=open(fnameR,'r')
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
n=3*npoin

x    =np.zeros([2,npoin],dtype=np.float64) # Coordinates of nodes
ae   =np.zeros([3,nsec],dtype=np.float64)  # Section characteristics
node =np.zeros([3,nele],dtype=np.int)      # Node-element relationship
fp   =np.zeros(3*npoin,dtype=np.float64)   # External force vector
mpfix=np.zeros([3,npoin],dtype=np.int)     # Ristrict conditions
rdis =np.zeros([3,npoin],dtype=np.float64) # Ristricted displacement
ir   =np.zeros(6,dtype=np.int)             # Work vector for matrix assembly
gk   =np.zeros([n,n],dtype=np.float64)     # Global stiffness matrix
fsec =np.zeros([6,nele],dtype=np.float64)  # Section force vector
work =np.zeros(6,dtype=np.float64)         # work vector for section force calculation

# section characteristics
for i in range(0,nsec):
    text=f.readline()
    text=text.strip()
    text=text.split()
    ae[0,i]=float(text[0]) #section area
    ae[1,i]=float(text[1]) #moment of inertia
    ae[2,i]=float(text[2]) #elastic modulus
# element-node
for i in range(0,nele):
    text=f.readline()
    text=text.strip()
    text=text.split()
    node[0,i]=int(text[0]) #node_1
    node[1,i]=int(text[1]) #node_2
    node[2,i]=int(text[2]) #section characteristic number
# node coordinates
for i in range(0,npoin):
    text=f.readline()
    text=text.strip()
    text=text.split()
    x[0,i]=float(text[0]) # x-coordinate
    x[1,i]=float(text[1]) # y-coordinate
# boundary conditions (0:free, 1:restricted)
for i in range(0,npfix):
    text=f.readline()
    text=text.strip()
    text=text.split()
    lp=int(text[0])             #fixed node
    mpfix[0,lp-1]=int(text[1])  #fixed in x-direction
    mpfix[1,lp-1]=int(text[2])  #fixed in y-direction
    mpfix[2,lp-1]=int(text[3])  #fixed in rotation
    rdis[0,lp-1]=float(text[4]) #fixed displacement in x-direction
    rdis[1,lp-1]=float(text[5]) #fixed displacement in y-direction
    rdis[2,lp-1]=float(text[6]) #fixed displacement in z-direction
# load
if 1<=nlod:
    for i in range(0,nlod):
        text=f.readline()
        text=text.strip()
        text=text.split()
        lp=int(text[0])           #loaded node
        fp[3*lp-3]=float(text[1]) #load in x-direction
        fp[3*lp-2]=float(text[2]) #load in y-direction
        fp[3*lp-1]=float(text[3]) #load in rotation
f.close()

print_inp(fnameW,npoin,nele,nsec,npfix,nlod,ae,x,fp,mpfix,rdis,node)

# assembly of stiffness matrix
for ne in range(0,nele):
    ek,tt,mpf,mps=STIFF(ne,node,x,ae)
    ck=np.dot(np.dot(tt.T,ek),tt)
    ir[5]=3*mps-1
    ir[4]=ir[5]-1
    ir[3]=ir[4]-1
    ir[2]=3*mpf-1
    ir[1]=ir[2]-1
    ir[0]=ir[1]-1
    for i in range(0,6):
        for j in range(0,6):
            it=ir[i]
            jt=ir[j]
            gk[it,jt]=gk[it,jt]+ck[i,j]

# treatment of boundary conditions
for i in range(0,npoin):
    for j in range(0,3):
        if mpfix[j,i]==1:
            iz=i*3+j
            fp[iz]=0.0
            for k in range(0,n):
                fp[k]=fp[k]-rdis[j,i]*gk[k,iz]
                gk[k,iz]=0.0
            gk[iz,iz]=1.0

# solution of simultaneous linear equations
disg = np.linalg.solve(gk, fp)

# recovery of restricted displacements
for i in range(0,npoin):
    for j in range(0,3):
        if mpfix[j,i]==1:
            iz=i*3+j
            for k in range(0,n):
                disg[iz]=rdis[j,i]

# calculation of section force
for ne in range(0,nele):
    ek,tt,mpf,mps=STIFF(ne,node,x,ae)
    work[0]=disg[3*mpf-3]; work[1]=disg[3*mpf-2]; work[2]=disg[3*mpf-1]
    work[3]=disg[3*mps-3]; work[4]=disg[3*mps-2]; work[5]=disg[3*mps-1]
    fsec[:,ne]=np.dot(ek,np.dot(tt,work))

print_out(fnameW,npoin,nele,disg,fsec)

dtime=time.time()-start
print('n={0}  time={1:.3f}'.format(n,dtime)+' sec')
fout=open(fnameW,'a')
print('n={0}  time={1:.3f}'.format(n,dtime)+' sec',file=fout)
fout.close()

Data creation and execution script

cat << EOT > inp_r0.txt
6 5 2 3 1
1.0 1.0 1.0
1.0 1.0 1.0
1 2 1
3 4 1
5 6 1
2 4 2
4 6 2
0 0
0 3.5
6 0
6 3.5
12 0
12 3.5
1 1 1 1 0 0 0
3 1 1 1 0 0 0
5 1 1 1 0 0 0
2 4 0 0
EOT

cat << EOT > inp_r1.txt
6 5 2 4 0
1.0 1.0 1.0
1.0 1.0 1.0
1 2 1
3 4 1
5 6 1
2 4 2
4 6 2
0 0
0 3.5
6 0
6 3.5
12 0
12 3.5
1 1 1 1 0 0 0
3 1 1 1 0 0 0
5 1 1 1 0 0 0
2 1 1 1 16.079284 2.3039125 -4.5858390
EOT

python3 py_frame_r.py inp_r0.txt out_r0.txt
python3 py_frame_r.py inp_r1.txt out_r1.txt

Execution example

Case-1 Specify load

`out_r0.txt`


npoin  nele  nsec npfix  nlod
    6     5     2     3     1
  sec               A               I               E
    1   1.0000000e+00   1.0000000e+00   1.0000000e+00
    2   1.0000000e+00   1.0000000e+00   1.0000000e+00
 node               x               y              fx              fy              fr   kox   koy   kor
    1   0.0000000e+00   0.0000000e+00   0.0000000e+00   0.0000000e+00   0.0000000e+00     1     1     1
    2   0.0000000e+00   3.5000000e+00   4.0000000e+00   0.0000000e+00   0.0000000e+00     0     0     0
    3   6.0000000e+00   0.0000000e+00   0.0000000e+00   0.0000000e+00   0.0000000e+00     1     1     1
    4   6.0000000e+00   3.5000000e+00   0.0000000e+00   0.0000000e+00   0.0000000e+00     0     0     0
    5   1.2000000e+01   0.0000000e+00   0.0000000e+00   0.0000000e+00   0.0000000e+00     1     1     1
    6   1.2000000e+01   3.5000000e+00   0.0000000e+00   0.0000000e+00   0.0000000e+00     0     0     0
 node   kox   koy   kor          rdis_x          rdis_y          rdis_r
    1     1     1     1   0.0000000e+00   0.0000000e+00   0.0000000e+00
    3     1     1     1   0.0000000e+00   0.0000000e+00   0.0000000e+00
    5     1     1     1   0.0000000e+00   0.0000000e+00   0.0000000e+00
 elem     i     j   sec
    1     1     2     1
    2     3     4     1
    3     5     6     1
    4     2     4     2
    5     4     6     2
 node           dis-x           dis-y           dis-r
    1   0.0000000e+00   0.0000000e+00   0.0000000e+00
    2   1.6079284e+01   2.3039125e+00  -4.5858390e+00
    3   0.0000000e+00   0.0000000e+00   0.0000000e+00
    4   5.6044784e+00  -1.4855500e+00  -6.2687943e-01
    5   0.0000000e+00   0.0000000e+00   0.0000000e+00
    6   2.6990174e+00  -8.1836247e-01  -5.5363182e-01
 elem             N_i             S_i             M_i             N_j             S_j             M_j
    1  -6.5826071e-01   2.2541991e+00   5.2550881e+00   6.5826071e-01  -2.2541991e+00   2.6346087e+00
    2   4.2444286e-01   1.2615574e+00   2.3868338e+00  -4.2444286e-01  -1.2615574e+00   2.0286170e+00
    3   2.3381785e-01   4.8424351e-01   1.0056067e+00  -2.3381785e-01  -4.8424351e-01   6.8924562e-01
    4   1.7458009e+00  -6.5826071e-01  -2.6346087e+00  -1.7458009e+00   6.5826071e-01  -1.3149555e+00
    5   4.8424351e-01  -2.3381785e-01  -7.1366148e-01  -4.8424351e-01   2.3381785e-01  -6.8924562e-01
n=18  time=0.005 sec

Case-2 Forced displacement designation

Input the displacement of the load specified point in Case-1 as the forced displacement

`out_r1.txt`


npoin  nele  nsec npfix  nlod
    6     5     2     4     0
  sec               A               I               E
    1   1.0000000e+00   1.0000000e+00   1.0000000e+00
    2   1.0000000e+00   1.0000000e+00   1.0000000e+00
 node               x               y              fx              fy              fr   kox   koy   kor
    1   0.0000000e+00   0.0000000e+00   0.0000000e+00   0.0000000e+00   0.0000000e+00     1     1     1
    2   0.0000000e+00   3.5000000e+00   0.0000000e+00   0.0000000e+00   0.0000000e+00     1     1     1
    3   6.0000000e+00   0.0000000e+00   0.0000000e+00   0.0000000e+00   0.0000000e+00     1     1     1
    4   6.0000000e+00   3.5000000e+00   0.0000000e+00   0.0000000e+00   0.0000000e+00     0     0     0
    5   1.2000000e+01   0.0000000e+00   0.0000000e+00   0.0000000e+00   0.0000000e+00     1     1     1
    6   1.2000000e+01   3.5000000e+00   0.0000000e+00   0.0000000e+00   0.0000000e+00     0     0     0
 node   kox   koy   kor          rdis_x          rdis_y          rdis_r
    1     1     1     1   0.0000000e+00   0.0000000e+00   0.0000000e+00
    2     1     1     1   1.6079284e+01   2.3039125e+00  -4.5858390e+00
    3     1     1     1   0.0000000e+00   0.0000000e+00   0.0000000e+00
    5     1     1     1   0.0000000e+00   0.0000000e+00   0.0000000e+00
 elem     i     j   sec
    1     1     2     1
    2     3     4     1
    3     5     6     1
    4     2     4     2
    5     4     6     2
 node           dis-x           dis-y           dis-r
    1   0.0000000e+00   0.0000000e+00   0.0000000e+00
    2   1.6079284e+01   2.3039125e+00  -4.5858390e+00
    3   0.0000000e+00   0.0000000e+00   0.0000000e+00
    4   5.6044785e+00  -1.4855500e+00  -6.2687945e-01
    5   0.0000000e+00   0.0000000e+00   0.0000000e+00
    6   2.6990174e+00  -8.1836249e-01  -5.5363183e-01
 elem             N_i             S_i             M_i             N_j             S_j             M_j
    1  -6.5826071e-01   2.2541992e+00   5.2550882e+00   6.5826071e-01  -2.2541992e+00   2.6346088e+00
    2   4.2444286e-01   1.2615574e+00   2.3868339e+00  -4.2444286e-01  -1.2615574e+00   2.0286170e+00
    3   2.3381785e-01   4.8424351e-01   1.0056067e+00  -2.3381785e-01  -4.8424351e-01   6.8924562e-01
    4   1.7458009e+00  -6.5826071e-01  -2.6346087e+00  -1.7458009e+00   6.5826071e-01  -1.3149555e+00
    5   4.8424351e-01  -2.3381785e-01  -7.1366150e-01  -4.8424351e-01   2.3381785e-01  -6.8924562e-01
n=18  time=0.007 sec

Sectional force diagram creation program

Since the output items of the analysis results were changed, the section force diagram creation program was also changed.

`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,npfix):
    text=f.readline()

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()

that's all