有限要素法(FEM) を python で作成する ～damyarouさんのいじくる編～

#シリーズ目次

- [有限要素法(FEM) を python で作成する ～vba→python 翻訳編～](http://qiita.com/sasaco/items/917429a497c6e9a08401)
- 有限要素法(FEM) を python で作成する ～damyarouさんのいじくる編～
- 有限要素法(FEM) を python で作成する ～行列の解法編～（予定）


---

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

[前回](http://qiita.com/sasaco/items/917429a497c6e9a08401)の記事を投稿したところ
有限要素法(FEM) を python で作成している方が既にいらっしゃいました
***[damyarou](http://qiita.com/damyarou/items/2b83b9b20181c1b8e35d)***さん のソースを ***いじくって*** みたので、記事にします。

# サンプルソース と サンプルデータの入手
damyarouさんの [ホームページ](http://civilyarou.web.fc2.com/py_FEM/sub_j_py_fem.html#3Dfrm)からサンプルソースとサンプルデータを入手します。

- サンプルソース:[py_fem_3dfrm.py]
(http://civilyarou.web.fc2.com/py_FEM/data_3dfrm/py_fem_3dfrm.py) 立体骨組構造解析プログラム

- サンプルデータ:[inp_3dfrm_grid.txt](http://civilyarou.web.fc2.com/py_FEM/data_3dfrm/inp_3dfrm_grid.txt) FEM解析入力データ(1)

![fig_gmt_0she[1].png](https://qiita-image-store.s3.amazonaws.com/0/142847/1b28d490-082f-4e32-0b7e-86980f1813bf.png)

使い方などは damyarouさんの [ホームページ](http://civilyarou.web.fc2.com/py_FEM/sub_j_py_fem.html#3Dfrm)を参照してください
サンプルデータの結果はこうなります

```text

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

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

# ちょっと手を加えてみる
ただ、コピペ して計算しただけだと面白くないので
下記のクラスに分けてみました。

```py: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 | アウトプットデータを書き込む |


## インプットデータを読み込む
```py:FramePython.py

fnameR= 'inp_grid.txt'
fnameW= 'out_grid.txt'

inp = inpData(fnameR)
inp.PRINP_3DFRM()

```

```py: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

```


## マトリックス作成

```py: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
```

### 要素剛性マトリックス

```py: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

```

### 全体剛性マトリックス

```py: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]

```

### 荷重剛性マトリックス
```py: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 にするなど

```py: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) で解いている．

```py:FramePython.py
# solution of simultaneous linear equations
result = np.linalg.solve(gmat.gk, fmat.fp)
```

## アウトプットデータを書き込む

```py: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()

```

```py: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`](https://github.com/sasaco/FramePython.git)


環境が整ったら、ソースコードのディレクトリに移動して、FramePython.pyを叩くと解析が始まります。

```text

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

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

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