Flugzeugskelettanalyse mit Python (4) Umgang mit erzwungener Verschiebung

Überblick

In diesem Programm werden Randbedingungen eingeführt, ohne die Größe der Steifigkeitsmatrix zu ändern, und simultane Gleichungen werden gelöst. Betrachten Sie bei Verwendung dieser Methode eine Lösungsmethode, wenn eine erzwungene Verschiebung ungleich Null angegeben wird. Dies ist das Ersetzen von Unbekannten und bekannten Zahlen in simultanen linearen Gleichungen und nicht so schwierig wie eine Denkweise.

Betrachten Sie als einfaches Beispiel die folgenden simultanen linearen Gleichungen. $ \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} $

Wenn die unbekannte Zahl $ x_1, x_3, f_2 $ und die bekannte Zahl $ f_1, f_3, x_2 $ ist, muss hier wie folgt umgeschrieben werden, um die Matrixoperation auszuführen. $ \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} $

Die obige Beziehung kann wie folgt erweitert und betrachtet werden.

Ursprüngliche Steifheitsgleichung

Starrheitsgleichung nach Einführung von Randbedingungen

Wenn $ \ delta_i $ und $ \ delta_j $ bekannt sind, sei $ k_ {ii} $ und $ k_ {jj} $ in der Steifheitsmatrix $ 1 $ und die anderen Elemente der Spalten $ i $ und $ j $. Wird auf Null gesetzt. Außerdem werden die Effekte der Spalten $ i $ und $ j $ in der Steifheitsmatrix auf die rechte Seite übertragen. $ \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()

Skript zur Datenerstellung und -ausführung

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

Ausführungsbeispiel

Fall 1 Geben Sie die Last an

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

Fall 2 Bezeichnung für erzwungene Verschiebung

Geben Sie die Verschiebung des in Fall 1 angegebenen Lastpunkts als erzwungene Verschiebung ein

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

Programm zur Erstellung eines Schnittkraftdiagramms

Da die Ausgabeelemente der Analyseergebnisse geändert wurden, wurde auch das Programm zur Erstellung der Querschnittskraftkarte geändert.

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

das ist alles