damyarouさんのいじくる編の概要

前回の記事を投稿したところ
有限要素法(FEM) を python で作成している方が既にいらっしゃいました
***damyarou***さんのソースを いじくって みたので、記事にします。

サンプルソースとサンプルデータの入手

damyarouさんのホームページからサンプルソースとサンプルデータを入手します。

サンプルソース:[py_fem_3dfrm.py]
(http://civilyarou.web.fc2.com/py_FEM/data_3dfrm/py_fem_3dfrm.py) 立体骨組構造解析プログラム
サンプルデータ:inp_3dfrm_grid.txt FEM解析入力データ(1)

使い方などは 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

・・・・・・<以下略>・・・・・・・

正しく解析が行われました。