概要
このプログラムでは，剛性マトリックスの大きさを変えずに境界条件を導入し，連立方程式を解いている．
この方法を用いながら，ゼロでない強制変位を与えた場合の解法を考える．これは，連立一次方程式における，未知数と既知数の入れ替えであり，考え方としてはそれほど難しいものではない．
簡単な事例として，以下の連立一次方程式を考える．
$$
\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()
以 上