3次元有限要素法をPythonで実装する(四面体要素)

概要

有限要素法といえば3次元での解析が一般的ですが、Qiitaに3次元解析の実装方法についての記事があまり無かったの投稿しました。
今回はオーソドックスな四面体要素を実装してみます。
簡単のため、まずは微小変形の弾性解析で実装してみます。
また、この記事では下記のサンプルを解く処理を実装してみます。

解くべき方程式

有限要素法では以下の式を解くことが目的となります。

\boldsymbol{K}\boldsymbol{\tilde{U}} = \boldsymbol{F} \\
\begin{align}
&\boldsymbol{K} : 剛性マトリクス \\
&\boldsymbol{\tilde{U}} : 節点の変位ベクトル \\
&\boldsymbol{F} : 節点に働く荷重のベクトル
\end{align}

ここで、剛性マトリクスは下記の式で表すことができます。

\boldsymbol{K} = \sum_e \boldsymbol{K}_e \\
\boldsymbol{K}_e : 要素剛性マトリクス

要素剛性マトリクスは要素毎の剛性マトリクスで、それを全て足し合わせることで剛性マトリクスになります。この手順はアセンブリングと呼ばれます。
このため、要素剛性マトリクスから剛性マトリクスを作成し、荷重ベクトルにその逆行列をかけることで変位ベクトルを求めることができます。

四面体要素

応力とひずみの関係式

3次元の応力とひずみの関係は下記のように表されます。

\left\{
\begin{array}{ll}
\epsilon_x = \frac{1}{E}(\sigma_x - \nu (\sigma_y + \sigma_z)) \\
\epsilon_y = \frac{1}{E}(\sigma_y - \nu (\sigma_x + \sigma_z)) \\
\epsilon_z = \frac{1}{E}(\sigma_z - \nu (\sigma_x + \sigma_y))　\\
\gamma_{xy} = \frac{2(1+\nu)}{E}\tau_{xy} \\
\gamma_{yz} = \frac{2(1+\nu)}{E}\tau_{yz} \\
\gamma_{xz} = \frac{2(1+\nu)}{E}\tau_{xz}
\end{array}
\right.
\\
\begin{align}
&\sigma : 垂直応力 \\
&\tau : せん断応力 \\
&\epsilon : 垂直ひずみ \\
&\gamma : 工学せん断ひずみ \\
&E : ヤング率 \\
&\nu : ポアソン比
\end{align}

これを行列形式にして逆行列をとると下記のように表すことができます。

\boldsymbol{\sigma} = \boldsymbol{D} \boldsymbol{\epsilon} \\
\boldsymbol{\sigma} = 
\begin{bmatrix}
\sigma_x \\
\sigma_y \\
\sigma_z \\
\tau_{yz} \\
\tau_{xz} \\
\tau_{xy}
\end{bmatrix}
,
\boldsymbol{\epsilon} =
\begin{bmatrix}
\epsilon_x \\
\epsilon_y \\
\epsilon_z \\
\gamma_{yz} \\
\gamma_{xz} \\
\gamma_{xy} 
\end{bmatrix}
\\
\boldsymbol{D} =
\frac{E}{(1+\nu)(1-2\nu)}
\begin{bmatrix}
1-\nu && \nu && \nu && 0 && 0 && 0 \\
\nu && 1-\nu && \nu && 0 && 0 && 0  \\
\nu && \nu && 1-\nu && 0 && 0 && 0  \\
0 && 0 && 0 && \frac{1-2\nu}{2} && 0 && 0 \\
0 && 0 && 0 && 0 && \frac{1-2\nu}{2} && 0  \\
0 && 0 && 0 && 0 && 0 && \frac{1-2\nu}{2}
\end{bmatrix}

ひずみと変位の関係式

ひずみは変位を偏微分した下記の式で表すことができます。

\boldsymbol{\epsilon} = 
\begin{bmatrix}
\epsilon_x \\
\epsilon_y \\
\epsilon_z \\
\gamma_{yz} \\ 
\gamma_{xz} \\
\gamma_{xy} 
\end{bmatrix}
=
\begin{bmatrix}
\frac{\partial u}{\partial x} \\
\frac{\partial v}{\partial y} \\
\frac{\partial w}{\partial z} \\
\frac{\partial w}{\partial y} + \frac{\partial v}{\partial z} \\
\frac{\partial u}{\partial z} + \frac{\partial w}{\partial x} \\
\frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} 
\end{bmatrix}
\\
\begin{align}
&u : x方向の変位 \\
&v : y方向の変位 \\
&w : z方向の変位 \\
\end{align}

ここで、四面体要素の変位の分布は形状関数を使って表すことができるという特徴があります。

\left\{
\begin{array}{ll}
u = N_1 u_1 + N_2 u_2 + N_3 u_3 + N_4 u_4 \\
v = N_1 v_1 + N_2 v_2 + N_3 v_3 + N_4 v_4 \\
w = N_1 w_1 + N_2 w_2 + N_3 w_3 + N_4 w_4 
\end{array}
\right.
\\
\left\{
\begin{array}{ll}
N_1 = 1-a-b-c\\
N_2 = a \\
N_3 = b \\
N_4 = c
\end{array}
\right.
\\
\begin{align}
&u_i : i節点のx方向の変位 \\
&v_i : i節点のy方向の変位 \\
&w_i : i節点のz方向の変位 \\
&N_i : 形状関数 \\
\end{align}

形状関数は正規化座標系a,b,cで表され、上の式はx,y,z座標系から正規化座標系a,b,cに変換することを表しています。
図にすると下記のようになります。

上記の形状関数を使うことにより変位の偏微分は下記の式で表すことができます。

\frac{\partial u}{\partial x} = \frac{\partial N_1}{\partial x}u_1 + \frac{\partial N_2}{\partial x}u_2 + \frac{\partial N_3}{\partial x}u_3 +
\frac{\partial N_4}{\partial x}u_4 \\

u,wに関する偏微分も同様の形式になります。
形状関数を使うことでひずみは節点の変位とBマトリクスによって表すことができます。

\boldsymbol{\epsilon} = \boldsymbol{B} \tilde{\boldsymbol{u}} \\
\tilde{\boldsymbol{u}} = 
\begin{bmatrix}
u_1 \\
v_1 \\
w_1 \\
u_2 \\
v_2 \\
w_2 \\
u_3 \\
v_3 \\
w_3 \\
u_4 \\
v_4 \\
w_4
\end{bmatrix}
\\
\boldsymbol{B} =
\begin{bmatrix}
\frac{\partial N_1}{\partial x} && 0 && 0 && \frac{\partial N_2}{\partial x} && 0 && 0 && \frac{\partial N_3}{\partial x} && 0 && 0 && \frac{\partial N_4}{\partial x} && 0 && 0 \\
0 && \frac{\partial N_1}{\partial y} && 0 && 0 && \frac{\partial N_2}{\partial y} && 0 && 0 && \frac{\partial N_3}{\partial y} && 0 && 0 && \frac{\partial N_4}{\partial y} && 0 \\
0 && 0 && \frac{\partial N_1}{\partial z} && 0 && 0 && \frac{\partial N_2}{\partial z} && 0 && 0 && \frac{\partial N_3}{\partial z} && 0 && 0 && \frac{\partial N_4}{\partial z} \\
0 && \frac{\partial N_1}{\partial z} && \frac{\partial N_1}{\partial y} && 0 && \frac{\partial N_2}{\partial z} && \frac{\partial N_2}{\partial y} && 0 && \frac{\partial N_3}{\partial z} && \frac{\partial N_3}{\partial y} && 0 && \frac{\partial N_4}{\partial z} && \frac{\partial N_4}{\partial y} \\
\frac{\partial N_1}{\partial z} && 0 && \frac{\partial N_1}{\partial x} && \frac{\partial N_2}{\partial z} && 0 && \frac{\partial N_2}{\partial x} &&
\frac{\partial N_3}{\partial z} && 0 && \frac{\partial N_3}{\partial x} &&
\frac{\partial N_4}{\partial z} && 0 && \frac{\partial N_4}{\partial x} \\
\frac{\partial N_1}{\partial y} && \frac{\partial N_1}{\partial x} && 0 && \frac{\partial N_2}{\partial y} && \frac{\partial N_2}{\partial x} && 0 && \frac{\partial N_3}{\partial y} && \frac{\partial N_3}{\partial x} && 0 && \frac{\partial N_4}{\partial y} && \frac{\partial N_4}{\partial x} && 0 

\end{bmatrix}

ここで、形状関数を偏微分した値を求めることができれば、Bマトリクスを求めることができます。形状関数を偏微分した値はヤコビ行列から求めることが可能です。
ヤコビ行列は下記の式で定義できます。

\boldsymbol{J}=
\begin{bmatrix}
\frac{\partial x}{\partial a} && \frac{\partial y}{\partial a} && \frac{\partial z}{\partial a} \\
\frac{\partial x}{\partial b} && \frac{\partial y}{\partial b} && \frac{\partial z}{\partial b} \\
\frac{\partial x}{\partial c} && \frac{\partial y}{\partial c} && \frac{\partial z}{\partial c} 
\end{bmatrix}
\\
\begin{align}
&\frac{\partial x}{\partial a} = \sum_{i=1}^4 \frac{\partial N_i}{\partial a}x_i & \frac{\partial x}{\partial b} = \sum_{i=1}^4 \frac{\partial N_i}{\partial b}x_i & & \frac{\partial x}{\partial c} = \sum_{i=1}^4 \frac{\partial N_i}{\partial c}x_i \\
&\frac{\partial y}{\partial a} = \sum_{i=1}^4 \frac{\partial N_i}{\partial a}y_i &  \frac{\partial y}{\partial b} = \sum_{i=1}^4 \frac{\partial N_i}{\partial b}y_i & & \frac{\partial y}{\partial c} = \sum_{i=1}^4 \frac{\partial N_i}{\partial c}y_i\\
&\frac{\partial z}{\partial a} = \sum_{i=1}^4 \frac{\partial N_i}{\partial a}z_i & \frac{\partial z}{\partial b} = \sum_{i=1}^4 \frac{\partial N_i}{\partial b}z_i & & \frac{\partial z}{\partial c} = \sum_{i=1}^4 \frac{\partial N_i}{\partial c}z_i\\
\end{align}

ヤコビ行列を使うと形状関数の偏微分は下記の式で表すことができます。

\begin{bmatrix}
\frac{\partial N_i}{\partial a} \\
\frac{\partial N_i}{\partial b} \\
\frac{\partial N_i}{\partial c} 
\end{bmatrix}
=
\boldsymbol{J}
\begin{bmatrix}
\frac{\partial N_i}{\partial x} \\
\frac{\partial N_i}{\partial y} \\
\frac{\partial N_i}{\partial z} 
\end{bmatrix}
\\
\begin{bmatrix}
\frac{\partial N_i}{\partial x} \\
\frac{\partial N_i}{\partial y} \\
\frac{\partial N_i}{\partial z} 
\end{bmatrix}
=
\boldsymbol{J}^{-1}
\begin{bmatrix}
\frac{\partial N_i}{\partial a} \\
\frac{\partial N_i}{\partial b} \\
\frac{\partial N_i}{\partial c} 
\end{bmatrix}

要素剛性マトリクス

四面体要素に仮想仕事の原理を適用すると下記のように表すことができます。

\int_v \delta \boldsymbol{\epsilon}^T \boldsymbol{\sigma} dV = 
\int_v \delta \boldsymbol{u}^T \boldsymbol{b} dV + \int_s \delta \boldsymbol{u}^T \boldsymbol{t} dS 
\\
\delta \boldsymbol{\epsilon} : 仮想ひずみ \\
\delta \boldsymbol{u} : 仮想変位 \\
\sigma : 応力 \\
\boldsymbol{b} : 物体力　\\
\boldsymbol{t} : 表面力 \\

仮想仕事の原理の導出については下記を参照ください

・仮想仕事の原理
・有限要素法解析Part 1：仮想仕事の原理 - YouTube

ここで、要素内の変位とひずみは下記のように表すことができます。

\boldsymbol{u} = \boldsymbol{N} \tilde{\boldsymbol{u}} \\
\boldsymbol{N} = 
\begin{bmatrix}
N_1 && 0 && N_2 && 0 && N_3 && 0 && N_4 && 0 \\
0 && N_1 && 0 && N_2 && 0 && N_3 && 0 && N_4 \\
\end{bmatrix}
\\
\boldsymbol{\epsilon} = \boldsymbol{B} \tilde{\boldsymbol{u}}

上の式を仮想仕事の式に代入すると、下記のように表すことができます。

\delta \boldsymbol{\epsilon} = \boldsymbol{B} \delta \tilde{\boldsymbol{u}} \\
\delta \boldsymbol{u} = \boldsymbol{N} \delta \tilde{\boldsymbol{u}} \\
\int_v \delta \boldsymbol{\epsilon}^T \boldsymbol{\sigma} dV = 
\int_v \delta \boldsymbol{u}^T \boldsymbol{b} dV + \int_s \delta \boldsymbol{u}^T \boldsymbol{t} dS \\
\int_v (\boldsymbol{B} \delta \tilde{\boldsymbol{u}} )^T \boldsymbol{D} \boldsymbol{B} \tilde{\boldsymbol{u}} dV = \int_v (\boldsymbol{N} \delta \tilde{\boldsymbol{u}})^T \boldsymbol{b} dV + \int_s (\boldsymbol{N} \delta \tilde{\boldsymbol{u}})^T \boldsymbol{t} dS \\ 
\delta \tilde{\boldsymbol{u}} : 要素の節点の仮想変位

ここで、$\delta \tilde{\boldsymbol{u}}$は要素内で一定値となるので、下記のように積分の外に出すことができます。

\delta \tilde{\boldsymbol{u}}^T \int_v \boldsymbol{B}^T \boldsymbol{D} \boldsymbol{B}  dV \tilde{\boldsymbol{u}} = \delta \tilde{\boldsymbol{u}}^T \int_v \boldsymbol{N}^T \boldsymbol{b} dV + \delta \tilde{\boldsymbol{u}}^T \int_s \boldsymbol{N}^T \boldsymbol{t} dS \\
\int_v \boldsymbol{B}^T \boldsymbol{D} \boldsymbol{B}  dV \tilde{\boldsymbol{u}} = \int_v \boldsymbol{N}^T \boldsymbol{b} dV + \int_s \boldsymbol{N}^T \boldsymbol{t} dS 

ここで、式を置き換えると下記のように表すことができます。

\int_v \boldsymbol{B}^T \boldsymbol{D} \boldsymbol{B} dV \tilde{\boldsymbol{u}} = \int_v \boldsymbol{N}^T \boldsymbol{b} dV + \int_s \boldsymbol{N}^T \boldsymbol{t} dS \\
\boldsymbol{K_e} \tilde{\boldsymbol{u}} = \boldsymbol{f_e} \\
\boldsymbol{K_e} = \int_v \boldsymbol{B}^T \boldsymbol{D} \boldsymbol{B} dV \\
\boldsymbol{f_e} = \int_v \boldsymbol{N}^T \boldsymbol{b} dV + \int_s \boldsymbol{N}^T \boldsymbol{t} dS \\
\boldsymbol{K_e} : 要素剛性マトリクス \\
\boldsymbol{f_e} : 等価節点力

上記より、要素剛性マトリクスは下記の式で求めることが可能です。

\boldsymbol{K_e} = \int_v \boldsymbol{B}^T \boldsymbol{D} \boldsymbol{B} dV 

ただ、上記のままでは計算しづらいので、正規座標系に変数変換します。

\begin{align}
\boldsymbol{K_e} &= \int_{0}^{1} \int_{0}^{1} \int_{0}^{1} \boldsymbol{B}^T\boldsymbol{D}\boldsymbol{B} |\boldsymbol{J}|dadbdc \\

\end{align}

正規化座標系では、各辺が0~1までの範囲のため上記のような積分範囲となります。そのため、ガウス求積を使用して積分を簡単な積の式に変換することができます。

\begin{align}
\boldsymbol{K_e} &= \int_v \boldsymbol{B}^T\boldsymbol{D}\boldsymbol{B} dV \\
&= \int_{0}^{1} \int_{0}^{1} \int_{0}^{1} \boldsymbol{B}^T\boldsymbol{D}\boldsymbol{B} |\boldsymbol{J}|dadbdc \\
&= w \boldsymbol{B}^T\boldsymbol{D}\boldsymbol{B} |\boldsymbol{J}|
\end{align}
\\
w = \frac{1}{6}: 重み係数

等価節点力

次に四面体要素の等価節点力を算出します。
等価節点力は物体力と表面力に別れているので、それぞれ求めてみます。

物体力

物体力はガウス求積を使うことで簡単に求めることが可能です。

\begin{align}
\int_v \boldsymbol{N}^T \boldsymbol{b} dV &= 
\int_{0}^{1} \int_{0}^{1} \int_{0}^{1} \boldsymbol{N}^T \boldsymbol{b} |\boldsymbol{J}| \\
&= w \boldsymbol{N}^T \boldsymbol{b} |\boldsymbol{J}| \\
\end{align}
\\
w = \frac{1}{6}: 重み係数

ここで、重力の場合は$\boldsymbol{b} = (密度)×(加速度ベクトル)$となります。

表面力

四面体要素は下記のように4面に分割することができます。

上記の図のように面番号を①~④に割り振り、各面の面積を$A_i$、各面に作用する表面力のベクトルを$\boldsymbol{p_i}$とおきます。
すると等価節点力の式は下記のように変形することができます。

\begin{align}
\int_s\boldsymbol{N}^T \boldsymbol{t} dS = &
\int_{s1} \boldsymbol{N}^T \boldsymbol{p_1}dS_1 + \int_{s2} \boldsymbol{N}^T \boldsymbol{p_2}dS_2 + \int_{s3} \boldsymbol{N}^T \boldsymbol{p_3}dS_3 + \int_{s4} \boldsymbol{N}^T \boldsymbol{p_4}dS_4 \\ 
=&
2A_1\int_0^1 \int_0^{1-a} \boldsymbol{N}^T(a,0,c) \boldsymbol{p_1}dcda + 2A_2\int_0^1 \int_0^{1-c} \boldsymbol{N}^T(0,b,c) \boldsymbol{p_2}dbdc + \\ 
& 2A_3 \int_0^1 \int_0^{1-b} \boldsymbol{N}^T(a,b,0) \boldsymbol{p_3}dadb + 2A_4 \int_0^1 \int_0^{1-c} \boldsymbol{N}^T(1-b-c,b,c) \boldsymbol{p_4}dbdc
\end{align}

ここで、$2A_i$は正規化座標系に変数変換するときのヤコビアンになります。
この式を計算すると下記のように表すことができます。

\int_s \boldsymbol{N}^T \boldsymbol{t} dS = 
\frac{A_1}{3}
\begin{bmatrix}
p_{1,x} \\
p_{1,y} \\
p_{1,z} \\
p_{1,x} \\
p_{1,y} \\
p_{1,z} \\
0 \\
0 \\
0 \\
p_{1,x} \\
p_{1,y} \\
p_{1,z}
\end{bmatrix}
+
\frac{A_2}{3}
\begin{bmatrix}
p_{2,x} \\
p_{2,y} \\
p_{2,z} \\
0 \\
0 \\
0 \\
p_{2,x} \\
p_{2,y} \\
p_{2,z} \\
p_{2,x} \\
p_{2,y} \\
p_{2,z}
\end{bmatrix}
+
\frac{A_3}{3}
\begin{bmatrix}
p_{3,x} \\
p_{3,y} \\
p_{3,z} \\
p_{3,x} \\
p_{3,y} \\
p_{3,z} \\
p_{3,x} \\
p_{3,y} \\
p_{3,z} \\
0 \\
0 \\
0 
\end{bmatrix}
+
\frac{A_4}{3}
\begin{bmatrix}
0 \\
0 \\
0 \\
p_{4,x} \\
p_{4,y} \\
p_{4,z} \\
p_{4,x} \\
p_{4,y} \\
p_{4,z} \\
p_{4,x} \\
p_{4,y} \\
p_{4,z}
\end{bmatrix}

Pythonで実装する

定式化が終わったのでPythonで要素の計算を行うクラスを実装してみます。
ここで、要素の計算を行うために、節点を格納するクラスNodeとDマトリクスの計算を行うクラスDmatrixを実装します。

Node.py

# 3次元節点格納用のクラス
class Node:

    # コンストラクタ
    # no : 節点番号
    # x  : x座標
    # y  : y座標
    # z  : z座標
    def __init__(self, no, x, y, z):
        self.no = no   # 節点番号
        self.x = x     # x座標
        self.y = y     # y座標
        self.z = z     # z座標

    # 節点の情報を表示する
    def printNode(self):
        print("Node No: %d, x: %f, y: %f, z: %f" % (self.no, self.x, self.y, self.z))

Dmatirx.py

import numpy as np

class Dmatrix:
    # コンストラクタ
    # young   : ヤング率
    # poisson : ポアソン比
    def __init__(self, young, poisson):
        self.young = young
        self.poisson = poisson

    # 弾性状態のDマトリクスを作成する
    def makeDematrix(self):

        tmp = self.young / ((1.0 + self.poisson) * (1.0 - 2.0 * self.poisson))
        matD = np.array([[1.0 - self.poisson, self.poisson, self.poisson, 0.0, 0.0, 0.0],
                         [self.poisson, 1.0 - self.poisson, self.poisson, 0.0, 0.0, 0.0],
                         [self.poisson, self.poisson, 1.0 - self.poisson, 0.0, 0.0, 0.0],
                         [0.0, 0.0, 0.0, 0.5 * (1.0 - 2.0 * self.poisson), 0.0, 0.0],
                         [0.0, 0.0, 0.0, 0.0, 0.5 * (1.0 - 2.0 * self.poisson), 0.0],
                         [0.0, 0.0, 0.0, 0.0, 0.0, 0.5 * (1.0 - 2.0 * self.poisson)]])
        matD = tmp * matD

        return matD

次に要素剛性マトリクス、等価節点力を計算するC3D4クラスを実装します。
今回は等価節点力は重力の項だけ実装します。

C3D4.py

import numpy as np
import numpy.linalg as LA
from Dmatrix import Dmatrix
from Node import Node

# 四面体4節点要素のクラス
class C3D4:
    # コンストラクタ
    # no           : 要素番号
    # nodes        : 節点の集合(Node型のリスト)
    # young        : ヤング率
    # poisson      : ポアソン比
    # density      : 密度
    # vecGravity   : 重力加速度のベクトル(np.array型)
    def __init__(self, no, nodes, young, poisson, density, vecGravity = None):

        # インスタンス変数を定義する
        self.nodeNum = 4               # 節点の数
        self.nodeDof = 3               # 節点の自由度
        self.no = no                   # 要素番号
        self.nodes = nodes             # nodesは反時計回りの順番になっている前提(Node2d型のリスト形式)
        self.young = young             # ヤング率
        self.poisson = poisson         # ポアソン比
        self.density = density         # 密度
        self.vecGravity = vecGravity   # 重力加速度のベクトル(np.array型)
        self.ipNum = 1                 # 積分点の数
        self.w = 1.0 / 6.0             # 積分点の重み係数
        self.ai = 1.0 / 4.0            # 積分点の座標(a,b,c座標系)
        self.bi = 1.0 / 4.0            # 積分点の座標(a,b,c座標系)
        self.ci = 1.0 / 4.0            # 積分点の座標(a,b,c座標系)

    # 要素剛性マトリクスKeを作成する
    def makeKematrix(self):

        # ヤコビ行列を計算する
        matJ = self.makeJmatrix()

        # Bマトリクスを計算する
        matB = self.makeBmatrix()

        # Dマトリクスを計算する
        matD = self.makeDmatrix()

        # Ketマトリクスをガウス積分で計算する
        matKet = self.w * matB.T @ matD @ matB * LA.det(matJ)

        return matKet

    # Dマトリクスを作成する
    def makeDmatrix(self):

        matD = Dmatrix(self.young, self.poisson).makeDematrix()
        return matD

    # ヤコビ行列を計算する
    def makeJmatrix(self):

        dxda = -self.nodes[0].x + self.nodes[1].x
        dyda = -self.nodes[0].y + self.nodes[1].y
        dzda = -self.nodes[0].z + self.nodes[1].z
        dxdb = -self.nodes[0].x + self.nodes[2].x
        dydb = -self.nodes[0].y + self.nodes[2].y
        dzdb = -self.nodes[0].z + self.nodes[2].z
        dxdc = -self.nodes[0].x + self.nodes[3].x
        dydc = -self.nodes[0].y + self.nodes[3].y
        dzdc = -self.nodes[0].z + self.nodes[3].z

        matJ = np.array([[dxda, dyda, dzda],
                         [dxdb, dydb, dzdb],
                         [dxdc, dydc, dzdc]])        

        # ヤコビアンが負にならないかチェックする
        if LA.det(matJ) < 0:
            raise ValueError("要素の計算に失敗しました")

        return matJ

    # Bマトリクスを作成する
    def makeBmatrix(self):

        # dNi/da, dNi/dbを計算する
        dN1da = -1.0
        dN2da = 1.0
        dN3da = 0.0
        dN4da = 0.0
        dN1db = -1.0
        dN2db = 0.0
        dN3db = 1.0
        dN4db = 0.0
        dN1dc = -1.0
        dN2dc = 0.0
        dN3dc = 0.0
        dN4dc = 1.0

        # dNi/dx, dNi/dyを計算する
        matdNdab = np.matrix([[dN1da, dN2da, dN3da, dN4da],
                              [dN1db, dN2db, dN3db, dN4db],
                              [dN1dc, dN2dc, dN3dc, dN4dc]])

         # ヤコビ行列を計算する
        matJ = self.makeJmatrix()

        #dNdxy = matJinv * matdNdab
        dNdxy = LA.solve(matJ, matdNdab)

        # Bマトリクスを計算する
        matB = np.array([[dNdxy[0, 0], 0.0, 0.0, dNdxy[0, 1], 0.0, 0.0, dNdxy[0, 2], 0.0, 0.0, dNdxy[0, 3], 0.0, 0.0],
                         [0.0, dNdxy[1, 0], 0.0, 0.0, dNdxy[1, 1], 0.0, 0.0, dNdxy[1, 2], 0.0, 0.0, dNdxy[1, 3], 0.0],
                         [0.0, 0.0, dNdxy[2, 0], 0.0, 0.0, dNdxy[2, 1], 0.0, 0.0, dNdxy[2, 2], 0.0, 0.0, dNdxy[2, 3]],
                         [0.0, dNdxy[2, 0], dNdxy[1, 0], 0.0, dNdxy[2, 1], dNdxy[1, 1], 0.0, dNdxy[2, 2], dNdxy[1, 2], 0.0, dNdxy[2, 3], dNdxy[1, 3]],
                         [dNdxy[2, 0], 0.0, dNdxy[0, 0], dNdxy[2, 1], 0.0, dNdxy[0, 1], dNdxy[2, 2], 0.0, dNdxy[0, 2], dNdxy[2, 3], 0.0, dNdxy[0, 3]],
                         [dNdxy[1, 0], dNdxy[0, 0], 0.0, dNdxy[1, 1], dNdxy[0, 1], 0.0, dNdxy[1, 2], dNdxy[0, 2], 0.0, dNdxy[1, 3], dNdxy[0, 3], 0.0]])

        return matB

    # 等価節点力の荷重ベクトルを作成する
    def makeEqNodeForceVector(self):

        # ヤコビ行列を計算する
        matJ = self.makeJmatrix()

        vecEqNodeForce = np.zeros(self.nodeNum * self.nodeDof)

        # 物体力による等価節点力を計算する
        vecBodyForce = np.zeros(self.nodeNum * self.nodeDof)
        if not self.vecGravity is None:
            vecb = self.density * self.vecGravity   # 単位体積あたりの物体力のベクトル
            N1 = 1 - self.ai - self.bi - self.ci
            N2 = self.ai
            N3 = self.bi
            N4 = self.ci
            matN = np.matrix([[N1, 0.0, 0.0, N2, 0.0, 0.0, N3, 0.0, 0.0, N4, 0.0, 0.0],
                              [0.0, N1, 0.0, 0.0, N2, 0.0, 0.0, N3, 0.0, 0.0, N4, 0.0],
                              [0.0, 0.0, N1, 0.0, 0.0, N2, 0.0, 0.0, N3, 0.0, 0.0, N4]])
            vecBodyForce = self.w * np.array(matN.T @ vecb).flatten() * LA.det(matJ)

        # 表面力と物体力の等価節点力を計算する
        vecEqNodeForce = vecBodyForce

        return vecEqNodeForce

剛性マトリクスの計算

上記より、要素剛性マトリクスが求まったので、それを足し合わせることによって、剛性マトリクスを計算することができます。
全体剛性マトリクスは3次元の場合、節点数を$N$とすると下記のように表されます。

\boldsymbol{K}=
\begin{bmatrix}
K_{11} && K_{12} && \cdots && K_{1, 3N} \\
K_{21} && K_{22} && \cdots && K_{2, 3N} \\
\vdots && \vdots && \ddots && \vdots \\
K_{3N, 1} && K_{3N, 2} && \cdots && K_{3N, 3N}
\end{bmatrix}
\begin{bmatrix}
u_1 \\
v_1 \\
\vdots \\
u_{N} \\
v_{N}
\end{bmatrix}

上記のように、$3N×3N$の行列になります。ここで、四面体要素の要素剛性マトリクスは節点数4個なので$12×12$の行列になります。そのため、全体剛性マトリクスに足し合わせるために対応する行列の位置を求める必要があります。
三角形要素の節点番号を$i$, $j$, $k$, $l$とすると、全体剛性マトリクスに格納する位置は下記のように表すことができます。

\boldsymbol{K}_e = 
\begin{bmatrix}
K_{3i-2, 3i-2} && K_{3i-2, 3i-1} && K_{3i-2, 3i} && \cdots && K_{3i-2, 3l-2} && K_{3i-2, 3l-1} && K_{3i-2, 3l} \\
K_{3i-1, 3i-2} && K_{3i-1, 3i-1} && K_{3i-1, 3i} && \cdots && K_{3i-1, 3l-2} && K_{3i-1, 3l-1} && K_{3i-1, 3l} \\
\vdots && \vdots && \vdots && \ddots && \vdots && \vdots && \vdots\\
K_{3l-1, 3i-2} && K_{3l-1, 3i-1} && K_{3l-1, 3i} && \cdots && K_{3l-1, 3l-2} && K_{3l-1, 3l-1} && K_{3l-1, 3l} \\
K_{3l, 3i-2} && K_{3l, 3i-1} && K_{3l, 3i} && \cdots && K_{3l, 3l-2} && K_{3l, 3l-1} && K_{3l, 3l}
\end{bmatrix}

境界条件の処理

剛性マトリクスが求まれば、いよいよ$\Delta \tilde{\boldsymbol{U}}$の計算に入れるわけですが、その前に境界条件を考慮する必要があります。
この境界条件を適用しないと剛性マトリクスは逆行列が無い状態なので、$\Delta \tilde{\boldsymbol{U}}$を計算することができません。
今回は$\boldsymbol{K} \Delta \tilde{\boldsymbol{U}}= \boldsymbol{R}$に単点拘束の境界条件を考慮した処理を行います。
単点拘束とは一つの節点自由度で構成される拘束のことで、例えば$u_1 = 0(u_1: 節点1の変位)$などの拘束を言います。
ここで、例として、剛性マトリクスが$3×3$で既知量$u_2 = \bar{u}$の場合で考えます。
すると、下記のように表すことができます。

\begin{bmatrix}
K_{11} && K_{12} && K_{13} \\
K_{21} && K_{22} && K_{23} \\
K_{31} && K_{32} && K_{33} \\
\end{bmatrix}
\begin{bmatrix}
u_1 \\
\bar{u} \\
u_3 \\
\end{bmatrix}
=
\begin{bmatrix}
r_1 \\
r_2 \\
r_3 \\
\end{bmatrix}
\\
K_{11} u_1 +  K_{12} \bar{u} + K_{13} u_3 = r_1 \\
K_{21} u_1 +  K_{22} \bar{u} + K_{23} u_3 = r_2 \\
K_{31} u_1 +  K_{32} \bar{u} + K_{33} u_3 = r_3 \\

既知量を右に移動すると

K_{11} u_1 + K_{13} u_3 = r_1 - K_{12} \bar{u} \\
K_{21} u_1 + K_{23} u_3 = r_2 - K_{22} \bar{u} \\
K_{31} u_1 + K_{33} u_3 = r_3 - K_{32} \bar{u}

となり、真ん中の式を省いて代わりに、$u_2 = \bar{u}$を代入すると、

K_{11} u_1 + K_{13} u_3 = r_1 - K_{12} \bar{u} \\
u_2 = \bar{u} \\
K_{31} u_1 + K_{33} u_3 = r_3 - K_{32} \bar{u}　\\
\begin{bmatrix}
K_{11} && 0 && K_{13} \\
0 && 1 && 0 \\
K_{31} && 0 && K_{33} 
\end{bmatrix}
\begin{bmatrix}
u_1 \\
\bar{u} \\
u_3 
\end{bmatrix}
=
\begin{bmatrix}
r_1 \\
\bar{u} \\
r_3
\end{bmatrix}
-
\begin{bmatrix}
K_{12}  \\
0 \\
K_{32}
\end{bmatrix}
\bar{u}

とすることができます。ここで、式を一つ省くことを行っていますが、一般に節点が拘束されている点の荷重値(反力)は未知数なので、未知の値$f_2$を行列から消して変位を求めることが可能となります。
このような処理をすることで、$\boldsymbol{K} \Delta \tilde{\boldsymbol{U}}= \boldsymbol{f}$の式を解くことが可能となります。

解析を行う処理をPythonで実装する

解析を行うクラスをPythonで実装します。
ここで、境界条件を格納するためのクラスBoundrayをまずは実装します。
Boundrayは基本的には境界条件を格納するだけのクラスです。
addSPC関数で単点拘束を追加、addForce関数で荷重を追加します。makeDispVector関数、makeForceVector関数は設定した拘束条件をベクトルの形で返します。

Boundary.py

import numpy as np

# 境界条件を格納するクラス
class Boundary:
    # コンストラクタ
    # nodeNum : 節点数
    def __init__(self, nodeNum):
        # インスタンス変数を定義する
        self.nodeNum = nodeNum                                     # 全節点数
        self.nodeDof = 3                                           # 節点の自由度
        self.vecDisp = np.array(nodeNum * self.nodeDof * [None])   # 単点拘束の強制変位
        self.vecForce = np.array(nodeNum * self.nodeDof * [0.0])   # 荷重ベクトル
        self.matC = np.empty((0, nodeNum * self.nodeDof))          # 多点拘束用のCマトリクス
        self.vecd = np.empty(0)                                    # 多点拘束用のdベクトル

    # 単点拘束を追加する
    # nodeNo : 節点番号
    # dispX  : x方向の強制変位
    # dispY  : y方向の強制変位
    # dispZ  : z方向の強制変位
    def addSPC(self, nodeNo, dispX, dispY, dispZ):

        self.vecDisp[self.nodeDof * (nodeNo - 1) + 0] = dispX
        self.vecDisp[self.nodeDof * (nodeNo - 1) + 1] = dispY
        self.vecDisp[self.nodeDof * (nodeNo - 1) + 2] = dispZ

    # 多点拘束を追加する
    # 条件式 : vecC x u = d
    def addMPC(self, vecC, d):

        self.matC = np.vstack((self.matC, vecC))
        self.vecd = np.hstack((self.vecd, d))

    # 単点拘束条件から変位ベクトルを作成する
    def makeDispVector(self):

        return self.vecDisp

    # 荷重を追加する
    def addForce(self, nodeNo, fx, fy, fz):

        self.vecForce[self.nodeDof * (nodeNo - 1) + 0] = fx
        self.vecForce[self.nodeDof * (nodeNo - 1) + 1] = fy
        self.vecForce[self.nodeDof * (nodeNo - 1) + 2] = fz

    # 境界条件から荷重ベクトルを作成する
    def makeForceVector(self):

        return self.vecForce

    # 多点拘束の境界条件を表すCマトリクス、dベクトルを作成する
    def makeMPCmatrixes(self):

        return self.matC, self.vecd

    # 拘束条件を出力する
    def printBoundary(self):
        print("Node Number: ", self.nodeNum)
        print("SPC Constraint Condition")
        print(self.vecDisp)
        print("Force Condition")
        print(self.vecForce)
        print("MPC Constraint Condition")
        print("C x u = d")
        print("C Matrix")
        print(self.matC)
        print("d vector")
        print(self.vecd)

次に解析を行うクラスFEMを実装します。
FEMクラスはコンストラクタで全節点、全要素、境界条件を引数にとります。
そして、analysis関数で解析を実行し、outputTxt関数で処理の結果をtxt形式で出力できます。

FEM.py

import numpy as np
import numpy.linalg as LA
from Boundary import Boundary

class FEM:
    # コンストラクタ
    # nodes    : 節点は1から始まる順番で並んでいる前提(Node型のリスト)
    # elements : 要素は種類ごとにソートされている前提(C3D4型のリスト)
    # bound    : 境界条件(d2Boundary型)
    def __init__(self, nodes, elements, bound):

        # インスタンス変数を定義する
        self.nodeDof = 3   # 節点の自由度
        self.nodes = nodes
        self.elements = elements
        self.bound = bound

    # 解析を行う
    def analysis(self):

        # 境界条件を考慮しないKマトリクスを作成する
        matK = self.makeKmatrix()

        # 荷重ベクトルを作成する
        vecf = self.makeForceVector()

        # 境界条件を考慮したKマトリクス、荷重ベクトルを作成する
        matKc, vecfc = self.setBoundCondition(matK, vecf)

        if np.isclose(LA.det(matKc), 0.0) :
            raise ValueError("有限要素法の計算に失敗しました。")

        # 変位ベクトルを計算する
        vecDisp = LA.solve(matKc, vecfc)
        self.vecDisp = vecDisp

        # 節点反力を計算する
        vecRF = np.array(matK @ vecDisp - vecf).flatten()
        self.vecRF = vecRF

        return vecDisp, vecRF

    # 節点に負荷する荷重、等価節点力を考慮した荷重ベクトルを作成する
    def makeForceVector(self):

        # 節点に負荷する荷重ベクトルを作成する
        vecCondiForce = self.bound.makeForceVector()

        # 等価節点力の荷重ベクトルを作成する
        vecEqNodeForce = np.zeros(len(self.nodes) * self.nodeDof)
        for elem in self.elements:
            vecElemEqNodeForce = elem.makeEqNodeForceVector()
            for i in range(len(elem.nodes)):
                for j in range(self.nodeDof):
                    vecEqNodeForce[self.nodeDof * (elem.nodes[i].no - 1) + j] += vecElemEqNodeForce[self.nodeDof * i + j]

        # 境界条件、等価節点力の荷重ベクトルを足し合わせる
        vecf = vecCondiForce + vecEqNodeForce

        return vecf

    # 境界条件を考慮しないKマトリクスを作成する
    def makeKmatrix(self):

        matK = np.matrix(np.zeros((len(self.nodes) * self.nodeDof, len(self.nodes) * self.nodeDof)))
        for elem in self.elements:

            # ketマトリクスを計算する
            matKe = elem.makeKematrix()

            # Ktマトリクスに代入する
            for c in range(len(elem.nodes) * self.nodeDof):
                ct = (elem.nodes[c // self.nodeDof].no - 1) * self.nodeDof + c % self.nodeDof
                for r in range(len(elem.nodes) * self.nodeDof):
                    rt = (elem.nodes[r // self.nodeDof].no - 1) * self.nodeDof + r % self.nodeDof
                    matK[ct, rt] += matKe[c, r]

        return matK

    # Kマトリクス、荷重ベクトルに境界条件を考慮する
    # matK         : 剛性マトリクス
    # vecf         : 荷重ベクトル
    # vecBoundDisp : 節点の境界条件の変位ベクトル
    # vecDisp      : 全節点の変位ベクトル(np.array型)
    def setBoundCondition(self, matKt, vecf):

        matKtc = np.copy(matKt)
        vecfc = np.copy(vecf)
        vecBoundDisp = self.bound.makeDispVector()

        # 単点拘束条件を考慮したKマトリクス、荷重ベクトルを作成する
        for i in range(len(vecBoundDisp)):
            if not vecBoundDisp[i] == None:
                # Kマトリクスからi列を抽出する
                vecx = np.array(matKt[:, i]).flatten()

                # 変位ベクトルi列の影響を荷重ベクトルに適用する
                vecfc = vecfc - (vecBoundDisp[i]) * vecx

                # Kマトリクスのi行、i列を全て0にし、i行i列の値を1にする
                matKtc[:, i] = 0.0
                matKtc[i, :] = 0.0
                matKtc[i, i] = 1.0
        for i in range(len(vecBoundDisp)):
            if not vecBoundDisp[i] == None:
                vecfc[i] = vecBoundDisp[i]

        return matKtc, vecfc

    # 解析結果をテキストファイルに出力する
    def outputTxt(self, filePath):

        # ファイルを作成し、開く
        f = open(filePath + ".txt", 'w')

        # 出力する文字の情報を定義する
        columNum = 20
        floatDigits = ".10g"

        # 入力データのタイトルを書きこむ
        f.write("*********************************\n")
        f.write("*          Input Data           *\n")
        f.write("*********************************\n")
        f.write("\n")

        # 節点情報を出力する
        f.write("***** Node Data ******\n")
        f.write("No".rjust(columNum) + "X".rjust(columNum) + "Y".rjust(columNum) + "Z".rjust(columNum) + "\n")
        f.write("-" * columNum * 4 + "\n")
        for node in self.nodes:
            strNo = str(node.no).rjust(columNum)
            strX = str(format(node.x, floatDigits).rjust(columNum))
            strY = str(format(node.y, floatDigits).rjust(columNum))
            strZ = str(format(node.z, floatDigits).rjust(columNum))
            f.write(strNo + strX + strY + strZ + "\n")
        f.write("\n")

        # 要素情報を出力する
        nodeNoColumNum = 36
        f.write("***** Element Data ******\n")
        f.write("No".rjust(columNum) + "Type".rjust(columNum) + "Node No".rjust(nodeNoColumNum) + 
                "Young".rjust(columNum) + "Poisson".rjust(columNum) + "Density".rjust(columNum) + "\n")
        f.write("-" * columNum * 5 + "-" * nodeNoColumNum + "\n")
        for elem in self.elements:
            strNo = str(elem.no).rjust(columNum)
            strType = str(elem.__class__.__name__ ).rjust(columNum)
            strNodeNo = ""
            for node in elem.nodes:
                strNodeNo += " " + str(node.no)
            strNodeNo = strNodeNo.rjust(nodeNoColumNum)
            strYoung = str(format(elem.young, floatDigits).rjust(columNum))
            strPoisson = "None".rjust(columNum)
            if hasattr(elem, 'poisson'):
                strPoisson = str(format(elem.poisson, floatDigits).rjust(columNum))
            strDensity = "None".rjust(columNum)
            if not elem.density is None:
                strDensity = str(format(elem.density, floatDigits).rjust(columNum))
            f.write(strNo + strType + strNodeNo + strYoung + strPoisson + strDensity + "\n")
        f.write("\n")

        # 単点拘束情報を出力する
        f.write("***** SPC Constraint Data ******\n")
        f.write("NodeNo".rjust(columNum) + "X Displacement".rjust(columNum) + "Y Displacement".rjust(columNum) + "Z Displacement".rjust(columNum) +"\n")
        f.write("-" * columNum * 4 + "\n")
        vecd = self.bound.makeDispVector()
        for i in range(len(self.nodes)):
            flg = False
            for j in range(self.nodeDof):
                if not vecd[self.nodeDof * i + j] == None:
                    flg = True
            if flg == True:
                strNo = str(i + 1).rjust(columNum)
                strXDisp = str(format(vecd[self.nodeDof * i], floatDigits).rjust(columNum))
                strYDisp = str(format(vecd[self.nodeDof * i + 1], floatDigits).rjust(columNum))
                strZDisp = str(format(vecd[self.nodeDof * i + 2], floatDigits).rjust(columNum))
                f.write(strNo + strXDisp + strYDisp + strZDisp + "\n")
        f.write("\n")

        # 荷重条件を出力する(等価節点力も含む)
        f.write("***** Nodal Force Data ******\n")
        f.write("NodeNo".rjust(columNum) + "X Force".rjust(columNum) + "Y Force".rjust(columNum) + "Z Force".rjust(columNum) +"\n")
        f.write("-" * columNum * 4 + "\n")
        vecf = self.makeForceVector()
        for i in range(len(self.nodes)):
            flg = False
            for j in range(self.nodeDof):
                if not vecf[self.nodeDof * i + j] == None:
                    flg = True
            if flg == True:
                strNo = str(i + 1).rjust(columNum)
                strXForce = str(format(vecf[self.nodeDof * i], floatDigits).rjust(columNum))
                strYForce = str(format(vecf[self.nodeDof * i + 1], floatDigits).rjust(columNum))
                strZForce = str(format(vecf[self.nodeDof * i + 2], floatDigits).rjust(columNum))
                f.write(strNo + strXForce + strYForce + strZForce + "\n")
        f.write("\n")

        # 結果データのタイトルを書きこむ
        f.write("**********************************\n")
        f.write("*          Result Data           *\n")
        f.write("**********************************\n")
        f.write("\n")

        # 変位のデータを出力する
        f.write("***** Displacement Data ******\n")
        f.write("NodeNo".rjust(columNum) + "Magnitude".rjust(columNum) + "X Displacement".rjust(columNum) +
                "Y Displacement".rjust(columNum) + "Z Displacement".rjust(columNum) + "\n")
        f.write("-" * columNum * 5 + "\n")
        for i in range(len(self.nodes)):
            strNo = str(i + 1).rjust(columNum)
            mag = np.linalg.norm(np.array((self.vecDisp[self.nodeDof * i], self.vecDisp[self.nodeDof * i + 1], self.vecDisp[self.nodeDof * i + 2])))
            strMag = str(format(mag, floatDigits).rjust(columNum))
            strXDisp = str(format(self.vecDisp[self.nodeDof * i], floatDigits).rjust(columNum))
            strYDisp = str(format(self.vecDisp[self.nodeDof * i + 1], floatDigits).rjust(columNum))
            strZDisp = str(format(self.vecDisp[self.nodeDof * i + 2], floatDigits).rjust(columNum))
            f.write(strNo + strMag + strXDisp + strYDisp + strZDisp + "\n")            
        f.write("\n")

        # 反力のデータを出力する
        f.write("***** Reaction Force Data ******\n")
        f.write("NodeNo".rjust(columNum) + "Magnitude".rjust(columNum) + "X Force".rjust(columNum) + "Y Force".rjust(columNum) + "Z Force".rjust(columNum) + "\n")
        f.write("-" * columNum * 5 + "\n")
        for i in range(len(self.nodes)):
            strNo = str(i + 1).rjust(columNum)
            mag = np.linalg.norm(np.array((self.vecRF[self.nodeDof * i], self.vecRF[self.nodeDof * i + 1], self.vecRF[self.nodeDof * i + 2])))
            strMag = str(format(mag, floatDigits).rjust(columNum))
            strXForce = str(format(self.vecRF[self.nodeDof * i], floatDigits).rjust(columNum))
            strYForce = str(format(self.vecRF[self.nodeDof * i + 1], floatDigits).rjust(columNum))
            strZForce = str(format(self.vecRF[self.nodeDof * i + 2], floatDigits).rjust(columNum))
            f.write(strNo + strMag + strXForce + strYForce + strZForce + "\n")            
        f.write("\n")

        # ファイルを閉じる
        f.close()

サンプルを解く

今回は下記のサンプルを解く処理を実装してみます。

メイン処理を実装する

メイン処理を実装してみます。
今回は読み込み機能を作らなかったので少々長くなってますが、基本的には入力して、analysis関数、outputTXT関数を実行しているだけです。

main.py

import numpy as np
from Node import Node
from C3D4 import C3D4
from Boundary import Boundary
from FEM import FEM

def main():

    # 節点を定義する
    node1=Node(1,0,-0.5,0)
    node2=Node(2,0,-0.443317,-0.231237)
    node3=Node(3,0,-0.287033,-0.409404)
    node4=Node(4,0,-0.066259,-0.49559)
    node5=Node(5,0,0.170342,-0.470089)
    node6=Node(6,0,0.371994,-0.334096)
    node7=Node(7,0,0.484751,-0.122543)
    node8=Node(8,0,0.485038,0.121402)
    node9=Node(9,0,0.37345,0.332468)
    node10=Node(10,0,0.17626,0.467902)
    node11=Node(11,0,-0.059932,0.496395)
    node12=Node(12,0,-0.283548,0.411826)
    node13=Node(13,0,-0.441818,0.234088)
    node14=Node(14,2,0.5,0)
    node15=Node(15,2,0.442529,-0.23274)
    node16=Node(16,2,0.284268,-0.411329)
    node17=Node(17,2,0.059082,-0.496497)
    node18=Node(18,2,-0.180417,-0.466315)
    node19=Node(19,2,-0.374624,-0.331144)
    node20=Node(20,2,-0.485183,-0.120822)
    node21=Node(21,2,-0.486407,0.115794)
    node22=Node(22,2,-0.378032,0.327249)
    node23=Node(23,2,-0.181838,0.465763)
    node24=Node(24,2,0.060095,0.496375)
    node25=Node(25,2,0.284233,0.411354)
    node26=Node(26,2,0.441545,0.234603)
    node27=Node(27,0.219345,-0.486987,-0.11333)
    node28=Node(28,0.219262,-0.180112,-0.466433)
    node29=Node(29,0.215686,0.281749,-0.413058)
    node30=Node(30,0.217934,0.5,-0.000152)
    node31=Node(31,0.218046,0.286409,0.409841)
    node32=Node(32,0.21857,-0.172211,0.469407)
    node33=Node(33,0.218914,-0.486701,0.114551)
    node34=Node(34,1.775186,0.486262,-0.116402)
    node35=Node(35,1.778671,0.174784,-0.468455)
    node36=Node(36,1.767065,-0.284339,-0.41128)
    node37=Node(37,1.772705,-0.499998,-0.001356)
    node38=Node(38,1.772127,-0.284499,0.41117)
    node39=Node(39,1.774461,0.170553,0.470012)
    node40=Node(40,1.770111,0.48506,0.121313)
    node41=Node(41,0.220275,-0.37675,-0.328723)
    node42=Node(42,0.21736,0.06268,-0.496056)
    node43=Node(43,0.216208,0.438843,-0.239618)
    node44=Node(44,0.217682,0.439287,0.238802)
    node45=Node(45,0.217807,0.055738,0.496884)
    node46=Node(46,0.219215,-0.374596,0.331177)
    node47=Node(47,1.771623,0.372089,-0.33399)
    node48=Node(48,1.77211,-0.063579,-0.495941)
    node49=Node(49,1.771088,-0.440641,-0.236296)
    node50=Node(50,1.772236,-0.441588,0.23452)
    node51=Node(51,1.774019,-0.05923,0.496479)
    node52=Node(52,1.771624,0.375264,0.330419)
    node53=Node(53,0.443641,-0.443765,-0.230374)
    node54=Node(54,0.433649,0.375589,-0.33005)
    node55=Node(55,0.438789,0.484293,0.124338)
    node56=Node(56,0.440291,0.175988,0.468005)
    node57=Node(57,0.440471,-0.286696,0.409641)
    node58=Node(58,1.549876,0.052894,-0.497194)
    node59=Node(59,1.551787,-0.371373,0.334786)
    node60=Node(60,1.542283,0.443115,-0.231623)
    node61=Node(61,1.544515,0.281706,-0.413088)
    node62=Node(62,1.542492,-0.485407,0.119916)
    node63=Node(63,1.550807,0.055197,0.496944)
    node64=Node(64,1.54722,0.443974,0.229972)
    node65=Node(65,0.433208,0.175573,-0.46816)
    node66=Node(66,1.541184,-0.485197,-0.120765)
    node67=Node(67,0.439717,-0.067249,-0.495457)
    node68=Node(68,0.439725,-0.499999,-0.000726)
    node69=Node(69,1.5443,0.499979,0.004533)
    node70=Node(70,1.548883,-0.372145,-0.333928)
    node71=Node(71,0.438245,0.484727,-0.122636)
    node72=Node(72,0.438095,0.372564,0.333461)
    node73=Node(73,0.437163,-0.057667,0.496663)
    node74=Node(74,0.442541,-0.443765,0.230375)
    node75=Node(75,0.443896,-0.287311,-0.409209)
    node76=Node(76,1.543119,0.283009,0.412196)
    node77=Node(77,1.553188,-0.179705,-0.46659)
    node78=Node(78,1.550653,-0.177171,0.467558)
    node79=Node(79,0.660277,-0.376071,-0.3295)
    node80=Node(80,0.657649,-0.373041,0.332927)
    node81=Node(81,1.340994,-0.283777,0.411668)
    node82=Node(82,0.655573,0.443299,0.23127)
    node83=Node(83,0.657914,0.058288,0.496591)
    node84=Node(84,1.336614,-0.286721,-0.409623)
    node85=Node(85,0.655531,0.44321,-0.23144)
    node86=Node(86,0.662015,-0.177323,0.4675)
    node87=Node(87,1.34018,-0.06555,-0.495685)
    node88=Node(88,0.661262,0.499998,0.001494)
    node89=Node(89,1.331122,0.374588,0.331185)
    node90=Node(90,1.329442,-0.499999,-0.000749)
    node91=Node(91,0.659764,0.053194,-0.497162)
    node92=Node(92,1.33323,0.176561,0.467789)
    node93=Node(93,0.658631,-0.485343,-0.120174)
    node94=Node(94,1.337409,-0.060666,0.496306)
    node95=Node(95,0.657166,0.285642,-0.410376)
    node96=Node(96,1.336979,-0.442835,0.232157)
    node97=Node(97,0.657073,-0.486696,0.114572)
    node98=Node(98,0.659412,-0.181359,-0.46595)
    node99=Node(99,0.655253,0.280314,0.414034)
    node100=Node(100,1.331314,0.37441,-0.331386)
    node101=Node(101,1.331121,-0.443209,-0.231442)
    node102=Node(102,1.334423,0.486047,-0.117298)
    node103=Node(103,1.342468,0.173141,-0.469065)
    node104=Node(104,1.330897,0.485077,0.121244)
    node105=Node(105,0.87575,-0.285777,-0.410282)
    node106=Node(106,0.879478,-0.063094,0.496003)
    node107=Node(107,1.103789,-0.184569,-0.464687)
    node108=Node(108,0.8798,0.375677,0.32995)
    node109=Node(109,0.874036,0.176788,0.467703)
    node110=Node(110,0.879604,-0.444914,-0.228147)
    node111=Node(111,1.116915,0.284625,0.411082)
    node112=Node(112,1.110382,-0.484984,-0.121616)
    node113=Node(113,0.88174,-0.44412,0.22969)
    node114=Node(114,0.879771,-0.066232,-0.495594)
    node115=Node(115,0.873607,0.486085,-0.117136)
    node116=Node(116,1.105422,-0.171848,0.46954)
    node117=Node(117,1.111634,-0.381301,-0.323434)
    node118=Node(118,0.879454,0.484141,0.124928)
    node119=Node(119,1.113052,-0.486369,0.115952)
    node120=Node(120,0.876175,-0.499977,-0.004841)
    node121=Node(121,1.112997,0.443324,-0.231223)
    node122=Node(122,1.111533,0.443535,0.230817)
    node123=Node(123,1.113221,0.061902,0.496153)
    node124=Node(124,1.10027,0.050111,-0.497483)
    node125=Node(125,1.109838,-0.37911,0.325999)
    node126=Node(126,0.876364,-0.282194,0.412755)
    node127=Node(127,0.878445,0.171812,-0.469554)
    node128=Node(128,1.109903,0.290881,-0.40668)
    node129=Node(129,0.871392,0.374921,-0.330809)
    node130=Node(130,1.112101,0.499919,0.008978)
    node131=Node(131,0,-0.26932,-0.071351)
    node132=Node(132,0,-0.15324,-0.253618)
    node133=Node(133,0,0.105699,-0.263277)
    node134=Node(134,0,0.276798,-0.063511)
    node135=Node(135,0,0.219089,0.184924)
    node136=Node(136,0,-0.031267,0.282196)
    node137=Node(137,0,-0.241217,0.157009)
    node138=Node(138,0,0.000335,0.000002)
    node139=Node(139,2,0.284315,-0.067786)
    node140=Node(140,2,0.160303,-0.241562)
    node141=Node(141,2,-0.105609,-0.259739)
    node142=Node(142,2,-0.270591,-0.065344)
    node143=Node(143,2,-0.211882,0.189364)
    node144=Node(144,2,0.037639,0.288793)
    node145=Node(145,2,0.263761,0.156163)
    node146=Node(146,2,0.019704,-0.005889)
    node147=Node(147,0.350091,0.06243,-0.082931)
    node148=Node(148,0.711505,-0.154336,-0.056662)
    node149=Node(149,1.482372,-0.113173,0.099924)
    node150=Node(150,1.100927,-0.144121,-0.030511)
    node151=Node(151,1.672873,-0.126232,-0.186987)
    node152=Node(152,0.428483,-0.070824,0.22419)
    node153=Node(153,1.662726,0.172345,0.10064)
    node154=Node(154,1.272937,0.214279,0.023313)
    node155=Node(155,0.86915,0.197796,0.049831)
    nodes=[node1,node2,node3,node4,node5,node6,node7,node8,node9,node10,
           node11,node12,node13,node14,node15,node16,node17,node18,node19,node20,
           node21,node22,node23,node24,node25,node26,node27,node28,node29,node30,
           node31,node32,node33,node34,node35,node36,node37,node38,node39,node40,
           node41,node42,node43,node44,node45,node46,node47,node48,node49,node50,
           node51,node52,node53,node54,node55,node56,node57,node58,node59,node60,
           node61,node62,node63,node64,node65,node66,node67,node68,node69,node70,
           node71,node72,node73,node74,node75,node76,node77,node78,node79,node80,
           node81,node82,node83,node84,node85,node86,node87,node88,node89,node90,
           node91,node92,node93,node94,node95,node96,node97,node98,node99,node100,
           node101,node102,node103,node104,node105,node106,node107,node108,node109,node110,
           node111,node112,node113,node114,node115,node116,node117,node118,node119,node120,
           node121,node122,node123,node124,node125,node126,node127,node128,node129,node130,
           node131,node132,node133,node134,node135,node136,node137,node138,node139,node140,
           node141,node142,node143,node144,node145,node146,node147,node148,node149,node150,
           node151,node152,node153,node154,node155]

    # 材料情報を定義する
    young = 210e9
    poisson = 0.3
    density = 7850.0
    vecGrav = np.array([0.0, 0.0, -9.81])

    # 要素を定義する
    elem1=C3D4(1,[node111, node89, node92, node154], young, poisson, density, vecGrav)
    elem2=C3D4(2,[node150,node96,node149,node81], young, poisson, density, vecGrav)
    elem3=C3D4(3,[node113,node150,node148,node126], young, poisson, density, vecGrav)
    elem4=C3D4(4,[node120,node113,node148,node97], young, poisson, density, vecGrav)
    elem5=C3D4(5,[node8,node30,node135,node134], young, poisson, density, vecGrav)
    elem6=C3D4(6,[node147,node133,node132,node138], young, poisson, density, vecGrav)
    elem7=C3D4(7,[node144,node153,node51,node39], young, poisson, density, vecGrav)
    elem8=C3D4(8,[node30,node8,node7,node134], young, poisson, density, vecGrav)
    elem9=C3D4(9,[node150,node106,node126,node116], young, poisson, density, vecGrav)
    elem10=C3D4(10,[node8,node44,node9,node135], young, poisson, density, vecGrav)
    elem11=C3D4(11,[node110,node120,node148,node93], young, poisson, density, vecGrav)
    elem12=C3D4(12,[node111,node155,node109,node108], young, poisson, density, vecGrav)
    elem13=C3D4(13,[node148,node68,node53,node93], young, poisson, density, vecGrav)
    elem14=C3D4(14,[node60,node154,node153,node69], young, poisson, density, vecGrav)
    elem15=C3D4(15,[node99,node155,node108,node109], young, poisson, density, vecGrav)
    elem16=C3D4(16,[node136,node45,node11,node32], young, poisson, density, vecGrav)
    elem17=C3D4(17,[node129,node155,node121,node115], young, poisson, density, vecGrav)
    elem18=C3D4(18,[node153,node89,node64,node154], young, poisson, density, vecGrav)
    elem19=C3D4(19,[node103,node100,node128,node154], young, poisson, density, vecGrav)
    elem20=C3D4(20,[node78,node81,node94,node149], young, poisson, density, vecGrav)
    elem21=C3D4(21,[node153,node89,node76,node64], young, poisson, density, vecGrav)
    elem22=C3D4(22,[node153,node52,node76,node39], young, poisson, density, vecGrav)
    elem23=C3D4(23,[node83,node73,node152,node86], young, poisson, density, vecGrav)
    elem24=C3D4(24,[node153,node89,node92,node76], young, poisson, density, vecGrav)
    elem25=C3D4(25,[node23,node144,node143,node51], young, poisson, density, vecGrav)
    elem26=C3D4(26,[node26,node52,node40,node145], young, poisson, density, vecGrav)
    elem27=C3D4(27,[node63,node153,node92,node76], young, poisson, density, vecGrav)
    elem28=C3D4(28,[node111,node122,node155,node108], young, poisson, density, vecGrav)
    elem29=C3D4(29,[node41,node147,node131,node27], young, poisson, density, vecGrav)
    elem30=C3D4(30,[node43,node29,node54,node147], young, poisson, density, vecGrav)
    elem31=C3D4(31,[node155,node106,node148,node150], young, poisson, density, vecGrav)
    elem32=C3D4(32,[node140,node17,node35,node16], young, poisson, density, vecGrav)
    elem33=C3D4(33,[node140,node15,node47,node139], young, poisson, density, vecGrav)
    elem34=C3D4(34,[node15,node34,node47,node139], young, poisson, density, vecGrav)
    elem35=C3D4(35,[node47,node151,node153,node61], young, poisson, density, vecGrav)
    elem36=C3D4(36,[node128,node155,node150,node154], young, poisson, density, vecGrav)
    elem37=C3D4(37,[node61,node151,node149,node103], young, poisson, density, vecGrav)
    elem38=C3D4(38,[node150,node125,node126,node113], young, poisson, density, vecGrav)
    elem39=C3D4(39,[node122,node155,node130,node154], young, poisson, density, vecGrav)
    elem40=C3D4(40,[node147,node91,node95,node148], young, poisson, density, vecGrav)
    elem41=C3D4(41,[node25,node144,node24,node39], young, poisson, density, vecGrav)
    elem42=C3D4(42,[node145,node39,node153,node144], young, poisson, density, vecGrav)
    elem43=C3D4(43,[node101,node90,node150,node112], young, poisson, density, vecGrav)
    elem44=C3D4(44,[node29,node133,node6,node5], young, poisson, density, vecGrav)
    elem45=C3D4(45,[node92,node111,node154,node123], young, poisson, density, vecGrav)
    elem46=C3D4(46,[node47,node151,node140,node153], young, poisson, density, vecGrav)
    elem47=C3D4(47,[node151,node77,node87,node58], young, poisson, density, vecGrav)
    elem48=C3D4(48,[node92,node154,node94,node123], young, poisson, density, vecGrav)
    elem49=C3D4(49,[node117,node101,node150,node112], young, poisson, density, vecGrav)
    elem50=C3D4(50,[node123,node155,node150,node106], young, poisson, density, vecGrav)
    elem51=C3D4(51,[node147,node136,node135,node138], young, poisson, density, vecGrav)
    elem52=C3D4(52,[node154,node150,node94,node123], young, poisson, density, vecGrav)
    elem53=C3D4(53,[node99,node155,node83,node152], young, poisson, density, vecGrav)
    elem54=C3D4(54,[node147,node91,node65,node95], young, poisson, density, vecGrav)
    elem55=C3D4(55,[node147,node91,node67,node65], young, poisson, density, vecGrav)
    elem56=C3D4(56,[node67,node147,node28,node75], young, poisson, density, vecGrav)
    elem57=C3D4(57,[node136,node32,node137,node152], young, poisson, density, vecGrav)
    elem58=C3D4(58,[node152,node82,node155,node147], young, poisson, density, vecGrav)
    elem59=C3D4(59,[node34,node153,node40,node69], young, poisson, density, vecGrav)
    elem60=C3D4(60,[node136,node147,node152,node137], young, poisson, density, vecGrav)
    elem61=C3D4(61,[node81,node96,node125,node150], young, poisson, density, vecGrav)
    elem62=C3D4(62,[node47,node151,node35,node140], young, poisson, density, vecGrav)
    elem63=C3D4(63,[node154,node150,node149,node94], young, poisson, density, vecGrav)
    elem64=C3D4(64,[node74,node80,node97,node152], young, poisson, density, vecGrav)
    elem65=C3D4(65,[node102,node121,node154,node130], young, poisson, density, vecGrav)
    elem66=C3D4(66,[node94,node150,node116,node123], young, poisson, density, vecGrav)
    elem67=C3D4(67,[node124,node103,node150,node87], young, poisson, density, vecGrav)
    elem68=C3D4(68,[node154,node122,node104,node130], young, poisson, density, vecGrav)
    elem69=C3D4(69,[node140,node48,node17,node141], young, poisson, density, vecGrav)
    elem70=C3D4(70,[node85,node95,node147,node54], young, poisson, density, vecGrav)
    elem71=C3D4(71,[node154,node103,node149,node150], young, poisson, density, vecGrav)
    elem72=C3D4(72,[node135,node44,node9,node31], young, poisson, density, vecGrav)
    elem73=C3D4(73,[node99,node155,node109,node83], young, poisson, density, vecGrav)
    elem74=C3D4(74,[node59,node78,node149,node81], young, poisson, density, vecGrav)
    elem75=C3D4(75,[node144,node153,node145,node146], young, poisson, density, vecGrav)
    elem76=C3D4(76,[node150,node106,node116,node123], young, poisson, density, vecGrav)
    elem77=C3D4(77,[node151,node143,node50,node142], young, poisson, density, vecGrav)
    elem78=C3D4(78,[node153,node139,node145,node146], young, poisson, density, vecGrav)
    elem79=C3D4(79,[node148,node147,node75,node53], young, poisson, density, vecGrav)
    elem80=C3D4(80,[node122,node118,node155,node108], young, poisson, density, vecGrav)
    elem81=C3D4(81,[node32,node12,node137,node46], young, poisson, density, vecGrav)
    elem82=C3D4(82,[node45,node136,node11,node10], young, poisson, density, vecGrav)
    elem83=C3D4(83,[node154,node61,node153,node149], young, poisson, density, vecGrav)
    elem84=C3D4(84,[node147,node148,node68,node53], young, poisson, density, vecGrav)
    elem85=C3D4(85,[node128,node103,node150,node124], young, poisson, density, vecGrav)
    elem86=C3D4(86,[node154,node60,node100,node102], young, poisson, density, vecGrav)
    elem87=C3D4(87,[node103,node128,node150,node154], young, poisson, density, vecGrav)
    elem88=C3D4(88,[node136,node45,node31,node10], young, poisson, density, vecGrav)
    elem89=C3D4(89,[node155,node121,node115,node130], young, poisson, density, vecGrav)
    elem90=C3D4(90,[node100,node121,node154,node102], young, poisson, density, vecGrav)
    elem91=C3D4(91,[node50,node149,node151,node62], young, poisson, density, vecGrav)
    elem92=C3D4(92,[node155,node121,node130,node154], young, poisson, density, vecGrav)
    elem93=C3D4(93,[node153,node34,node40,node139], young, poisson, density, vecGrav)
    elem94=C3D4(94,[node80,node74,node57,node152], young, poisson, density, vecGrav)
    elem95=C3D4(95,[node152,node72,node147,node31], young, poisson, density, vecGrav)
    elem96=C3D4(96,[node151,node37,node49,node142], young, poisson, density, vecGrav)
    elem97=C3D4(97,[node154,node60,node61,node100], young, poisson, density, vecGrav)
    elem98=C3D4(98,[node113,node80,node148,node97], young, poisson, density, vecGrav)
    elem99=C3D4(99,[node63,node153,node76,node39], young, poisson, density, vecGrav)
    elem100=C3D4(100,[node117,node110,node150,node105], young, poisson, density, vecGrav)
    elem101=C3D4(101,[node153,node89,node154,node92], young, poisson, density, vecGrav)
    elem102=C3D4(102,[node40,node153,node139,node145], young, poisson, density, vecGrav)
    elem103=C3D4(103,[node79,node148,node53,node93], young, poisson, density, vecGrav)
    elem104=C3D4(104,[node128,node155,node127,node150], young, poisson, density, vecGrav)
    elem105=C3D4(105,[node154,node155,node121,node128], young, poisson, density, vecGrav)
    elem106=C3D4(106,[node30,node43,node71,node147], young, poisson, density, vecGrav)
    elem107=C3D4(107,[node30,node8,node135,node44], young, poisson, density, vecGrav)
    elem108=C3D4(108,[node132,node41,node3,node28], young, poisson, density, vecGrav)
    elem109=C3D4(109,[node136,node135,node31,node147], young, poisson, density, vecGrav)
    elem110=C3D4(110,[node151,node36,node141,node49], young, poisson, density, vecGrav)
    elem111=C3D4(111,[node135,node44,node31,node147], young, poisson, density, vecGrav)
    elem112=C3D4(112,[node135,node136,node31,node10], young, poisson, density, vecGrav)
    elem113=C3D4(113,[node154,node60,node153,node61], young, poisson, density, vecGrav)
    elem114=C3D4(114,[node148,node68,node97,node152], young, poisson, density, vecGrav)
    elem115=C3D4(115,[node144,node23,node24,node51], young, poisson, density, vecGrav)
    elem116=C3D4(116,[node33,node147,node131,node137], young, poisson, density, vecGrav)
    elem117=C3D4(117,[node145,node39,node52,node153], young, poisson, density, vecGrav)
    elem118=C3D4(118,[node68,node74,node97,node152], young, poisson, density, vecGrav)
    elem119=C3D4(119,[node37,node21,node142,node50], young, poisson, density, vecGrav)
    elem120=C3D4(120,[node41,node147,node27,node53], young, poisson, density, vecGrav)
    elem121=C3D4(121,[node133,node43,node6,node134], young, poisson, density, vecGrav)
    elem122=C3D4(122,[node142,node37,node49,node20], young, poisson, density, vecGrav)
    elem123=C3D4(123,[node48,node151,node58,node77], young, poisson, density, vecGrav)
    elem124=C3D4(124,[node118,node155,node88,node115], young, poisson, density, vecGrav)
    elem125=C3D4(125,[node21,node143,node142,node50], young, poisson, density, vecGrav)
    elem126=C3D4(126,[node155,node106,node83,node148], young, poisson, density, vecGrav)
    elem127=C3D4(127,[node84,node149,node151,node87], young, poisson, density, vecGrav)
    elem128=C3D4(128,[node100,node121,node128,node154], young, poisson, density, vecGrav)
    elem129=C3D4(129,[node85,node54,node147,node71], young, poisson, density, vecGrav)
    elem130=C3D4(130,[node153,node34,node60,node69], young, poisson, density, vecGrav)
    elem131=C3D4(131,[node152,node82,node99,node155], young, poisson, density, vecGrav)
    elem132=C3D4(132,[node153,node47,node61,node60], young, poisson, density, vecGrav)
    elem133=C3D4(133,[node30,node135,node134,node147], young, poisson, density, vecGrav)
    elem134=C3D4(134,[node148,node110,node79,node105], young, poisson, density, vecGrav)
    elem135=C3D4(135,[node47,node151,node61,node35], young, poisson, density, vecGrav)
    elem136=C3D4(136,[node70,node151,node49,node36], young, poisson, density, vecGrav)
    elem137=C3D4(137,[node149,node66,node151,node62], young, poisson, density, vecGrav)
    elem138=C3D4(138,[node19,node36,node49,node141], young, poisson, density, vecGrav)
    elem139=C3D4(139,[node114,node150,node124,node127], young, poisson, density, vecGrav)
    elem140=C3D4(140,[node55,node30,node71,node147], young, poisson, density, vecGrav)
    elem141=C3D4(141,[node150,node114,node105,node148], young, poisson, density, vecGrav)
    elem142=C3D4(142,[node129,node85,node155,node115], young, poisson, density, vecGrav)
    elem143=C3D4(143,[node97,node120,node93,node148], young, poisson, density, vecGrav)
    elem144=C3D4(144,[node150,node155,node127,node148], young, poisson, density, vecGrav)
    elem145=C3D4(145,[node51,node63,node153,node149], young, poisson, density, vecGrav)
    elem146=C3D4(146,[node153,node140,node146,node151], young, poisson, density, vecGrav)
    elem147=C3D4(147,[node151,node103,node87,node149], young, poisson, density, vecGrav)
    elem148=C3D4(148,[node45,node73,node32,node152], young, poisson, density, vecGrav)
    elem149=C3D4(149,[node57,node80,node152,node86], young, poisson, density, vecGrav)
    elem150=C3D4(150,[node21,node37,node142,node20], young, poisson, density, vecGrav)
    elem151=C3D4(151,[node87,node103,node150,node149], young, poisson, density, vecGrav)
    elem152=C3D4(152,[node96,node59,node149,node81], young, poisson, density, vecGrav)
    elem153=C3D4(153,[node98,node148,node79,node105], young, poisson, density, vecGrav)
    elem154=C3D4(154,[node151,node48,node36,node77], young, poisson, density, vecGrav)
    elem155=C3D4(155,[node10,node135,node9,node31], young, poisson, density, vecGrav)
    elem156=C3D4(156,[node117,node110,node112,node150], young, poisson, density, vecGrav)
    elem157=C3D4(157,[node133,node29,node42,node5], young, poisson, density, vecGrav)
    elem158=C3D4(158,[node154,node61,node103,node100], young, poisson, density, vecGrav)
    elem159=C3D4(159,[node13,node33,node137,node46], young, poisson, density, vecGrav)
    elem160=C3D4(160,[node96,node59,node62,node149], young, poisson, density, vecGrav)
    elem161=C3D4(161,[node103,node151,node58,node61], young, poisson, density, vecGrav)
    elem162=C3D4(162,[node50,node149,node143,node151], young, poisson, density, vecGrav)
    elem163=C3D4(163,[node33,node13,node137,node1], young, poisson, density, vecGrav)
    elem164=C3D4(164,[node149,node92,node63,node153], young, poisson, density, vecGrav)
    elem165=C3D4(165,[node43,node30,node134,node147], young, poisson, density, vecGrav)
    elem166=C3D4(166,[node84,node149,node101,node151], young, poisson, density, vecGrav)
    elem167=C3D4(167,[node133,node43,node134,node147], young, poisson, density, vecGrav)
    elem168=C3D4(168,[node7,node30,node134,node43], young, poisson, density, vecGrav)
    elem169=C3D4(169,[node77,node151,node87,node84], young, poisson, density, vecGrav)
    elem170=C3D4(170,[node63,node78,node94,node149], young, poisson, density, vecGrav)
    elem171=C3D4(171,[node140,node48,node35,node17], young, poisson, density, vecGrav)
    elem172=C3D4(172,[node149,node90,node150,node101], young, poisson, density, vecGrav)
    elem173=C3D4(173,[node153,node52,node64,node76], young, poisson, density, vecGrav)
    elem174=C3D4(174,[node149,node63,node78,node51], young, poisson, density, vecGrav)
    elem175=C3D4(175,[node153,node52,node40,node64], young, poisson, density, vecGrav)
    elem176=C3D4(176,[node52,node26,node25,node145], young, poisson, density, vecGrav)
    elem177=C3D4(177,[node149,node92,node94,node63], young, poisson, density, vecGrav)
    elem178=C3D4(178,[node145,node39,node25,node52], young, poisson, density, vecGrav)
    elem179=C3D4(179,[node48,node151,node36,node141], young, poisson, density, vecGrav)
    elem180=C3D4(180,[node149,node50,node143,node38], young, poisson, density, vecGrav)
    elem181=C3D4(181,[node131,node33,node137,node1], young, poisson, density, vecGrav)
    elem182=C3D4(182,[node145,node39,node144,node25], young, poisson, density, vecGrav)
    elem183=C3D4(183,[node143,node151,node146,node142], young, poisson, density, vecGrav)
    elem184=C3D4(184,[node86,node80,node148,node126], young, poisson, density, vecGrav)
    elem185=C3D4(185,[node117,node84,node150,node101], young, poisson, density, vecGrav)
    elem186=C3D4(186,[node144,node51,node24,node39], young, poisson, density, vecGrav)
    elem187=C3D4(187,[node114,node150,node127,node148], young, poisson, density, vecGrav)
    elem188=C3D4(188,[node69,node154,node102,node60], young, poisson, density, vecGrav)
    elem189=C3D4(189,[node153,node140,node47,node139], young, poisson, density, vecGrav)
    elem190=C3D4(190,[node143,node151,node149,node153], young, poisson, density, vecGrav)
    elem191=C3D4(191,[node128,node155,node129,node127], young, poisson, density, vecGrav)
    elem192=C3D4(192,[node106,node86,node148,node126], young, poisson, density, vecGrav)
    elem193=C3D4(193,[node153,node151,node146,node143], young, poisson, density, vecGrav)
    elem194=C3D4(194,[node110,node120,node112,node150], young, poisson, density, vecGrav)
    elem195=C3D4(195,[node136,node45,node32,node152], young, poisson, density, vecGrav)
    elem196=C3D4(196,[node111,node155,node122,node154], young, poisson, density, vecGrav)
    elem197=C3D4(197,[node89,node111,node122,node154], young, poisson, density, vecGrav)
    elem198=C3D4(198,[node33,node147,node68,node27], young, poisson, density, vecGrav)
    elem199=C3D4(199,[node78,node51,node38,node149], young, poisson, density, vecGrav)
    elem200=C3D4(200,[node56,node45,node31,node152], young, poisson, density, vecGrav)
    elem201=C3D4(201,[node150,node120,node112,node119], young, poisson, density, vecGrav)
    elem202=C3D4(202,[node69,node154,node153,node64], young, poisson, density, vecGrav)
    elem203=C3D4(203,[node139,node40,node145,node14], young, poisson, density, vecGrav)
    elem204=C3D4(204,[node153,node143,node51,node149], young, poisson, density, vecGrav)
    elem205=C3D4(205,[node45,node136,node31,node152], young, poisson, density, vecGrav)
    elem206=C3D4(206,[node147,node133,node29,node42], young, poisson, density, vecGrav)
    elem207=C3D4(207,[node47,node140,node35,node16], young, poisson, density, vecGrav)
    elem208=C3D4(208,[node143,node22,node38,node50], young, poisson, density, vecGrav)
    elem209=C3D4(209,[node52,node153,node40,node145], young, poisson, density, vecGrav)
    elem210=C3D4(210,[node23,node143,node38,node51], young, poisson, density, vecGrav)
    elem211=C3D4(211,[node81,node125,node116,node150], young, poisson, density, vecGrav)
    elem212=C3D4(212,[node140,node48,node151,node35], young, poisson, density, vecGrav)
    elem213=C3D4(213,[node55,node44,node147,node72], young, poisson, density, vecGrav)
    elem214=C3D4(214,[node40,node26,node145,node14], young, poisson, density, vecGrav)
    elem215=C3D4(215,[node81,node150,node116,node94], young, poisson, density, vecGrav)
    elem216=C3D4(216,[node155,node118,node82,node108], young, poisson, density, vecGrav)
    elem217=C3D4(217,[node34,node153,node47,node139], young, poisson, density, vecGrav)
    elem218=C3D4(218,[node32,node46,node137,node152], young, poisson, density, vecGrav)
    elem219=C3D4(219,[node107,node114,node105,node150], young, poisson, density, vecGrav)
    elem220=C3D4(220,[node4,node133,node42,node5], young, poisson, density, vecGrav)
    elem221=C3D4(221,[node147,node91,node148,node67], young, poisson, density, vecGrav)
    elem222=C3D4(222,[node148,node147,node68,node152], young, poisson, density, vecGrav)
    elem223=C3D4(223,[node151,node48,node58,node35], young, poisson, density, vecGrav)
    elem224=C3D4(224,[node4,node132,node28,node42], young, poisson, density, vecGrav)
    elem225=C3D4(225,[node67,node148,node98,node91], young, poisson, density, vecGrav)
    elem226=C3D4(226,[node152,node72,node56,node99], young, poisson, density, vecGrav)
    elem227=C3D4(227,[node147,node41,node75,node53], young, poisson, density, vecGrav)
    elem228=C3D4(228,[node118,node155,node115,node130], young, poisson, density, vecGrav)
    elem229=C3D4(229,[node150,node128,node124,node127], young, poisson, density, vecGrav)
    elem230=C3D4(230,[node147,node44,node31,node72], young, poisson, density, vecGrav)
    elem231=C3D4(231,[node155,node88,node147,node82], young, poisson, density, vecGrav)
    elem232=C3D4(232,[node150,node113,node148,node120], young, poisson, density, vecGrav)
    elem233=C3D4(233,[node4,node133,node132,node42], young, poisson, density, vecGrav)
    elem234=C3D4(234,[node155,node85,node88,node115], young, poisson, density, vecGrav)
    elem235=C3D4(235,[node51,node143,node38,node149], young, poisson, density, vecGrav)
    elem236=C3D4(236,[node143,node22,node50,node21], young, poisson, density, vecGrav)
    elem237=C3D4(237,[node151,node70,node84,node77], young, poisson, density, vecGrav)
    elem238=C3D4(238,[node132,node147,node28,node42], young, poisson, density, vecGrav)
    elem239=C3D4(239,[node152,node72,node31,node56], young, poisson, density, vecGrav)
    elem240=C3D4(240,[node46,node32,node57,node152], young, poisson, density, vecGrav)
    elem241=C3D4(241,[node133,node147,node132,node42], young, poisson, density, vecGrav)
    elem242=C3D4(242,[node123,node154,node150,node155], young, poisson, density, vecGrav)
    elem243=C3D4(243,[node84,node149,node150,node101], young, poisson, density, vecGrav)
    elem244=C3D4(244,[node147,node136,node138,node137], young, poisson, density, vecGrav)
    elem245=C3D4(245,[node56,node99,node83,node152], young, poisson, density, vecGrav)
    elem246=C3D4(246,[node151,node37,node62,node66], young, poisson, density, vecGrav)
    elem247=C3D4(247,[node147,node41,node131,node132], young, poisson, density, vecGrav)
    elem248=C3D4(248,[node131,node147,node138,node137], young, poisson, density, vecGrav)
    elem249=C3D4(249,[node99,node155,node82,node108], young, poisson, density, vecGrav)
    elem250=C3D4(250,[node152,node82,node72,node99], young, poisson, density, vecGrav)
    elem251=C3D4(251,[node2,node131,node1,node27], young, poisson, density, vecGrav)
    elem252=C3D4(252,[node88,node85,node147,node71], young, poisson, density, vecGrav)
    elem253=C3D4(253,[node114,node107,node124,node150], young, poisson, density, vecGrav)
    elem254=C3D4(254,[node95,node155,node148,node127], young, poisson, density, vecGrav)
    elem255=C3D4(255,[node147,node68,node27,node53], young, poisson, density, vecGrav)
    elem256=C3D4(256,[node69,node154,node64,node104], young, poisson, density, vecGrav)
    elem257=C3D4(257,[node147,node132,node131,node138], young, poisson, density, vecGrav)
    elem258=C3D4(258,[node147,node55,node72,node82], young, poisson, density, vecGrav)
    elem259=C3D4(259,[node131,node33,node1,node27], young, poisson, density, vecGrav)
    elem260=C3D4(260,[node96,node62,node90,node149], young, poisson, density, vecGrav)
    elem261=C3D4(261,[node37,node151,node49,node66], young, poisson, density, vecGrav)
    elem262=C3D4(262,[node146,node151,node141,node142], young, poisson, density, vecGrav)
    elem263=C3D4(263,[node83,node86,node152,node148], young, poisson, density, vecGrav)
    elem264=C3D4(264,[node30,node44,node147,node55], young, poisson, density, vecGrav)
    elem265=C3D4(265,[node15,node140,node47,node16], young, poisson, density, vecGrav)
    elem266=C3D4(266,[node102,node154,node104,node130], young, poisson, density, vecGrav)
    elem267=C3D4(267,[node154,node61,node149,node103], young, poisson, density, vecGrav)
    elem268=C3D4(268,[node23,node22,node38,node143], young, poisson, density, vecGrav)
    elem269=C3D4(269,[node38,node50,node59,node149], young, poisson, density, vecGrav)
    elem270=C3D4(270,[node153,node144,node143,node146], young, poisson, density, vecGrav)
    elem271=C3D4(271,[node147,node33,node152,node137], young, poisson, density, vecGrav)
    elem272=C3D4(272,[node68,node148,node97,node93], young, poisson, density, vecGrav)
    elem273=C3D4(273,[node147,node67,node148,node75], young, poisson, density, vecGrav)
    elem274=C3D4(274,[node41,node131,node2,node27], young, poisson, density, vecGrav)
    elem275=C3D4(275,[node58,node151,node35,node61], young, poisson, density, vecGrav)
    elem276=C3D4(276,[node153,node63,node51,node39], young, poisson, density, vecGrav)
    elem277=C3D4(277,[node73,node56,node83,node152], young, poisson, density, vecGrav)
    elem278=C3D4(278,[node149,node66,node101,node151], young, poisson, density, vecGrav)
    elem279=C3D4(279,[node74,node46,node57,node152], young, poisson, density, vecGrav)
    elem280=C3D4(280,[node30,node44,node135,node147], young, poisson, density, vecGrav)
    elem281=C3D4(281,[node19,node142,node49,node20], young, poisson, density, vecGrav)
    elem282=C3D4(282,[node88,node55,node147,node82], young, poisson, density, vecGrav)
    elem283=C3D4(283,[node41,node147,node28,node132], young, poisson, density, vecGrav)
    elem284=C3D4(284,[node118,node155,node82,node88], young, poisson, density, vecGrav)
    elem285=C3D4(285,[node41,node132,node2,node131], young, poisson, density, vecGrav)
    elem286=C3D4(286,[node33,node46,node152,node137], young, poisson, density, vecGrav)
    elem287=C3D4(287,[node114,node148,node98,node105], young, poisson, density, vecGrav)
    elem288=C3D4(288,[node132,node41,node2,node3], young, poisson, density, vecGrav)
    elem289=C3D4(289,[node150,node106,node148,node126], young, poisson, density, vecGrav)
    elem290=C3D4(290,[node73,node57,node152,node86], young, poisson, density, vecGrav)
    elem291=C3D4(291,[node155,node83,node152,node148], young, poisson, density, vecGrav)
    elem292=C3D4(292,[node150,node125,node113,node119], young, poisson, density, vecGrav)
    elem293=C3D4(293,[node4,node132,node3,node28], young, poisson, density, vecGrav)
    elem294=C3D4(294,[node125,node150,node126,node116], young, poisson, density, vecGrav)
    elem295=C3D4(295,[node150,node81,node149,node94], young, poisson, density, vecGrav)
    elem296=C3D4(296,[node147,node41,node28,node75], young, poisson, density, vecGrav)
    elem297=C3D4(297,[node69,node154,node104,node102], young, poisson, density, vecGrav)
    elem298=C3D4(298,[node155,node106,node109,node83], young, poisson, density, vecGrav)
    elem299=C3D4(299,[node18,node48,node141,node17], young, poisson, density, vecGrav)
    elem300=C3D4(300,[node45,node56,node73,node152], young, poisson, density, vecGrav)
    elem301=C3D4(301,[node147,node155,node152,node148], young, poisson, density, vecGrav)
    elem302=C3D4(302,[node152,node82,node147,node72], young, poisson, density, vecGrav)
    elem303=C3D4(303,[node148,node67,node98,node75], young, poisson, density, vecGrav)
    elem304=C3D4(304,[node155,node154,node111,node123], young, poisson, density, vecGrav)
    elem305=C3D4(305,[node80,node86,node148,node152], young, poisson, density, vecGrav)
    elem306=C3D4(306,[node114,node91,node148,node127], young, poisson, density, vecGrav)
    elem307=C3D4(307,[node90,node112,node119,node150], young, poisson, density, vecGrav)
    elem308=C3D4(308,[node123,node155,node109,node111], young, poisson, density, vecGrav)
    elem309=C3D4(309,[node96,node125,node150,node119], young, poisson, density, vecGrav)
    elem310=C3D4(310,[node110,node148,node79,node93], young, poisson, density, vecGrav)
    elem311=C3D4(311,[node149,node92,node154,node94], young, poisson, density, vecGrav)
    elem312=C3D4(312,[node124,node107,node87,node150], young, poisson, density, vecGrav)
    elem313=C3D4(313,[node151,node70,node66,node101], young, poisson, density, vecGrav)
    elem314=C3D4(314,[node149,node66,node90,node101], young, poisson, density, vecGrav)
    elem315=C3D4(315,[node70,node151,node84,node101], young, poisson, density, vecGrav)
    elem316=C3D4(316,[node128,node155,node121,node129], young, poisson, density, vecGrav)
    elem317=C3D4(317,[node33,node74,node68,node152], young, poisson, density, vecGrav)
    elem318=C3D4(318,[node29,node147,node42,node65], young, poisson, density, vecGrav)
    elem319=C3D4(319,[node148,node98,node79,node75], young, poisson, density, vecGrav)
    elem320=C3D4(320,[node34,node153,node60,node47], young, poisson, density, vecGrav)
    elem321=C3D4(321,[node97,node80,node148,node152], young, poisson, density, vecGrav)
    elem322=C3D4(322,[node12,node13,node137,node46], young, poisson, density, vecGrav)
    elem323=C3D4(323,[node147,node67,node42,node65], young, poisson, density, vecGrav)
    elem324=C3D4(324,[node67,node147,node42,node28], young, poisson, density, vecGrav)
    elem325=C3D4(325,[node148,node79,node53,node75], young, poisson, density, vecGrav)
    elem326=C3D4(326,[node136,node12,node137,node32], young, poisson, density, vecGrav)
    elem327=C3D4(327,[node142,node151,node141,node49], young, poisson, density, vecGrav)
    elem328=C3D4(328,[node34,node15,node14,node139], young, poisson, density, vecGrav)
    elem329=C3D4(329,[node151,node153,node61,node149], young, poisson, density, vecGrav)
    elem330=C3D4(330,[node133,node147,node134,node138], young, poisson, density, vecGrav)
    elem331=C3D4(331,[node74,node33,node46,node152], young, poisson, density, vecGrav)
    elem332=C3D4(332,[node43,node7,node6,node134], young, poisson, density, vecGrav)
    elem333=C3D4(333,[node147,node33,node131,node27], young, poisson, density, vecGrav)
    elem334=C3D4(334,[node103,node151,node87,node58], young, poisson, density, vecGrav)
    elem335=C3D4(335,[node133,node29,node43,node147], young, poisson, density, vecGrav)
    elem336=C3D4(336,[node136,node147,node31,node152], young, poisson, density, vecGrav)
    elem337=C3D4(337,[node151,node70,node49,node66], young, poisson, density, vecGrav)
    elem338=C3D4(338,[node19,node18,node36,node141], young, poisson, density, vecGrav)
    elem339=C3D4(339,[node147,node135,node134,node138], young, poisson, density, vecGrav)
    elem340=C3D4(340,[node151,node140,node146,node141], young, poisson, density, vecGrav)
    elem341=C3D4(341,[node140,node48,node141,node151], young, poisson, density, vecGrav)
    elem342=C3D4(342,[node89,node122,node104,node154], young, poisson, density, vecGrav)
    elem343=C3D4(343,[node90,node96,node150,node119], young, poisson, density, vecGrav)
    elem344=C3D4(344,[node89,node154,node104,node64], young, poisson, density, vecGrav)
    elem345=C3D4(345,[node118,node122,node155,node130], young, poisson, density, vecGrav)
    elem346=C3D4(346,[node32,node73,node57,node152], young, poisson, density, vecGrav)
    elem347=C3D4(347,[node91,node95,node148,node127], young, poisson, density, vecGrav)
    elem348=C3D4(348,[node85,node155,node147,node95], young, poisson, density, vecGrav)
    elem349=C3D4(349,[node70,node151,node36,node77], young, poisson, density, vecGrav)
    elem350=C3D4(350,[node153,node140,node139,node146], young, poisson, density, vecGrav)
    elem351=C3D4(351,[node18,node48,node36,node141], young, poisson, density, vecGrav)
    elem352=C3D4(352,[node43,node54,node71,node147], young, poisson, density, vecGrav)
    elem353=C3D4(353,[node29,node147,node65,node54], young, poisson, density, vecGrav)
    elem354=C3D4(354,[node12,node136,node11,node32], young, poisson, density, vecGrav)
    elem355=C3D4(355,[node85,node129,node155,node95], young, poisson, density, vecGrav)
    elem356=C3D4(356,[node133,node29,node6,node43], young, poisson, density, vecGrav)
    elem357=C3D4(357,[node129,node155,node95,node127], young, poisson, density, vecGrav)
    elem358=C3D4(358,[node110,node150,node148,node120], young, poisson, density, vecGrav)
    elem359=C3D4(359,[node153,node40,node69,node64], young, poisson, density, vecGrav)
    elem360=C3D4(360,[node34,node40,node139,node14], young, poisson, density, vecGrav)
    elem361=C3D4(361,[node155,node85,node147,node88], young, poisson, density, vecGrav)
    elem362=C3D4(362,[node80,node113,node148,node126], young, poisson, density, vecGrav)
    elem363=C3D4(363,[node107,node117,node150,node105], young, poisson, density, vecGrav)
    elem364=C3D4(364,[node123,node155,node106,node109], young, poisson, density, vecGrav)
    elem365=C3D4(365,[node149,node92,node153,node154], young, poisson, density, vecGrav)
    elem366=C3D4(366,[node153,node144,node51,node143], young, poisson, density, vecGrav)
    elem367=C3D4(367,[node95,node54,node65,node147], young, poisson, density, vecGrav)
    elem368=C3D4(368,[node55,node88,node147,node71], young, poisson, density, vecGrav)
    elem369=C3D4(369,[node110,node150,node105,node148], young, poisson, density, vecGrav)
    elem370=C3D4(370,[node59,node78,node38,node149], young, poisson, density, vecGrav)
    elem371=C3D4(371,[node50,node149,node62,node59], young, poisson, density, vecGrav)
    elem372=C3D4(372,[node151,node37,node50,node62], young, poisson, density, vecGrav)
    elem373=C3D4(373,[node83,node106,node86,node148], young, poisson, density, vecGrav)
    elem374=C3D4(374,[node149,node96,node150,node90], young, poisson, density, vecGrav)
    elem375=C3D4(375,[node149,node66,node62,node90], young, poisson, density, vecGrav)
    elem376=C3D4(376,[node107,node84,node150,node117], young, poisson, density, vecGrav)
    elem377=C3D4(377,[node37,node151,node50,node142], young, poisson, density, vecGrav)
    elem378=C3D4(378,[node84,node149,node87,node150], young, poisson, density, vecGrav)
    elem379=C3D4(379,[node150,node113,node120,node119], young, poisson, density, vecGrav)
    elem380=C3D4(380,[node84,node107,node150,node87], young, poisson, density, vecGrav)
    elem381=C3D4(381,[node142,node19,node49,node141], young, poisson, density, vecGrav)
    elem382=C3D4(382,[node95,node155,node147,node148], young, poisson, density, vecGrav)
    elem383=C3D4(383,[node148,node114,node98,node91], young, poisson, density, vecGrav)
    elem384=C3D4(384,[node147,node33,node68,node152], young, poisson, density, vecGrav)
    elems=[elem1,elem2,elem3,elem4,elem5,elem6,elem7,elem8,elem9,elem10,
           elem11,elem12,elem13,elem14,elem15,elem16,elem17,elem18,elem19,elem20,
           elem21,elem22,elem23,elem24,elem25,elem26,elem27,elem28,elem29,elem30,
           elem31,elem32,elem33,elem34,elem35,elem36,elem37,elem38,elem39,elem40,
           elem41,elem42,elem43,elem44,elem45,elem46,elem47,elem48,elem49,elem50,
           elem51,elem52,elem53,elem54,elem55,elem56,elem57,elem58,elem59,elem60,
           elem61,elem62,elem63,elem64,elem65,elem66,elem67,elem68,elem69,elem70,
           elem71,elem72,elem73,elem74,elem75,elem76,elem77,elem78,elem79,elem80,
           elem81,elem82,elem83,elem84,elem85,elem86,elem87,elem88,elem89,elem90,
           elem91,elem92,elem93,elem94,elem95,elem96,elem97,elem98,elem99,elem100,
           elem101,elem102,elem103,elem104,elem105,elem106,elem107,elem108,elem109,elem110,
           elem111,elem112,elem113,elem114,elem115,elem116,elem117,elem118,elem119,elem120,
           elem121,elem122,elem123,elem124,elem125,elem126,elem127,elem128,elem129,elem130,
           elem131,elem132,elem133,elem134,elem135,elem136,elem137,elem138,elem139,elem140,
           elem141,elem142,elem143,elem144,elem145,elem146,elem147,elem148,elem149,elem150,
           elem151,elem152,elem153,elem154,elem155,elem156,elem157,elem158,elem159,elem160,
           elem161,elem162,elem163,elem164,elem165,elem166,elem167,elem168,elem169,elem170,
           elem171,elem172,elem173,elem174,elem175,elem176,elem177,elem178,elem179,elem180,
           elem181,elem182,elem183,elem184,elem185,elem186,elem187,elem188,elem189,elem190,
           elem191,elem192,elem193,elem194,elem195,elem196,elem197,elem198,elem199,elem200,
           elem201,elem202,elem203,elem204,elem205,elem206,elem207,elem208,elem209,elem210,
           elem211,elem212,elem213,elem214,elem215,elem216,elem217,elem218,elem219,elem220,
           elem221,elem222,elem223,elem224,elem225,elem226,elem227,elem228,elem229,elem230,
           elem231,elem232,elem233,elem234,elem235,elem236,elem237,elem238,elem239,elem240,
           elem241,elem242,elem243,elem244,elem245,elem246,elem247,elem248,elem249,elem250,
           elem251,elem252,elem253,elem254,elem255,elem256,elem257,elem258,elem259,elem260,
           elem261,elem262,elem263,elem264,elem265,elem266,elem267,elem268,elem269,elem270,
           elem271,elem272,elem273,elem274,elem275,elem276,elem277,elem278,elem279,elem280,
           elem281,elem282,elem283,elem284,elem285,elem286,elem287,elem288,elem289,elem290,
           elem291,elem292,elem293,elem294,elem295,elem296,elem297,elem298,elem299,elem300,
           elem301,elem302,elem303,elem304,elem305,elem306,elem307,elem308,elem309,elem310,
           elem311,elem312,elem313,elem314,elem315,elem316,elem317,elem318,elem319,elem320,
           elem321,elem322,elem323,elem324,elem325,elem326,elem327,elem328,elem329,elem330,
           elem331,elem332,elem333,elem334,elem335,elem336,elem337,elem338,elem339,elem340,
           elem341,elem342,elem343,elem344,elem345,elem346,elem347,elem348,elem349,elem350,
           elem351,elem352,elem353,elem354,elem355,elem356,elem357,elem358,elem359,elem360,
           elem361,elem362,elem363,elem364,elem365,elem366,elem367,elem368,elem369,elem370,
           elem371,elem372,elem373,elem374,elem375,elem376,elem377,elem378,elem379,elem380,
           elem381,elem382,elem383,elem384]

    # 境界条件を設定する
    bound = Boundary(len(nodes))
    bound.addSPC(1, 0.0, 0.0, 0.0)
    bound.addSPC(2, 0.0, 0.0, 0.0)
    bound.addSPC(3, 0.0, 0.0, 0.0)
    bound.addSPC(4, 0.0, 0.0, 0.0)
    bound.addSPC(5, 0.0, 0.0, 0.0)
    bound.addSPC(6, 0.0, 0.0, 0.0)
    bound.addSPC(7, 0.0, 0.0, 0.0)
    bound.addSPC(8, 0.0, 0.0, 0.0)
    bound.addSPC(9, 0.0, 0.0, 0.0)
    bound.addSPC(10, 0.0, 0.0, 0.0)
    bound.addSPC(11, 0.0, 0.0, 0.0)
    bound.addSPC(12, 0.0, 0.0, 0.0)
    bound.addSPC(13, 0.0, 0.0, 0.0)
    bound.addSPC(131, 0.0, 0.0, 0.0)
    bound.addSPC(132, 0.0, 0.0, 0.0)
    bound.addSPC(133, 0.0, 0.0, 0.0)
    bound.addSPC(134, 0.0, 0.0, 0.0)
    bound.addSPC(135, 0.0, 0.0, 0.0)
    bound.addSPC(136, 0.0, 0.0, 0.0)
    bound.addSPC(137, 0.0, 0.0, 0.0)
    bound.addSPC(138, 0.0, 0.0, 0.0)
    bound.addForce(146, 0.0, 0.0, -1e9)

    # 解析を行う
    fem = FEM(nodes, elems, bound)
    fem.analysis()
    fem.outputTxt("test")

main()

Abaqus2021 Student Edition2021と比較する

出力された結果が正しいか確認するためにAbaqusと比較してみます。

Abaqusの入力ファイル

*node
1, 0, -0.5, 0
2, 0, -0.443317, -0.231237
3, 0, -0.287033, -0.409404
4, 0, -0.066259, -0.49559
5, 0, 0.170342, -0.470089
6, 0, 0.371994, -0.334096
7, 0, 0.484751, -0.122543
8, 0, 0.485038, 0.121402
9, 0, 0.37345, 0.332468
10, 0, 0.17626, 0.467902
11, 0, -0.059932, 0.496395
12, 0, -0.283548, 0.411826
13, 0, -0.441818, 0.234088
14, 2, 0.5, 0
15, 2, 0.442529, -0.23274
16, 2, 0.284268, -0.411329
17, 2, 0.059082, -0.496497
18, 2, -0.180417, -0.466315
19, 2, -0.374624, -0.331144
20, 2, -0.485183, -0.120822
21, 2, -0.486407, 0.115794
22, 2, -0.378032, 0.327249
23, 2, -0.181838, 0.465763
24, 2, 0.060095, 0.496375
25, 2, 0.284233, 0.411354
26, 2, 0.441545, 0.234603
27, 0.219345, -0.486987, -0.11333
28, 0.219262, -0.180112, -0.466433
29, 0.215686, 0.281749, -0.413058
30, 0.217934, 0.5, -0.000152
31, 0.218046, 0.286409, 0.409841
32, 0.21857, -0.172211, 0.469407
33, 0.218914, -0.486701, 0.114551
34, 1.775186, 0.486262, -0.116402
35, 1.778671, 0.174784, -0.468455
36, 1.767065, -0.284339, -0.41128
37, 1.772705, -0.499998, -0.001356
38, 1.772127, -0.284499, 0.41117
39, 1.774461, 0.170553, 0.470012
40, 1.770111, 0.48506, 0.121313
41, 0.220275, -0.37675, -0.328723
42, 0.21736, 0.06268, -0.496056
43, 0.216208, 0.438843, -0.239618
44, 0.217682, 0.439287, 0.238802
45, 0.217807, 0.055738, 0.496884
46, 0.219215, -0.374596, 0.331177
47, 1.771623, 0.372089, -0.33399
48, 1.77211, -0.063579, -0.495941
49, 1.771088, -0.440641, -0.236296
50, 1.772236, -0.441588, 0.23452
51, 1.774019, -0.05923, 0.496479
52, 1.771624, 0.375264, 0.330419
53, 0.443641, -0.443765, -0.230374
54, 0.433649, 0.375589, -0.33005
55, 0.438789, 0.484293, 0.124338
56, 0.440291, 0.175988, 0.468005
57, 0.440471, -0.286696, 0.409641
58, 1.549876, 0.052894, -0.497194
59, 1.551787, -0.371373, 0.334786
60, 1.542283, 0.443115, -0.231623
61, 1.544515, 0.281706, -0.413088
62, 1.542492, -0.485407, 0.119916
63, 1.550807, 0.055197, 0.496944
64, 1.54722, 0.443974, 0.229972
65, 0.433208, 0.175573, -0.46816
66, 1.541184, -0.485197, -0.120765
67, 0.439717, -0.067249, -0.495457
68, 0.439725, -0.499999, -0.000726
69, 1.5443, 0.499979, 0.004533
70, 1.548883, -0.372145, -0.333928
71, 0.438245, 0.484727, -0.122636
72, 0.438095, 0.372564, 0.333461
73, 0.437163, -0.057667, 0.496663
74, 0.442541, -0.443765, 0.230375
75, 0.443896, -0.287311, -0.409209
76, 1.543119, 0.283009, 0.412196
77, 1.553188, -0.179705, -0.46659
78, 1.550653, -0.177171, 0.467558
79, 0.660277, -0.376071, -0.3295
80, 0.657649, -0.373041, 0.332927
81, 1.340994, -0.283777, 0.411668
82, 0.655573, 0.443299, 0.23127
83, 0.657914, 0.058288, 0.496591
84, 1.336614, -0.286721, -0.409623
85, 0.655531, 0.44321, -0.23144
86, 0.662015, -0.177323, 0.4675
87, 1.34018, -0.06555, -0.495685
88, 0.661262, 0.499998, 0.001494
89, 1.331122, 0.374588, 0.331185
90, 1.329442, -0.499999, -0.000749
91, 0.659764, 0.053194, -0.497162
92, 1.33323, 0.176561, 0.467789
93, 0.658631, -0.485343, -0.120174
94, 1.337409, -0.060666, 0.496306
95, 0.657166, 0.285642, -0.410376
96, 1.336979, -0.442835, 0.232157
97, 0.657073, -0.486696, 0.114572
98, 0.659412, -0.181359, -0.46595
99, 0.655253, 0.280314, 0.414034
100, 1.331314, 0.37441, -0.331386
101, 1.331121, -0.443209, -0.231442
102, 1.334423, 0.486047, -0.117298
103, 1.342468, 0.173141, -0.469065
104, 1.330897, 0.485077, 0.121244
105, 0.87575, -0.285777, -0.410282
106, 0.879478, -0.063094, 0.496003
107, 1.103789, -0.184569, -0.464687
108, 0.8798, 0.375677, 0.32995
109, 0.874036, 0.176788, 0.467703
110, 0.879604, -0.444914, -0.228147
111, 1.116915, 0.284625, 0.411082
112, 1.110382, -0.484984, -0.121616
113, 0.88174, -0.44412, 0.22969
114, 0.879771, -0.066232, -0.495594
115, 0.873607, 0.486085, -0.117136
116, 1.105422, -0.171848, 0.46954
117, 1.111634, -0.381301, -0.323434
118, 0.879454, 0.484141, 0.124928
119, 1.113052, -0.486369, 0.115952
120, 0.876175, -0.499977, -0.004841
121, 1.112997, 0.443324, -0.231223
122, 1.111533, 0.443535, 0.230817
123, 1.113221, 0.061902, 0.496153
124, 1.10027, 0.050111, -0.497483
125, 1.109838, -0.37911, 0.325999
126, 0.876364, -0.282194, 0.412755
127, 0.878445, 0.171812, -0.469554
128, 1.109903, 0.290881, -0.40668
129, 0.871392, 0.374921, -0.330809
130, 1.112101, 0.499919, 0.008978
131, 0, -0.26932, -0.071351
132, 0, -0.15324, -0.253618
133, 0, 0.105699, -0.263277
134, 0, 0.276798, -0.063511
135, 0, 0.219089, 0.184924
136, 0, -0.031267, 0.282196
137, 0, -0.241217, 0.157009
138, 0, 0.000335, 0.000002
139, 2, 0.284315, -0.067786
140, 2, 0.160303, -0.241562
141, 2, -0.105609, -0.259739
142, 2, -0.270591, -0.065344
143, 2, -0.211882, 0.189364
144, 2, 0.037639, 0.288793
145, 2, 0.263761, 0.156163
146, 2, 0.019704, -0.005889
147, 0.350091, 0.06243, -0.082931
148, 0.711505, -0.154336, -0.056662
149, 1.482372, -0.113173, 0.099924
150, 1.100927, -0.144121, -0.030511
151, 1.672873, -0.126232, -0.186987
152, 0.428483, -0.070824, 0.22419
153, 1.662726, 0.172345, 0.10064
154, 1.272937, 0.214279, 0.023313
155, 0.86915, 0.197796, 0.049831
*element, type=C3D4, elset=elems
1, 111, 89, 92, 154
2, 150, 96, 149, 81
3, 113, 150, 148, 126
4, 120, 113, 148, 97
5, 8, 30, 135, 134
6, 147, 133, 132, 138
7, 144, 153, 51, 39
8, 30, 8, 7, 134
9, 150, 106, 126, 116
10, 8, 44, 9, 135
11, 110, 120, 148, 93
12, 111, 155, 109, 108
13, 148, 68, 53, 93
14, 60, 154, 153, 69
15, 99, 155, 108, 109
16, 136, 45, 11, 32
17, 129, 155, 121, 115
18, 153, 89, 64, 154
19, 103, 100, 128, 154
20, 78, 81, 94, 149
21, 153, 89, 76, 64
22, 153, 52, 76, 39
23, 83, 73, 152, 86
24, 153, 89, 92, 76
25, 23, 144, 143, 51
26, 26, 52, 40, 145
27, 63, 153, 92, 76
28, 111, 122, 155, 108
29, 41, 147, 131, 27
30, 43, 29, 54, 147
31, 155, 106, 148, 150
32, 140, 17, 35, 16
33, 140, 15, 47, 139
34, 15, 34, 47, 139
35, 47, 151, 153, 61
36, 128, 155, 150, 154
37, 61, 151, 149, 103
38, 150, 125, 126, 113
39, 122, 155, 130, 154
40, 147, 91, 95, 148
41, 25, 144, 24, 39
42, 145, 39, 153, 144
43, 101, 90, 150, 112
44, 29, 133, 6, 5
45, 92, 111, 154, 123
46, 47, 151, 140, 153
47, 151, 77, 87, 58
48, 92, 154, 94, 123
49, 117, 101, 150, 112
50, 123, 155, 150, 106
51, 147, 136, 135, 138
52, 154, 150, 94, 123
53, 99, 155, 83, 152
54, 147, 91, 65, 95
55, 147, 91, 67, 65
56, 67, 147, 28, 75
57, 136, 32, 137, 152
58, 152, 82, 155, 147
59, 34, 153, 40, 69
60, 136, 147, 152, 137
61, 81, 96, 125, 150
62, 47, 151, 35, 140
63, 154, 150, 149, 94
64, 74, 80, 97, 152
65, 102, 121, 154, 130
66, 94, 150, 116, 123
67, 124, 103, 150, 87
68, 154, 122, 104, 130
69, 140, 48, 17, 141
70, 85, 95, 147, 54
71, 154, 103, 149, 150
72, 135, 44, 9, 31
73, 99, 155, 109, 83
74, 59, 78, 149, 81
75, 144, 153, 145, 146
76, 150, 106, 116, 123
77, 151, 143, 50, 142
78, 153, 139, 145, 146
79, 148, 147, 75, 53
80, 122, 118, 155, 108
81, 32, 12, 137, 46
82, 45, 136, 11, 10
83, 154, 61, 153, 149
84, 147, 148, 68, 53
85, 128, 103, 150, 124
86, 154, 60, 100, 102
87, 103, 128, 150, 154
88, 136, 45, 31, 10
89, 155, 121, 115, 130
90, 100, 121, 154, 102
91, 50, 149, 151, 62
92, 155, 121, 130, 154
93, 153, 34, 40, 139
94, 80, 74, 57, 152
95, 152, 72, 147, 31
96, 151, 37, 49, 142
97, 154, 60, 61, 100
98, 113, 80, 148, 97
99, 63, 153, 76, 39
100, 117, 110, 150, 105
101, 153, 89, 154, 92
102, 40, 153, 139, 145
103, 79, 148, 53, 93
104, 128, 155, 127, 150
105, 154, 155, 121, 128
106, 30, 43, 71, 147
107, 30, 8, 135, 44
108, 132, 41, 3, 28
109, 136, 135, 31, 147
110, 151, 36, 141, 49
111, 135, 44, 31, 147
112, 135, 136, 31, 10
113, 154, 60, 153, 61
114, 148, 68, 97, 152
115, 144, 23, 24, 51
116, 33, 147, 131, 137
117, 145, 39, 52, 153
118, 68, 74, 97, 152
119, 37, 21, 142, 50
120, 41, 147, 27, 53
121, 133, 43, 6, 134
122, 142, 37, 49, 20
123, 48, 151, 58, 77
124, 118, 155, 88, 115
125, 21, 143, 142, 50
126, 155, 106, 83, 148
127, 84, 149, 151, 87
128, 100, 121, 128, 154
129, 85, 54, 147, 71
130, 153, 34, 60, 69
131, 152, 82, 99, 155
132, 153, 47, 61, 60
133, 30, 135, 134, 147
134, 148, 110, 79, 105
135, 47, 151, 61, 35
136, 70, 151, 49, 36
137, 149, 66, 151, 62
138, 19, 36, 49, 141
139, 114, 150, 124, 127
140, 55, 30, 71, 147
141, 150, 114, 105, 148
142, 129, 85, 155, 115
143, 97, 120, 93, 148
144, 150, 155, 127, 148
145, 51, 63, 153, 149
146, 153, 140, 146, 151
147, 151, 103, 87, 149
148, 45, 73, 32, 152
149, 57, 80, 152, 86
150, 21, 37, 142, 20
151, 87, 103, 150, 149
152, 96, 59, 149, 81
153, 98, 148, 79, 105
154, 151, 48, 36, 77
155, 10, 135, 9, 31
156, 117, 110, 112, 150
157, 133, 29, 42, 5
158, 154, 61, 103, 100
159, 13, 33, 137, 46
160, 96, 59, 62, 149
161, 103, 151, 58, 61
162, 50, 149, 143, 151
163, 33, 13, 137, 1
164, 149, 92, 63, 153
165, 43, 30, 134, 147
166, 84, 149, 101, 151
167, 133, 43, 134, 147
168, 7, 30, 134, 43
169, 77, 151, 87, 84
170, 63, 78, 94, 149
171, 140, 48, 35, 17
172, 149, 90, 150, 101
173, 153, 52, 64, 76
174, 149, 63, 78, 51
175, 153, 52, 40, 64
176, 52, 26, 25, 145
177, 149, 92, 94, 63
178, 145, 39, 25, 52
179, 48, 151, 36, 141
180, 149, 50, 143, 38
181, 131, 33, 137, 1
182, 145, 39, 144, 25
183, 143, 151, 146, 142
184, 86, 80, 148, 126
185, 117, 84, 150, 101
186, 144, 51, 24, 39
187, 114, 150, 127, 148
188, 69, 154, 102, 60
189, 153, 140, 47, 139
190, 143, 151, 149, 153
191, 128, 155, 129, 127
192, 106, 86, 148, 126
193, 153, 151, 146, 143
194, 110, 120, 112, 150
195, 136, 45, 32, 152
196, 111, 155, 122, 154
197, 89, 111, 122, 154
198, 33, 147, 68, 27
199, 78, 51, 38, 149
200, 56, 45, 31, 152
201, 150, 120, 112, 119
202, 69, 154, 153, 64
203, 139, 40, 145, 14
204, 153, 143, 51, 149
205, 45, 136, 31, 152
206, 147, 133, 29, 42
207, 47, 140, 35, 16
208, 143, 22, 38, 50
209, 52, 153, 40, 145
210, 23, 143, 38, 51
211, 81, 125, 116, 150
212, 140, 48, 151, 35
213, 55, 44, 147, 72
214, 40, 26, 145, 14
215, 81, 150, 116, 94
216, 155, 118, 82, 108
217, 34, 153, 47, 139
218, 32, 46, 137, 152
219, 107, 114, 105, 150
220, 4, 133, 42, 5
221, 147, 91, 148, 67
222, 148, 147, 68, 152
223, 151, 48, 58, 35
224, 4, 132, 28, 42
225, 67, 148, 98, 91
226, 152, 72, 56, 99
227, 147, 41, 75, 53
228, 118, 155, 115, 130
229, 150, 128, 124, 127
230, 147, 44, 31, 72
231, 155, 88, 147, 82
232, 150, 113, 148, 120
233, 4, 133, 132, 42
234, 155, 85, 88, 115
235, 51, 143, 38, 149
236, 143, 22, 50, 21
237, 151, 70, 84, 77
238, 132, 147, 28, 42
239, 152, 72, 31, 56
240, 46, 32, 57, 152
241, 133, 147, 132, 42
242, 123, 154, 150, 155
243, 84, 149, 150, 101
244, 147, 136, 138, 137
245, 56, 99, 83, 152
246, 151, 37, 62, 66
247, 147, 41, 131, 132
248, 131, 147, 138, 137
249, 99, 155, 82, 108
250, 152, 82, 72, 99
251, 2, 131, 1, 27
252, 88, 85, 147, 71
253, 114, 107, 124, 150
254, 95, 155, 148, 127
255, 147, 68, 27, 53
256, 69, 154, 64, 104
257, 147, 132, 131, 138
258, 147, 55, 72, 82
259, 131, 33, 1, 27
260, 96, 62, 90, 149
261, 37, 151, 49, 66
262, 146, 151, 141, 142
263, 83, 86, 152, 148
264, 30, 44, 147, 55
265, 15, 140, 47, 16
266, 102, 154, 104, 130
267, 154, 61, 149, 103
268, 23, 22, 38, 143
269, 38, 50, 59, 149
270, 153, 144, 143, 146
271, 147, 33, 152, 137
272, 68, 148, 97, 93
273, 147, 67, 148, 75
274, 41, 131, 2, 27
275, 58, 151, 35, 61
276, 153, 63, 51, 39
277, 73, 56, 83, 152
278, 149, 66, 101, 151
279, 74, 46, 57, 152
280, 30, 44, 135, 147
281, 19, 142, 49, 20
282, 88, 55, 147, 82
283, 41, 147, 28, 132
284, 118, 155, 82, 88
285, 41, 132, 2, 131
286, 33, 46, 152, 137
287, 114, 148, 98, 105
288, 132, 41, 2, 3
289, 150, 106, 148, 126
290, 73, 57, 152, 86
291, 155, 83, 152, 148
292, 150, 125, 113, 119
293, 4, 132, 3, 28
294, 125, 150, 126, 116
295, 150, 81, 149, 94
296, 147, 41, 28, 75
297, 69, 154, 104, 102
298, 155, 106, 109, 83
299, 18, 48, 141, 17
300, 45, 56, 73, 152
301, 147, 155, 152, 148
302, 152, 82, 147, 72
303, 148, 67, 98, 75
304, 155, 154, 111, 123
305, 80, 86, 148, 152
306, 114, 91, 148, 127
307, 90, 112, 119, 150
308, 123, 155, 109, 111
309, 96, 125, 150, 119
310, 110, 148, 79, 93
311, 149, 92, 154, 94
312, 124, 107, 87, 150
313, 151, 70, 66, 101
314, 149, 66, 90, 101
315, 70, 151, 84, 101
316, 128, 155, 121, 129
317, 33, 74, 68, 152
318, 29, 147, 42, 65
319, 148, 98, 79, 75
320, 34, 153, 60, 47
321, 97, 80, 148, 152
322, 12, 13, 137, 46
323, 147, 67, 42, 65
324, 67, 147, 42, 28
325, 148, 79, 53, 75
326, 136, 12, 137, 32
327, 142, 151, 141, 49
328, 34, 15, 14, 139
329, 151, 153, 61, 149
330, 133, 147, 134, 138
331, 74, 33, 46, 152
332, 43, 7, 6, 134
333, 147, 33, 131, 27
334, 103, 151, 87, 58
335, 133, 29, 43, 147
336, 136, 147, 31, 152
337, 151, 70, 49, 66
338, 19, 18, 36, 141
339, 147, 135, 134, 138
340, 151, 140, 146, 141
341, 140, 48, 141, 151
342, 89, 122, 104, 154
343, 90, 96, 150, 119
344, 89, 154, 104, 64
345, 118, 122, 155, 130
346, 32, 73, 57, 152
347, 91, 95, 148, 127
348, 85, 155, 147, 95
349, 70, 151, 36, 77
350, 153, 140, 139, 146
351, 18, 48, 36, 141
352, 43, 54, 71, 147
353, 29, 147, 65, 54
354, 12, 136, 11, 32
355, 85, 129, 155, 95
356, 133, 29, 6, 43
357, 129, 155, 95, 127
358, 110, 150, 148, 120
359, 153, 40, 69, 64
360, 34, 40, 139, 14
361, 155, 85, 147, 88
362, 80, 113, 148, 126
363, 107, 117, 150, 105
364, 123, 155, 106, 109
365, 149, 92, 153, 154
366, 153, 144, 51, 143
367, 95, 54, 65, 147
368, 55, 88, 147, 71
369, 110, 150, 105, 148
370, 59, 78, 38, 149
371, 50, 149, 62, 59
372, 151, 37, 50, 62
373, 83, 106, 86, 148
374, 149, 96, 150, 90
375, 149, 66, 62, 90
376, 107, 84, 150, 117
377, 37, 151, 50, 142
378, 84, 149, 87, 150
379, 150, 113, 120, 119
380, 84, 107, 150, 87
381, 142, 19, 49, 141
382, 95, 155, 147, 148
383, 148, 114, 98, 91
384, 147, 33, 68, 152
*material, name=mat
*elastic
210e9, 0.3
*density
7850.0
*boundary
1, 1, 3
2, 1, 3
3, 1, 3
4, 1, 3
5, 1, 3
6, 1, 3
7, 1, 3
8, 1, 3
9, 1, 3
10, 1, 3
11, 1, 3
12, 1, 3
13, 1, 3
131, 1, 3
132, 1, 3
133, 1, 3
134, 1, 3
135, 1, 3
136, 1, 3
137, 1, 3
138, 1, 3
*solid section, material=mat, elset=elems
*step
*static
*cload
146, 3, -1e9
*dload
elems, GRAV, 9.81, 0, 0, -1
*end step

一番節点の変位量が大きい点を比較してみます。

一致することが確認できました。