概要

このプログラムでは，剛性マトリックスの大きさを変えずに境界条件を導入し，連立方程式を解いている．
この方法を用いながら，ゼロでない強制変位を与えた場合の解法を考える．これは，連立一次方程式における，未知数と既知数の入れ替えであり，考え方としてはそれほど難しいものではない．

簡単な事例として，以下の連立一次方程式を考える．
$$
\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}
$$

ここで，未知数が $x_1, x_3, f_2$，既知数が $f_1, f_3, x_2$だとすれば，行列演算を行うためには，以下のように書き換える必要がある．
$$
\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}
$$

上記関係を，拡張して考えると，以下の通りとなる．

元の剛性方程式

$$
\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}
$$

境界条件導入後の剛性方程式

$\delta_i$ および $\delta_j$ が既知とすれば，剛性マトリックスの $k_{ii}$ および $k_{jj}$ を $1$ とし，それ以外の $i$ 列および $j$ 列の要素をゼロとする処理を行う．また剛性マトリックスにおける $i$ 列および $j$ 列の効果を右辺に移項する．
$$
\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}
$$

修正プログラム

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

データ作成と実行スクリプト

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

実行事例

Case-1 荷重を指定

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 強制変位指定

Case-1における荷重指定点の変位を強制変位として入力

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

断面力図作成プログラム

解析結果の出力項目を変更したので，断面力図作成プログラムの変更も行った．

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

以　上