シリーズ目次
- 有限要素法(FEM) を python で作成する ～vba→python 翻訳編～
- 有限要素法(FEM) を python で作成する ～damyarouさんのいじくる編～
- 有限要素法(FEM) を python で作成する ～行列の解法編～（予定）
damyarouさんのいじくる編 の概要
前回の記事を投稿したところ
有限要素法(FEM) を python で作成している方が既にいらっしゃいました
damyarouさん のソースを いじくって みたので、記事にします。
サンプルソース と サンプルデータの入手
damyarouさんの ホームページからサンプルソースとサンプルデータを入手します。
サンプルソース:py_fem_3dfrm.py 立体骨組構造解析プログラム
サンプルデータ:inp_3dfrm_grid.txt FEM解析入力データ(1)
![fig_gmt_0she[1].png](https://qiita-user-contents.imgix.net/https%3A%2F%2Fqiita-image-store.s3.amazonaws.com%2F0%2F142847%2F1b28d490-082f-4e32-0b7e-86980f1813bf.png?ixlib=rb-4.0.0&auto=format&gif-q=60&q=75&s=181d1a7776b1665b25794603dc58eca5)
使い方などは damyarouさんの ホームページを参照してください
サンプルデータの結果はこうなります
npoin nele nsec npfix nlod
392 759 4 16 392
sec E po A J Iy Iz theta
sec alpha gamma gkX gkY gkZ
1 2.5500000e+06 2.0000000e-01 6.0000000e-01 0.0000000e+00 1.8000000e-02 5.0000000e-02 0.0000000e+00
・・・・・・<中略>・・・・・・・
758 337 321 2
759 321 303 2
node dis-x dis-y dis-z rot-x rot-y rot-z
1 0.0000000e+00 0.0000000e+00 0.0000000e+00 -1.6243780e-04 1.9410889e-04 0.0000000e+00
2 0.0000000e+00 0.0000000e+00 -1.4268342e-04 -1.1792525e-04 3.1402105e-04 0.0000000e+00
3 0.0000000e+00 0.0000000e+00 -2.2831328e-04 -5.0801073e-05 4.0898102e-04 0.0000000e+00
4 0.0000000e+00 0.0000000e+00 -2.4248770e-04 2.1716711e-05 4.5111583e-04 0.0000000e+00
5 0.0000000e+00 0.0000000e+00 -1.8982394e-04 7.9402252e-05 4.3238604e-04 0.0000000e+00
6 0.0000000e+00 0.0000000e+00 -9.5023442e-05 1.0306705e-04 3.6727765e-04 0.0000000e+00
・・・・・・<中略>・・・・・・・
elem nodei N_i Sy_i Sz_i Mx_i My_i Mz_i
elem nodej N_j Sy_j Sz_j Mx_j My_j Mz_j
1 1 0.0000000e+00 0.0000000e+00 9.9678698e+01 -7.9374360e+01 9.7946784e+01 0.0000000e+00
1 2 0.0000000e+00 0.0000000e+00 -9.9678698e+01 7.9374360e+01 -1.9762548e+02 0.0000000e+00
2 2 0.0000000e+00 0.0000000e+00 5.0466963e+01 -6.2857564e+01 1.9762548e+02 0.0000000e+00
2 3 0.0000000e+00 0.0000000e+00 -5.0466963e+01 6.2857564e+01 -2.4809244e+02 0.0000000e+00
3 3 0.0000000e+00 0.0000000e+00 -1.4652296e+01 -2.7890607e+01 2.4809244e+02 0.0000000e+00
3 4 0.0000000e+00 0.0000000e+00 1.4652296e+01 2.7890607e+01 -2.3344015e+02 0.0000000e+00
・・・・・・<以下略>・・・・・・・
ちょっと手を加えてみる
ただ、コピペ して計算しただけだと面白くないので
下記のクラスに分けてみました。
FramePython.py
from inpData import inpData
from kMatrix import kMatrix
from gMatrix import gMatrix
from fMatrix import fMatrix
from outData import outData
| クラス | 用途 |
|---|---|
| inpData | インプットデータを読み込む |
| kMatrix | 要素毎の剛性マトリックス |
| gMatrix | 全体系の剛性マトリックス |
| fMatrix | 荷重剛性マトリックス |
| outData | アウトプットデータを書き込む |
インプットデータを読み込む
FramePython.py
fnameR= 'inp_grid.txt'
fnameW= 'out_grid.txt'
inp = inpData(fnameR)
inp.PRINP_3DFRM()
inpData.py
def __init__(self, fnameR):
f=open(fnameR,'r')
self.setdata(f)
f.close()
def setdata(self, f):
text=f.readline()
text=text.strip()
text=text.split()
self.npoin=int(text[0]) # Number of nodes
self.nele =int(text[1]) # Number of elements
self.nsec =int(text[2]) # Number of sections
self.npfix=int(text[3]) # Number of restricted nodes
self.nlod =int(text[4]) # Number of loaded nodes
self.nod=2 # Number of nodes per element
self.nfree=6 # Degree of freedom per node
self.n12=self.nod*self.nfree
self.n = self.nfree*self.npoin
self.x =np.zeros([3, self.npoin], dtype=np.float64) # Coordinates of nodes
self.deltaT=np.zeros(self.npoin, dtype=np.float64) # Temperature change of nodes
self.ae =np.zeros([12, self.nsec], dtype=np.float64) # Section characteristics
self.node =np.zeros([self.nod+1, self.nele], dtype=np.int) # Node-element relationship
self.fp =np.zeros(self.nfree*self.npoin,dtype=np.float64) # External force vector
self.mpfix =np.zeros([self.nfree, self.npoin],dtype=np.int) # Ristrict conditions
self.rdis =np.zeros([self.nfree, self.npoin],dtype=np.float64) # Ristricted displacement
# section characteristics
for i in range(0, self.nsec):
text=f.readline()
text=text.strip()
text=text.split()
self.ae[0,i] =float(text[0]) # E : elastic modulus
self.ae[1,i] =float(text[1]) # po : Poisson's ratio
self.ae[2,i] =float(text[2]) # aa : section area
self.ae[3,i] =float(text[3]) # aix : tortional constant
self.ae[4,i] =float(text[4]) # aiy : moment of inertia aroutd y-axis
self.ae[5,i] =float(text[5]) # aiz : moment of inertia around z-axis
self.ae[6,i] =float(text[6]) # theta : chord angle
self.ae[7,i] =float(text[7]) # alpha : thermal expansion coefficient
self.ae[8,i] =float(text[8]) # gamma : unit weight of material
self.ae[9,i] =float(text[9]) # gkX : acceleration in X-direction
self.ae[10,i]=float(text[10]) # gkY : acceleration in Y-direction
self.ae[11,i]=float(text[11]) # gkZ : acceleration in Z-direction
# element-node
for i in range(0, self.nele):
text=f.readline()
text=text.strip()
text=text.split()
self.node[0,i]=int(text[0]) #node_1
self.node[1,i]=int(text[1]) #node_2
self.node[2,i]=int(text[2]) #section characteristic number
# node coordinates
for i in range(0, self.npoin):
text=f.readline()
text=text.strip()
text=text.split()
self.x[0,i]=float(text[0]) # x-coordinate
self.x[1,i]=float(text[1]) # y-coordinate
self.x[2,i]=float(text[2]) # z-coordinate
self.deltaT[i]=float(text[3]) # Temperature change
# boundary conditions (0:free, 1:restricted)
for i in range(0, self.npfix):
text=f.readline()
text=text.strip()
text=text.split()
lp=int(text[0]) #fixed node
self.mpfix[0,lp-1]=int(text[1]) #fixed in x-direction
self.mpfix[1,lp-1]=int(text[2]) #fixed in y-direction
self.mpfix[2,lp-1]=int(text[3]) #fixed in z-direction
self.mpfix[3,lp-1]=int(text[4]) #fixed in rotation around x-axis
self.mpfix[4,lp-1]=int(text[5]) #fixed in rotation around y-axis
self.mpfix[5,lp-1]=int(text[6]) #fixed in rotation around z-axis
self.rdis[0,lp-1]=float(text[7]) #fixed displacement in x-direction
self.rdis[1,lp-1]=float(text[8]) #fixed displacement in y-direction
self.rdis[2,lp-1]=float(text[9]) #fixed displacement in z-direction
self.rdis[3,lp-1]=float(text[10]) #fixed rotation around x-axis
self.rdis[4,lp-1]=float(text[11]) #fixed rotation around y-axis
self.rdis[5,lp-1]=float(text[12]) #fixed rotation around z-axis
# load
if 0<self.nlod:
for i in range(0, self.nlod):
text=f.readline()
text=text.strip()
text=text.split()
lp=int(text[0]) #loaded node
self.fp[6*lp-6]=float(text[1]) #load in x-direction
self.fp[6*lp-5]=float(text[2]) #load in y-direction
self.fp[6*lp-4]=float(text[3]) #load in z-direction
self.fp[6*lp-3]=float(text[4]) #moment around x-axis
self.fp[6*lp-2]=float(text[5]) #moment around y-axis
self.fp[6*lp-1]=float(text[6]) #moment around z-axis
マトリックス作成
FramePython.py
gmat = gMatrix(inp.n)
fmat = fMatrix(inp)
for ne in range(0,inp.nele):
k= kMatrix(ne, inp)
gmat.set_kMatrix(k) # assembly of stiffness matrix
fmat.set_fMatrix(k) # assembly of load vectors
要素剛性マトリックス
kMatrix.py
def __init__(self, ne, input):
self.ne = ne
self.inp = input
self.iNo = self.inp.node[0,self.ne]
self.jNo = self.inp.node[1,self.ne]
self.eNo = self.inp.node[2,self.ne]
i = self.iNo-1
j = self.jNo-1
self.x1= self.inp.x[0,i]
self.y1= self.inp.x[1,i]
self.z1= self.inp.x[2,i]
self.x2= self.inp.x[0,j]
self.y2= self.inp.x[1,j]
self.z2= self.inp.x[2,j]
xx = self.x2-self.x1
yy = self.y2-self.y1
zz = self.z2-self.z1
self.el = np.sqrt(xx**2+yy**2+zz**2)
m = self.eNo-1
self.ee = self.inp.ae[0,m] # elastic modulus
self.po = self.inp.ae[1,m] # Poisson's ratio
self.aa = self.inp.ae[2,m] # section area
self.aix = self.inp.ae[3,m] # tortional constant
self.aiy = self.inp.ae[4,m] # moment of inertia around y-axis
self.aiz = self.inp.ae[5,m] # moment of inertia around z-axis
self.theta = self.inp.ae[6,m] # theta : chord angle
self.alpha = self.inp.ae[7,m] # alpha : thermal expansion coefficient
self.gamma = self.inp.ae[8,m] # gamma : unit weight of material
self.gkX = self.inp.ae[9,m] # gkX : acceleration in X-direction
self.gkY = self.inp.ae[10,m] # gkY : acceleration in Y-direction
self.gkZ = self.inp.ae[11,m] # gkZ : acceleration in Z-direction
self.GJ=self.ee/2/(1+self.po)*self.aix
self.EA=self.ee*self.aa
self.EIy=self.ee*self.aiy
self.EIz=self.ee*self.aiz
# Work vector for matrix assembly
ir = np.zeros(12, dtype=np.int)
ir[11]=6*j+5; ir[10]=ir[11]-1; ir[9]=ir[10]-1; ir[8]=ir[9]-1; ir[7]=ir[8]-1; ir[6]=ir[7]-1
ir[5] =6*i+5; ir[4] =ir[5]-1 ; ir[3]=ir[4]-1 ; ir[2]=ir[3]-1; ir[1]=ir[2]-1; ir[0]=ir[1]-1
self.ir = ir
self.ek = self.SM_3DFRM()
self.tt = self.TM_3DFRM()
def SM_3DFRM(self):
ek = np.zeros([12,12],dtype=np.float64) # local stiffness matrix
ek[ 0, 0]= self.EA/self.el
ek[ 0, 6]=-self.EA/self.el
ek[ 1, 1]= 12*self.EIz/self.el**3
ek[ 1, 5]= 6*self.EIz/self.el**2
ek[ 1, 7]=-12*self.EIz/self.el**3
ek[ 1,11]= 6*self.EIz/self.el**2
ek[ 2, 2]= 12*self.EIy/self.el**3
ek[ 2, 4]= -6*self.EIy/self.el**2
ek[ 2, 8]=-12*self.EIy/self.el**3
ek[ 2,10]= -6*self.EIy/self.el**2
ek[ 3, 3]= self.GJ/self.el
ek[ 3, 9]=-self.GJ/self.el
ek[ 4, 2]= -6*self.EIy/self.el**2
ek[ 4, 4]= 4*self.EIy/self.el
ek[ 4, 8]= 6*self.EIy/self.el**2
ek[ 4,10]= 2*self.EIy/self.el
ek[ 5, 1]= 6*self.EIz/self.el**2
ek[ 5, 5]= 4*self.EIz/self.el
ek[ 5, 7]= -6*self.EIz/self.el**2
ek[ 5,11]= 2*self.EIz/self.el
ek[ 6, 0]=-self.EA/self.el
ek[ 6, 6]= self.EA/self.el
ek[ 7, 1]=-12*self.EIz/self.el**3
ek[ 7, 5]= -6*self.EIz/self.el**2
ek[ 7, 7]= 12*self.EIz/self.el**3
ek[ 7,11]= -6*self.EIz/self.el**2
ek[ 8, 2]=-12*self.EIy/self.el**3
ek[ 8, 4]= 6*self.EIy/self.el**2
ek[ 8, 8]= 12*self.EIy/self.el**3
ek[ 8,10]= 6*self.EIy/self.el**2
ek[ 9, 3]=-self.GJ/self.el
ek[ 9, 9]= self.GJ/self.el
ek[10, 2]= -6*self.EIy/self.el**2
ek[10, 4]= 2*self.EIy/self.el
ek[10, 8]= 6*self.EIy/self.el**2
ek[10,10]= 4*self.EIy/self.el
ek[11, 1]= 6*self.EIz/self.el**2
ek[11, 5]= 2*self.EIz/self.el
ek[11, 7]= -6*self.EIz/self.el**2
ek[11,11]= 4*self.EIz/self.el
return ek
def TM_3DFRM(self):
tt=np.zeros([12,12],dtype=np.float64) # transformation matrix
t1=np.zeros([3,3],dtype=np.float64)
t2=np.zeros([3,3],dtype=np.float64)
theta=np.radians(self.theta)#(self.inp.ae[6,m]) # chord angle
t1[0,0]=1
t1[1,1]= np.cos(theta)
t1[1,2]= np.sin(theta)
t1[2,1]=-np.sin(theta)
t1[2,2]= np.cos(theta)
ll=(self.x2-self.x1)/self.el
mm=(self.y2-self.y1)/self.el
nn=(self.z2-self.z1)/self.el
if self.x2-self.x1==0.0 and self.y2-self.y1==0.0:
t2[0,2]=nn
t2[1,0]=nn
t2[2,1]=1.0
else:
qq=np.sqrt(ll**2+mm**2)
t2[0,0]=ll
t2[0,1]=mm
t2[0,2]=nn
t2[1,0]=-mm/qq
t2[1,1]= ll/qq
t2[2,0]=-ll*nn/qq
t2[2,1]=-mm*nn/qq
t2[2,2]=qq
t3=np.dot(t1,t2)
tt[0:3,0:3] =t3[0:3,0:3]
tt[3:6,3:6] =t3[0:3,0:3]
tt[6:9,6:9] =t3[0:3,0:3]
tt[9:12,9:12]=t3[0:3,0:3]
return tt
def ck(self):
return np.dot(np.dot(self.tt.T,self.ek),self.tt) # Stiffness matrix in global coordinate
全体剛性マトリックス
gMatrix.py
def __init__(self, n):
self.gk = np.zeros([n,n], dtype=np.float64) # Global stiffness matrix
self.kmat = {}
def set_kMatrix(self, k):
self.kmat[k.ne] = k
ck = k.ck()
for i in range(0,12):
it=k.ir[i]
for j in range(0,12):
jt=k.ir[j]
self.gk[it,jt]=self.gk[it,jt]+ck[i,j]
荷重剛性マトリックス
fMatrix.py
def __init__(self, input):
self.inp = input
self.fp = self.inp.fp.copy()
def set_fMatrix(self, k):
tfe_l=self.TFVEC_3DFRM(k)
tt = k.tt
tfe =np.dot(tt.T,tfe_l) # Thermal load vector in global coordinate
bfe =self.BFVEC_3DFRM( k) # Body force vector in global coordinate
for i in range(0,12):
it=k.ir[i]
self.fp[it]=self.fp[it]+tfe[i]+bfe[i]
def TFVEC_3DFRM(self, k):
# Thermal load vector in local coordinate system
tfe_l=np.zeros(12,dtype=np.float64)
i= k.iNo-1
j= k.jNo-1
E = k.ee # elastic modulus
AA = k.aa # section area
alpha = k.alpha # thermal expansion coefficient
tempe=0.5*(self.inp.deltaT[i]+self.inp.deltaT[j])
tfe_l[0]=-E*AA*alpha*tempe
tfe_l[6]= E*AA*alpha*tempe
return tfe_l
def BFVEC_3DFRM(self, k):
# Body force vector in global coordinate system
bfe_g =np.zeros(12,dtype=np.float64)
AA = k.aa # section area
gamma = k.gamma # unit weight of material
gkX = k.gkX # seismic coefficient in X-direction
gkY = k.gkY # seismic coefficient in Y-direction
gkZ = k.gkZ # seismic coefficient in Z-direction
el = k.el
bfe_g[0]=0.5*gamma*AA*el*gkX
bfe_g[1]=0.5*gamma*AA*el*gkY
bfe_g[2]=0.5*gamma*AA*el*gkZ
bfe_g[6]=0.5*gamma*AA*el*gkX
bfe_g[7]=0.5*gamma*AA*el*gkY
bfe_g[8]=0.5*gamma*AA*el*gkZ
return bfe_g
固定支点の変位量を 0 にするなど
FramePython.py
# treatment of boundary conditions
for i in range(0,inp.npoin):
for j in range(0,inp.nfree):
if inp.mpfix[j,i]==1:
iz=i*inp.nfree+j
fmat.fp[iz]=0.0
for k in range(0,inp.n):
fmat.fp[k]=fmat.fp[k]-inp.rdis[j,i]*gmat.gk[k,iz]
gmat.gk[k,iz]=0.0
gmat.gk[iz,iz]=1.0
ソルバー
連立一次方程式は，numpy.linalg.solve(A, b) で解いている．
FramePython.py
# solution of simultaneous linear equations
result = np.linalg.solve(gmat.gk, fmat.fp)
アウトプットデータを書き込む
FramePython.py
out = outData(inp, result, gmat)
out.PROUT_3DFRM(fnameW)
dtime=time.time()-start
print('n={0} time={1:.3f}'.format(inp.n,dtime)+' sec')
fout=open(fnameW,'a')
print('n={0} time={1:.3f}'.format(inp.n,dtime)+' sec',file=fout)
fout.close()
outData.py
def __init__(self, input, result, gmat):
self.inp = input
self.disg = result
self.fsec = np.zeros([self.inp.n12, self.inp.nele], dtype=np.float64) # Section force vector
# recovery of restricted displacements
for i in range(0,self.inp.npoin):
for j in range(0,self.inp.nfree):
if self.inp.mpfix[j,i]==1:
iz=i*self.inp.nfree+j
self.disg[iz]=self.inp.rdis[j,i]
# calculation of section force
work =np.zeros(self.inp.n12, dtype=np.float64) # work vector for section force calculation
for ne in range(0,self.inp.nele):
k = gmat.kmat[ne]
i = k.iNo-1
j = k.jNo-1
ek = k.ek # Stiffness matrix in local coordinate
tt = k.tt # Transformation matrix
E = k.ee # elastic modulus
AA = k.aa # section area
alpha = k.alpha # thermal expansion coefficient
tempe = 0.5*(self.inp.deltaT[i]+self.inp.deltaT[j])
work[0]=self.disg[6*i] ; work[1] =self.disg[6*i+1]; work[2]= self.disg[6*i+2]
work[3]=self.disg[6*i+3]; work[4] =self.disg[6*i+4]; work[5]= self.disg[6*i+5]
work[6]=self.disg[6*j] ; work[7] =self.disg[6*j+1]; work[8]= self.disg[6*j+2]
work[9]=self.disg[6*j+3]; work[10]=self.disg[6*j+4]; work[11]=self.disg[6*j+5]
self.fsec[:,ne]=np.dot(ek,np.dot(tt,work))
self.fsec[0,ne]=self.fsec[0,ne]+E*AA*alpha*tempe
self.fsec[6,ne]=self.fsec[6,ne]-E*AA*alpha*tempe
結果
手を加えてたソースはここにおいておきます。
$ git clone -b "#2" https://github.com/sasaco/FramePython.git
環境が整ったら、ソースコードのディレクトリに移動して、FramePython.pyを叩くと解析が始まります。
npoin nele nsec npfix nlod
392 759 4 16 392
sec E po A J Iy Iz theta
sec alpha gamma gkX gkY gkZ
1 2.5500000e+06 2.0000000e-01 6.0000000e-01 0.0000000e+00 1.8000000e-02 5.0000000e-02 0.0000000e+00
・・・・・・<中略>・・・・・・・
758 337 321 2
759 321 303 2
node dis-x dis-y dis-z rot-x rot-y rot-z
1 0.0000000e+00 0.0000000e+00 0.0000000e+00 -1.6243780e-04 1.9410889e-04 0.0000000e+00
2 0.0000000e+00 0.0000000e+00 -1.4268342e-04 -1.1792525e-04 3.1402105e-04 0.0000000e+00
3 0.0000000e+00 0.0000000e+00 -2.2831328e-04 -5.0801073e-05 4.0898102e-04 0.0000000e+00
4 0.0000000e+00 0.0000000e+00 -2.4248770e-04 2.1716711e-05 4.5111583e-04 0.0000000e+00
5 0.0000000e+00 0.0000000e+00 -1.8982394e-04 7.9402252e-05 4.3238604e-04 0.0000000e+00
6 0.0000000e+00 0.0000000e+00 -9.5023442e-05 1.0306705e-04 3.6727765e-04 0.0000000e+00
・・・・・・<中略>・・・・・・・
elem nodei N_i Sy_i Sz_i Mx_i My_i Mz_i
elem nodej N_j Sy_j Sz_j Mx_j My_j Mz_j
1 1 0.0000000e+00 0.0000000e+00 9.9678698e+01 -7.9374360e+01 9.7946784e+01 0.0000000e+00
1 2 0.0000000e+00 0.0000000e+00 -9.9678698e+01 7.9374360e+01 -1.9762548e+02 0.0000000e+00
2 2 0.0000000e+00 0.0000000e+00 5.0466963e+01 -6.2857564e+01 1.9762548e+02 0.0000000e+00
2 3 0.0000000e+00 0.0000000e+00 -5.0466963e+01 6.2857564e+01 -2.4809244e+02 0.0000000e+00
3 3 0.0000000e+00 0.0000000e+00 -1.4652296e+01 -2.7890607e+01 2.4809244e+02 0.0000000e+00
3 4 0.0000000e+00 0.0000000e+00 1.4652296e+01 2.7890607e+01 -2.3344015e+02 0.0000000e+00
・・・・・・<以下略>・・・・・・・
正しく解析が行われました。