Analyse du squelette de plan avec Python (4) Gestion du déplacement forcé

Aperçu

Dans ce programme, les conditions aux limites sont introduites sans changer la taille de la matrice de rigidité, et les équations simultanées sont résolues. Lors de l'utilisation de cette méthode, considérez une méthode de résolution lorsqu'un déplacement forcé non nul est donné. Il s'agit du remplacement d'inconnues et de nombres connus dans des équations linéaires simultanées, et ce n'est pas si difficile qu'une façon de penser.

À titre d'exemple simple, considérons les équations linéaires simultanées suivantes. $ \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} $

Ici, si le nombre inconnu est $ x_1, x_3, f_2 $ et le nombre connu est $ f_1, f_3, x_2 , il est nécessaire de réécrire comme suit pour effectuer l'opération de matrice. $ \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} $$

La relation ci-dessus peut être développée et considérée comme suit.

Équation de rigidité originale

Équation de rigidité après introduction des conditions aux limites

Si $ \ delta_i $ et $ \ delta_j $ sont connus, soit $ k_ {ii} $ et $ k_ {jj} $ dans la matrice de rigidité $ 1 $, et les autres éléments des colonnes $ i $ et $ j $. Est mis à zéro. Il transfère également les effets des colonnes $ i $ et $ j $ de la matrice de rigidité vers le côté droit. $ \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} $

Correctif

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

Script de création et d'exécution de données

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

Exemple d'exécution

Cas-1 Spécifier la charge

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

Cas-2 Désignation du déplacement forcé

Entrez le déplacement du point spécifié de la charge dans le cas 1 comme déplacement forcé

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

Programme de création de diagramme de force sectionnelle

Étant donné que les éléments de sortie des résultats d'analyse ont été modifiés, le programme de création de cartes de forces transversales a également été modifié.

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

c'est tout