Pythonによる２次元応力解析プログラム（四角形要素）の理論を少し詳しく解説する（改訂版）

• 2017.12.24 投稿
• 2018.08.05 修正（プログラム改良）

（情報）初期のプログラムから以下の点を変更しています．

• 疎行列処理導入による高速化
• 関数名を小文字にした
• データ入力部を関数にした
• 計算用関数（データ入力・プリントアウト部を除く）の引数を全てスカラーとした
• 無意味な配列初期化の部分をカットした
• 配列初期化時の形状をタプルで指定するようにした（何故かこれまでリストで渡していた）
• メインルーチンも関数化した．これによりグローバル変数を関数の引数で読み込むなどのおかしな点はなくなったと思います．

はじめに

２次元応力解析プログラムの四角形要素版の解説です．
基本となる剛性方程式，節点の変位拘束の処理方法などは２次元応力解析（三角形一次要素版）と同じです．

基本的な事項は，以下の記事にまとめてあります．

• 有限要素法とPythonプログラム（三角形一次要素を用いた２次元応力解析プログラム）

理論展開およびプログラムには，以下の条件が入っています．

• 線形弾性材料の微小変形を仮定している
• 物体力として慣性力を考慮する
• 初期ひずみとして温度変化によるひずみを考慮する

大規模連立方程式の処理を考慮し，疎行列処理を導入しています．

このシリーズの投稿のリンク

• Python：FEMでの疎行列計算利用による高速化（３次元骨組構造解析での事例）
• Pythonによる２次元応力解析プログラム（四角形要素）の理論を少し詳しく解説する（改訂版）
• Pythonによる平面骨組構造解析プログラムの理論を少し詳しく解説する（改訂版）
• Pythonによる２次元飽和ー不飽和定常浸透流解析プログラム（三角形要素）の理論を少し詳しく解説する（改訂版）
• Pythonによる２次元飽和ー不飽和定常浸透流解析プログラム（四角形要素）の理論を少し詳しく解説する（改訂版）

四角形要素を用いた２次元応力解析の理論

要素剛性方程式

変位を求めるための方程式

\begin{equation}
[\boldsymbol{k}]\{\boldsymbol{u_{nd}}\}=\{\boldsymbol{f}\}+\{\boldsymbol{f_t}\}+\{\boldsymbol{f_b}\}
\end{equation}

\begin{align}
&[\boldsymbol{k}]=\int_V[\boldsymbol{B}]^T[\boldsymbol{D_e}][\boldsymbol{B}]dV  & & \text{（要素剛性行列）}\\
&\{\boldsymbol{f}\}=\int_A [\boldsymbol{N}]^T\{\boldsymbol{S}\}dA  & & \text{（節点外力ベクトル）}\\
&\{\boldsymbol{f_t}\}=\int_V[\boldsymbol{B}]^T[\boldsymbol{D_e}]\{\boldsymbol{\epsilon_0}\}dV  & & \text{（節点初期歪ベクトル）}\\
&\{\boldsymbol{f_b}\}=\gamma\left(\int_V [\boldsymbol{N}]^T [\boldsymbol{N}]dV\right)\{\boldsymbol{w_{nd}}\} & & \text{（節点慣性力ベクトル）}
\end{align}

要素内部応力計算式

\begin{equation}
\{\boldsymbol{\sigma}\}=[\boldsymbol{D_e}]\{\boldsymbol{\epsilon}-\boldsymbol{\epsilon_0}\}
\end{equation}

\begin{align}
&\{\boldsymbol{\epsilon}\}  & &\text{：剛性方程式を解いて求めた変位から求めた要素のひずみベクトル}\\
&\{\boldsymbol{\epsilon_0}\} & &\text{：温度ひずみなどの初期ひずみ}\\
&[\boldsymbol{D_e}]          & &\text{：応力ーひずみ関係行列}
\end{align}

要素剛性方程式の基本形

直交座標系$(x-y)$平面上の四角形要素を考える．
２次元応力解析のための四角形要素は１要素４節点であり，平面問題のため，１節点２自由度（x方向変位およびy方向変位）を有する．
よって基本的な要素剛性方程式は，以下の形をとる．ここに，四角形要素を構成する節点番号を反時計回りに，$i$，$j$，$k$，$l$ とする．

\begin{equation}
\{\boldsymbol{f}\}=[\boldsymbol{k}]\{\boldsymbol{u_{nd}}\}
\end{equation}

\begin{equation}
\begin{Bmatrix} f_{x,i} \\ f_{y,i} \\ f_{x,j} \\ f_{y,j} \\ f_{x,k} \\ f_{y,k} \\ f_{x,l} \\ f_{y,l} \end{Bmatrix}
=\begin{bmatrix}
k_{11} & k_{12} & k_{13} & k_{14} & k_{15} & k_{16} & k_{17} & k_{18} \\
k_{21} & k_{22} & k_{23} & k_{24} & k_{25} & k_{26} & k_{27} & k_{28} \\
k_{31} & k_{32} & k_{33} & k_{34} & k_{35} & k_{36} & k_{37} & k_{38} \\
k_{41} & k_{42} & k_{43} & k_{44} & k_{45} & k_{46} & k_{47} & k_{48} \\
k_{51} & k_{52} & k_{53} & k_{54} & k_{55} & k_{56} & k_{57} & k_{58} \\
k_{61} & k_{62} & k_{63} & k_{64} & k_{65} & k_{66} & k_{67} & k_{68} \\
k_{71} & k_{72} & k_{73} & k_{74} & k_{75} & k_{76} & k_{77} & k_{78} \\
k_{81} & k_{82} & k_{83} & k_{84} & k_{85} & k_{86} & k_{87} & k_{88} 
\end{bmatrix}
\begin{Bmatrix} u_i \\ v_i \\ u_j \\ v_j \\ u_k \\ v_k \\ u_l \\ v_l \end{Bmatrix}
\end{equation}

\begin{align}
&f_{x,i} & &\text{：節点$i$の$x$方向荷重} & & u_i & &\text{：節点$i$の$x$方向変位}\\
&f_{y,i} & &\text{：節点$i$の$y$方向荷重} & & v_i & &\text{：節点$i$の$y$方向変位}\\
&f_{x,j} & &\text{：節点$j$の$x$方向荷重} & & u_i & &\text{：節点$j$の$x$方向変位}\\
&f_{y,j} & &\text{：節点$j$の$y$方向荷重} & & v_i & &\text{：節点$j$の$y$方向変位}\\
&f_{x,k} & &\text{：節点$k$の$x$方向荷重} & & u_i & &\text{：節点$k$の$x$方向変位}\\
&f_{y,k} & &\text{：節点$k$の$y$方向荷重} & & v_i & &\text{：節点$k$の$y$方向変位}\\
&f_{x,l} & &\text{：節点$l$の$x$方向荷重} & & u_l & &\text{：節点$l$の$x$方向変位}\\
&f_{y,l} & &\text{：節点$l$の$y$方向荷重} & & v_l & &\text{：節点$l$の$y$方向変位}
\end{align}

線形弾性論による基本式

平面応力状態での応力ーひずみ関係式

\begin{equation}
\begin{Bmatrix} \sigma_x \\ \sigma_y \\ \tau_{xy} \end{Bmatrix}
=\frac{E}{1-\nu^2}\begin{bmatrix} 1 & \nu & 0 \\ \nu & 1 & 0 \\ 0 & 0 & \cfrac{1-\nu}{2} \end{bmatrix}
\begin{Bmatrix} \epsilon_x-\alpha T \\ \epsilon_y-\alpha T \\ \gamma_{xy} \end{Bmatrix}
\end{equation}

平面ひずみ状態での応力ーひずみ関係式

\begin{equation}
\begin{Bmatrix} \sigma_x \\ \sigma_y \\ \tau_{xy} \end{Bmatrix}
=\frac{E}{(1+\nu)(1-2\nu)}\begin{bmatrix} 1-\nu & \nu & 0 \\ \nu & 1-\nu & 0 \\ 0 & 0 & \cfrac{1-2\nu}{2} \end{bmatrix}
\begin{Bmatrix} \epsilon_x-(1+\nu)\alpha T \\ \epsilon_y-(1+\nu)\alpha T \\ \gamma_{xy} \end{Bmatrix}
\end{equation}

ひずみー変位関係式

\begin{equation}
\{\boldsymbol{\epsilon}\}
=\begin{Bmatrix} \epsilon_x \\ \epsilon_y \\ \gamma_{xy} \end{Bmatrix}
=\begin{Bmatrix} \cfrac{\partial u}{\partial x} \\ \cfrac{\partial v}{\partial y} \\ \cfrac{\partial u}{\partial y}+\cfrac{\partial v}{\partial x} \end{Bmatrix}
\end{equation}

記号の定義

\begin{align}
&\sigma_x     & &\text{：x方向直応力}    & & \epsilon_x  & &\text{：x方向直ひずみ}\\
&\sigma_y     & &\text{：y方向直応力}    & & \epsilon_y  & &\text{：y方向直ひずみ}\\
&\tau_{xy}    & &\text{：せん断応力}     & & \gamma_{xy} & &\text{：せん断ひずみ}\\
&E            & &\text{：要素弾性係数}   & & \nu         & &\text{：要素ポアソン比}\\
&\alpha       & &\text{：要素線膨張係数} & & T           & &\text{：要素温度上昇量}\\
&u            & &\text{：x方向変位}      & & v           & &\text{：y方向変位}
\end{align}

以下に，プログラムを作成するために必要な，具体的な要素を求めていく．

ひずみ-変位関係行列

四角形要素内任意点の変位 $u$, $v$ を，内挿関数 $N$ を用いて以下のようにおく．
$u_{i,j,k,l}$ および $v_{i,j,k,l}$ は要素を構成する節点の変位を示す．

\begin{align}
u=&N_1\cdot u_i+N_2\cdot u_j+N_3\cdot u_k+N_4\cdot u_l \\
v=&N_1\cdot v_i+N_2\cdot v_j+N_3\cdot v_k+N_4\cdot v_l
\end{align}

要素内任意点の変位を${\boldsymbol{u}}$，節点変位を${\boldsymbol{u_{nd}}}$とすれば，

\begin{equation}
\{\boldsymbol{u}\}=\begin{Bmatrix} u \\ v \end{Bmatrix}
=\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}
\begin{Bmatrix}
u_i \\ v_i \\ u_j \\ v_j \\ u_k \\ v_k \\ u_l \\ v_l
\end{Bmatrix}
=[\boldsymbol{N}]\{\boldsymbol{u_{nd}}\}
\end{equation}

要素内任意点のひずみ${\boldsymbol{\epsilon}}$は，節点変位${\boldsymbol{u_{nd}}}$を用いて，

\begin{equation}
\{\boldsymbol{\epsilon}\}=\begin{Bmatrix} \cfrac{\partial u}{\partial x} \\ \cfrac{\partial v}{\partial y} \\ \cfrac{\partial u}{\partial y}+\cfrac{\partial v}{\partial x} \end{Bmatrix}
=\begin{bmatrix}
\cfrac{\partial N_1}{\partial x} & 0 & \cfrac{\partial N_2}{\partial x} & 0 & \cfrac{\partial N_3}{\partial x} & 0 & \cfrac{\partial N_4}{\partial x} & 0 \\
0 & \cfrac{\partial N_1}{\partial y} & 0 & \cfrac{\partial N_2}{\partial y} & 0 & \cfrac{\partial N_3}{\partial y} & 0 & \cfrac{\partial N_4}{\partial y} \\
\cfrac{\partial N_1}{\partial y} & \cfrac{\partial N_1}{\partial x} & \cfrac{\partial N_2}{\partial y} & \cfrac{\partial N_2}{\partial x} & \cfrac{\partial N_3}{\partial y} & \cfrac{\partial N_3}{\partial x} & \cfrac{\partial N_4}{\partial y} & \cfrac{\partial N_4}{\partial x}
\end{bmatrix}
\begin{Bmatrix}
u_i \\ v_i \\ u_j \\ v_j \\ u_k \\ v_k \\ u_l \\ v_l
\end{Bmatrix}
=[\boldsymbol{B}]\{\boldsymbol{u_{nd}}\}
\end{equation}

内挿関数

実際の計算においては，内挿関数の数式表現とその微分が必要となるため，これらを求めておくことにする．

以下，$(a,b)$ は，$a=[−1,1]$，$b=[−1,1]$ の範囲を有する正規化パラメータ座標である．

内挿関数の数式表現

２次元問題の四角形要素は，１要素４節点，１節点２自由度であるため，１要素あたり８自由度を有する．そこで，要素内任意点の変位を，４節点の正規化座標 $(a,b)$ および８個の未定係数 $\alpha_{1\sim 8}$ を用いて以下のように仮定する．

\begin{equation}
\begin{Bmatrix} u \\ v \end{Bmatrix}
=\begin{bmatrix}
1 & a & b & ab & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 1 & a & b & ab \\
\end{bmatrix}
\begin{Bmatrix} \alpha_1 \\ \alpha_2 \\ \alpha_3 \\ \alpha_4 \\ \alpha_5 \\ \alpha_6 \\ \alpha_7 \\ \alpha_8 \end{Bmatrix}
=[\boldsymbol{A}]\{\boldsymbol{\alpha}\}
\end{equation}

正規化座標を用いた節点座標を反時計回りに以下のようにおく．

\begin{align}
\text{節点}          & & \text{節点 $i$} & & \text{節点 $j$} & & \text{節点 $k$} & & \text{節点 $l$} \\
\text{座標値$(a,b)$} & & (-1,-1) & & (1,-1) & & (1,1) & & (-1,1)
\end{align}

各節点の変位を $\{\alpha\}$ を用いて表すと，

\begin{equation}
\begin{Bmatrix} u_i \\ v_i \\ u_j \\ v_j \\ u_k \\ v_k \\ u_l \\ v_l \end{Bmatrix}
=\begin{bmatrix}
1 & -1 & -1 &  1 &  0 &  0 &  0 &  0 \\
0 &  0 &  0 &  0 &  1 & -1 & -1 &  1 \\
1 &  1 & -1 & -1 &  0 &  0 &  0 &  0 \\
0 &  0 &  0 &  0 &  1 &  1 & -1 & -1 \\
1 &  1 &  1 &  1 &  0 &  0 &  0 &  0 \\
0 &  0 &  0 &  0 &  1 &  1 &  1 &  1 \\
1 & -1 &  1 & -1 &  0 &  0 &  0 &  0 \\
0 &  0 &  0 &  0 &  1 & -1 &  1 & -1
\end{bmatrix}
\begin{Bmatrix} \alpha_1 \\ \alpha_2 \\ \alpha_3 \\ \alpha_4 \\ \alpha_5 \\ \alpha_6 \\ \alpha_7 \\ \alpha_8 \end{Bmatrix}
\end{equation}

これを $\{\alpha\}$ について解けば，

\begin{equation}
\begin{Bmatrix} \alpha_1 \\ \alpha_2 \\ \alpha_3 \\ \alpha_4 \\ \alpha_5 \\ \alpha_6 \\ \alpha_7 \\ \alpha_8 \end{Bmatrix}
=\left[\begin{matrix}
\frac{1}{4} & 0 & \frac{1}{4} & 0 & \frac{1}{4} & 0 & \frac{1}{4} & 0\\
- \frac{1}{4} & 0 & \frac{1}{4} & 0 & \frac{1}{4} & 0 & - \frac{1}{4} & 0\\
- \frac{1}{4} & 0 & - \frac{1}{4} & 0 & \frac{1}{4} & 0 & \frac{1}{4} & 0\\
\frac{1}{4} & 0 & - \frac{1}{4} & 0 & \frac{1}{4} & 0 & - \frac{1}{4} & 0\\
0 & \frac{1}{4} & 0 & \frac{1}{4} & 0 & \frac{1}{4} & 0 & \frac{1}{4}\\
0 & - \frac{1}{4} & 0 & \frac{1}{4} & 0 & \frac{1}{4} & 0 & - \frac{1}{4}\\
0 & - \frac{1}{4} & 0 & - \frac{1}{4} & 0 & \frac{1}{4} & 0 & \frac{1}{4}\\
0 & \frac{1}{4} & 0 & - \frac{1}{4} & 0 & \frac{1}{4} & 0 & - \frac{1}{4}
\end{matrix}\right]
\begin{Bmatrix} u_i \\ v_i \\ u_j \\ v_j \\ u_k \\ v_k \\ u_l \\ v_l \end{Bmatrix}
=[\boldsymbol{C}]\{\boldsymbol{u_{nd}}\}
\end{equation}

よって，内挿関数 $[\boldsymbol{N}]$ は，

\begin{equation}
[\boldsymbol{N}]=[\boldsymbol{A}][\boldsymbol{C}]
=\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}
\end{equation}

\begin{align}
N_1=&\frac{1}{4}(1-a)(1-b) \qquad N_2=\frac{1}{4}(1+a)(1-b) \\
N_3=&\frac{1}{4}(1+a)(1+b) \qquad N_4=\frac{1}{4}(1-a)(1+b)
\end{align}

内挿関数の微分

\begin{align}
&\cfrac{\partial N_1}{\partial a}=-\cfrac{1}{4}(1-b)  &\cfrac{\partial N_2}{\partial a}=+\cfrac{1}{4}(1-b) & &\cfrac{\partial N_3}{\partial a}=+\cfrac{1}{4}(1+b) & &\cfrac{\partial N_4}{\partial a}=-\cfrac{1}{4}(1+b) \\
&\cfrac{\partial N_1}{\partial b}=-\cfrac{1}{4}(1-a)  &\cfrac{\partial N_2}{\partial b}=-\cfrac{1}{4}(1+a) & &\cfrac{\partial N_3}{\partial b}=+\cfrac{1}{4}(1+a) & &\cfrac{\partial N_4}{\partial b}=+\cfrac{1}{4}(1-a)
\end{align}

ここで，内挿関数 $N_i$ に関する微分は，Jacobi行列 $[J]$ を用いて以下の通り表せる．

\begin{equation}
\begin{Bmatrix} \cfrac{\partial N_i}{\partial a} \\ \cfrac{\partial N_i}{\partial b} \end{Bmatrix}
=[J] \begin{Bmatrix} \cfrac{\partial N_i}{\partial x} \\ \cfrac{\partial N_i}{\partial y} \end{Bmatrix} \qquad
\begin{Bmatrix} \cfrac{\partial N_i}{\partial x} \\ \cfrac{\partial N_i}{\partial y} \end{Bmatrix}
=[J]^{-1} \begin{Bmatrix} \cfrac{\partial N_i}{\partial a} \\ \cfrac{\partial N_i}{\partial b} \end{Bmatrix}
\end{equation}

\begin{equation}
[J]=\begin{bmatrix}
\cfrac{\partial x}{\partial a} & \cfrac{\partial y}{\partial a} \\
\cfrac{\partial x}{\partial b} & \cfrac{\partial y}{\partial b}
\end{bmatrix}
=\begin{bmatrix}
\displaystyle \sum_{i=1}^4\left(\cfrac{\partial N_i}{\partial a}x_i\right) & \displaystyle \sum_{i=1}^4\left(\cfrac{\partial N_i}{\partial a}y_i\right) \\
\displaystyle \sum_{i=1}^4\left(\cfrac{\partial N_i}{\partial b}x_i\right) & \displaystyle \sum_{i=1}^4\left(\cfrac{\partial N_i}{\partial b}y_i\right)
\end{bmatrix}
=\begin{bmatrix} J_{11} & J_{12} \\\ J_{21} & J_{22} \end{bmatrix}
\end{equation}

\begin{equation}
[J]^{-1}=\frac{1}{\det(J)}\begin{bmatrix} J_{22} & -J_{12} \\ -J_{21} & J_{11} \end{bmatrix} \qquad
\det(J)=J_{11}\cdot J_{22}-J_{12}\cdot J_{21}
\end{equation}

上記より，

\begin{equation}
\cfrac{\partial N_i}{\partial x}=\cfrac{1}{\det(J)}\left\{J_{22}\cfrac{\partial N_i}{\partial a}-J_{12}\cfrac{\partial N_i}{\partial b}\right\} \qquad
\cfrac{\partial N_i}{\partial y}=\cfrac{1}{\det(J)}\left\{-J_{21}\cfrac{\partial N_i}{\partial a}+J_{11}\cfrac{\partial N_i}{\partial b}\right\}
\end{equation}

以下に $[J]$ の要素を示す．ここに $x_{i,j,k,l}$ および $y_{i,j,k,l}$ は要素を構成する節点座標である．

\begin{align}
&J_{11}=\frac{\partial N_1}{\partial a} x_i+\frac{\partial N_2}{\partial a} x_j+\frac{\partial N_3}{\partial a} x_k+\frac{\partial N_4}{\partial a} x_l \\
&J_{12}=\frac{\partial N_1}{\partial a} y_i+\frac{\partial N_2}{\partial a} y_j+\frac{\partial N_3}{\partial a} y_k+\frac{\partial N_4}{\partial a} y_l \\
&J_{21}=\frac{\partial N_1}{\partial b} x_i+\frac{\partial N_2}{\partial b} x_j+\frac{\partial N_3}{\partial b} x_k+\frac{\partial N_4}{\partial b} x_l \\
&J_{22}=\frac{\partial N_1}{\partial b} y_i+\frac{\partial N_2}{\partial b} y_j+\frac{\partial N_3}{\partial b} y_k+\frac{\partial N_4}{\partial b} y_l
\end{align}

以上により，$[J]$ の各要素が定まり，$[\boldsymbol{N}]$ の微分値も定めることができる．

なお，この段階では $[J]$ および $[\boldsymbol{N}]$ とその微分値は，変数 $a$ および $b$ の関数として定義されていることに注意が必要である．

応力-ひずみ関係式

応力ーひずみ関係マトリックスは，以下の通り．

\begin{align}
&\text{平面応力：}& &[\boldsymbol{D_e}]=\frac{E}{1-\nu^2}\begin{bmatrix} 1 & \nu & 0 \\ \nu & 1 & 0 \\ 0 & 0 & \cfrac{1-\nu}{2} \end{bmatrix} \qquad
\begin{matrix} E & \text{：弾性係数} \\ \nu & \text{：ポアソン比} \end{matrix} \\
&\text{平面ひずみ：}& &[\boldsymbol{D_e}]=\frac{E}{(1+\nu)(1-2\nu)}\begin{bmatrix} 1-\nu & \nu & 0 \\ \nu & 1-\nu & 0 \\ 0 & 0 & \cfrac{1-2\nu}{2} \end{bmatrix}
\end{align}

ガウス積分

以降の要素剛性マトリックスや荷重ベクトル算定で必要となる要素に対する面積積分は，以下に示すガウス積分により行う．

積分点を 4 点 ($n=2$) とすれば，$a$，$b$ および重み $H$ の値は下表の通りであり，$a$，$b$ を変化させた 4 回の値の総和が積分の近似値となる．また，重みが $H=1$ であることから，計算は更に単純になる．

\begin{equation}
\int_{-1}^{1}\int_{-1}^{1} f(a,b)\cdot da\cdot db
\doteqdot \sum_{i=1}^n\sum_{j=1}^n H_i\cdot H_j\cdot f(a_i,b_j)
= \sum_{kk=1}^4 f(a_{kk}, b_{kk})
\end{equation}

\begin{align}
& i  & j & &  a              & &  b              & & H              & & kk \\
& 1  & 1 & & -0.57735 02692  & & -0.57735 02692  & & 1.00000 00000  & & 1  \\\
& 1  & 2 & & +0.57735 02692  & & -0.57735 02692  & & 1.00000 00000  & & 2  \\\
& 2  & 1 & & +0.57735 02692  & & +0.57735 02692  & & 1.00000 00000  & & 3  \\\
& 2  & 2 & & -0.57735 02692  & & +0.57735 02692  & & 1.00000 00000  & & 4
\end{align}

要素剛性マトリックス

4節点パラメトリック要素の剛性行列$[\boldsymbol{k}]$は，要素の板厚を一定値$t$とし，応力-ひずみ関係行列を$[\boldsymbol{D_e}]$とすれば以下の通りとなる．

\begin{align}
[\boldsymbol{k}]=&\int_V[\boldsymbol{B}]^T[\boldsymbol{D_e}][\boldsymbol{B}]dV \\
=&t\int_{-1}^{1}\int_{-1}^{1}[\boldsymbol{B}]^{T}[\boldsymbol{D_e}][\boldsymbol{B}]\cdot\det(J)\cdot da\cdot db \\
=&t\cdot \sum_{kk=1}^4 \left\{[\boldsymbol{B}]^{T}[\boldsymbol{D_e}][\boldsymbol{B}]\cdot\det(J)\right\}_{kk}
\end{align}

節点温度荷重ベクトル

初期ひずみとして温度ひずみを考える．
温度変化によるひずみ${\boldsymbol{\epsilon_0}}$による節点力ベクトル${\boldsymbol{f_t}}$は以下の通りとなる．

\begin{align}
\{\boldsymbol{f_t}\}=&\int_V[\boldsymbol{B}]^T[\boldsymbol{D_e}]\{\boldsymbol{\epsilon_0}\}dV \\
=&t\int_{-1}^{1}\int_{-1}^{1}[\boldsymbol{B}]^{T}[\boldsymbol{D_e}]\{\boldsymbol{\epsilon_0}\}\cdot\det(J)\cdot da\cdot db \\
=&t\cdot \sum_{kk=1}^4 \left\{[\boldsymbol{B}]^{T}[\boldsymbol{D_e}]\{\boldsymbol{\epsilon_0}\}\cdot\det(J)\right\}_{kk}
\end{align}

温度変化によるひずみの成分は以下のとおりである．

\begin{equation}
T_p=N_1\cdot\phi_i+N_2\cdot\phi_j+N_3\cdot\phi_k+N_4\cdot\phi_l
\end{equation}

\begin{align}
\text{平面応力：}&\{\boldsymbol{\epsilon_0}\}=\begin{Bmatrix} \alpha\cdot T_p \\ \alpha\cdot T_p \\ 0 \end{Bmatrix}
\\
\text{平面ひずみ：}&\{\boldsymbol{\epsilon_0}\}=\begin{Bmatrix} (1+\nu)\cdot \alpha\cdot T_p \\ (1+\nu)\cdot \alpha\cdot T_p \\ 0 \end{Bmatrix}
\end{align}

ここに，$\alpha$は線膨張係数，$\nu$はポアソン比，$T_p$は要素任意点の温度変化量，${\phi_i~~\phi_j~~\phi_k~~\phi_l}^T$は節点温度変化量，$[N_1N_2N_3~~N_4]$は$[\boldsymbol{N}]$の要素である．また，温度変化量の符号は，温度上昇を正とする．

節点慣性力ベクトル

\begin{align}
\{\boldsymbol{f_b}\}=&\gamma \cdot \left(\int_V [\boldsymbol{N}]^T [\boldsymbol{N}]dV\right)\{\boldsymbol{w_{nd}}\} \\
=&\gamma\cdot t\cdot\int_{A}[\boldsymbol{N}]^T[\boldsymbol{N}]dA\cdot\{\boldsymbol{w_{nd}}\} \\
=&\gamma\cdot t \cdot\sum_{kk=1}^4 \left\{[\boldsymbol{N}]^T[\boldsymbol{N}]\cdot\det(J)\right\}_ {kk}\cdot\{\boldsymbol{w_{nd}}\}
\end{align}

\begin{equation}
\{\boldsymbol{w_{nd}}\}=\begin{Bmatrix}k_h & k_v & k_h & k_v & k_h & k_v & k_h & k_v \end{Bmatrix}^T
\end{equation}

ここに$\gamma$は要素の単位体積重量（質量$\times$重力加速度），$k_h$および$k_v$は水平方向・鉛直方向加速度（重力加速度に対する比率）である．

節点外力ベクトル

節点外力ベクトルについては，ここで紹介するプログラムでは，直接，等価節点外力を入力データファイルから入力することとしており，特別な計算処理は行っていない．

要素応力ベクトル

要素応力ベクトル ${\boldsymbol{\sigma}}$ は，以下のように計算できる．

\begin{align}
&\{\boldsymbol{\epsilon}\}=[\boldsymbol{B}]\{\boldsymbol{u_{nd}}\} & &\text{（節点変位から求まるひずみ）}\\
&T_p=N_1\cdot\phi_i+N_2\cdot\phi_j+N_3\cdot\phi_k+N_4\cdot\phi_l & &\text{（要素内任意点の温度変化量）}\\
&\{\boldsymbol{\epsilon_0}\}=\begin{Bmatrix} \alpha\cdot T_p \\ \alpha\cdot T_p \\ 0 \end{Bmatrix} & &\text{（平面応力）}\\
&\{\boldsymbol{\epsilon_0}\}=\begin{Bmatrix} (1+\nu)\cdot \alpha\cdot T_p \\ (1+\nu)\cdot \alpha\cdot T_p \\ 0 \end{Bmatrix} & &\text{（平面ひずみ）}\\
&\{\boldsymbol{\sigma}\}=[\boldsymbol{D_e}]\{\boldsymbol{\epsilon}-\boldsymbol{\epsilon_0}\} & &\text{（要素応力ベクトル）}
\end{align}

ここに，$N_1 \sim N_4$ は内挿関数，$\phi_i$，$\phi_j$，$\phi_k$，$\phi_l$ は要素節点 $i$，$j$，$k$，$l$ の温度変化量である．

四角形要素を用いた２次元応力解析プログラム

プログラムの概要

• プログラミング言語はPythonとする
• 等方弾性体の２次元応力解析を行う
• 平面ひずみおよび平面応力問題を扱える
• 要素は，１要素４節点４ガウス積分点を有する四角形要素とする
• 座標システムは，平面上の直行座標系右方向を x 軸正方向，上方向を y 軸正方向とする．
• 荷重として，節点外力，節点強制変位（ゼロ変位含む），節点温度変化，要素への慣性力を扱える
• 連立一次方程式は，その大きさ（１節点自由度ｘ総節点数）を変えずに，節点強制変位導入による処理を行い解いている．
• 連立一次方程式は，大きな自由度での計算効率を考慮し，疎行列計算ライブラリの scipy.sparse.linalg.spsolve(A,b) で解く
• 出力は，節点変位および要素平均応力とする
• 入力データファイル，出力データファイルの文字区切りは空白とする
• 入力データファイル名，出力データファイル名は，コマンドライン引数として入力する

扱える境界条件の説明

項目	説明
節点外力	載荷節点数，積荷節点番号，荷重値を指定
節点温度変化	全節点の温度変化量を指定．温度上昇を正値とする
要素慣性力	材料特性データの中で水平方向および鉛直方向の加速度を重力加速度に対する比で指定
節点強制変位	強制変位を与える節点数，節点番号，変位値を指定

プログラムの構成（プログラム名：py_fem_pl4x.py）

1. モジュール読み込み
2. 関数：データ入力 (def inpdata_pl4)
   (1) 基本数値読み込み
   (2) 配列宣言
   (3) 材料特性データ読み込み
   (4) 要素構成節点番号および材料番号読み込み
   (5) 節点座標および節点温度変化量読み込み
   (6) 境界条件読み込み
   (7) 節点荷重データ読み込み
3. 関数：入力データ書き出し (def prinp_pl4)
4. 関数：計算結果書き出し (def prout_pl4)
5. 関数：Dマトリックス作成 (def dmat_pl4)
6. 関数：Bマトリックスおよびヤコビアン作成 (def bmat_pl4)
7. 関数：ガウス積分点指定 (defab_pl4)
8. 関数：要素剛性マトリックス作成 (def sm_pl4)
9. 関数：要素応力計算 (def calst_pl4)
10. 関数：主応力計算 (def pst_pl4)
11. 関数：節点温度荷重ベクトル作成 (def tfvec_pl4)
12. 関数：節点慣性力ベクトル作成 (def bfvec_pl4)
13. 関数：メインルーチン
    (1) コマンドラインより入出力ファイル名取得
    (2) データ入力
    (3) 入力データ書き出し
    (4) 剛性マトリックス組み立て
    (5) 境界条件処理（変位既知自由度の処理）
    (6) 剛性方程式を解く（変位計算）
    (7) 境界条件処理（既知変位の組み込み）
    (8) 要素応力計算
    (9) 計算結果書き出し
    (10) 計算時間等情報表示
14. メインルーチン実行

実行コマンドと入出力データ

解析実行コマンド

python3 py_fem_pl4x.py inp.txt out.txt

変数	説明
py_fem_pl4.py	Python による FEM プログラム名
inp.txt	入力データファイル名（空白区切りデータ）
out.txt	出力データファイル名（空白区切りデータ）

入力データファイルフォーマット

npoin  nele  nsec  npfix  nlod  NSTR  # Basic values for analysis
t  E  po  alpha  gamma  gkh  gkv      # Material properties
    ..... (1 to nsec) .....           #
node1  node2  node3  node4  isec      # Element connectivity, material set number
    ..... (1 to nele) .....           #   (counter-clockwise order of node number)
x  y  deltaT                          # Coordinate, Temperature change of node
    ..... (1 to npoin) .....          #
node  kox  koy  rdisx  rdisy          # Restricted node and restrict conditions
    ..... (1 to npfix) .....          #
node  fx  fy                          # Loaded node and loading conditions
    ..... (1 to nlod) .....           #   (omit data input if nlod=0)

変数	説明
npoin, nele, nsec	節点数，要素数，材料特性数
npfix, nlod	拘束節点数，載荷節点数
NSTR	応力状態（平面歪: 0，平面応力: 1）
t, E, po, alpha	板厚，弾性係数，ポアソン比，線膨張係数
gamma, gkh, gkv	単位体積重量，水平及び鉛直方向加速度（gの比）
node1, node2, node3, node4, isec	節点1，節点2，節点3，節点4，材料特性番号
x, y, deltaT	節点x座標，節点y座標，節点温度変化
node, kox, koy	拘束節点番号，x及びy方向拘束の有無（拘束: 1，自由: 0）
rdisx, rdisy	x及びy方向変位（無拘束でも0を入力）
node, fx, fy	載荷重節点番号，x方向荷重，y方向荷重

出力データファイルフォーマット

npoin  nele  nsec npfix  nlod   NSTR
    (Each value of above)
sec  t  E  po  alpha  gamma  gkh  gkv
    sec   : Material number
    t     : Element thickness
    E     : Elastic modulus
    po    : Poisson's ratio
    alpha : Thermal expansion coefficient
    gamma : Unit weight
    gkh   : Acceleration in x-direction (ratio to g)
    gkv   : Acceleration in y-direction (ratio to g)
    ..... (1 to nsec) .....
node  x  y  fx  fy  deltaT  kox  koy
    node   : Node number
    x      : x-coordinate
    y      : y-coordinate
    fx     : Load in x-direction
    fy     : Load in y-direction
    deltaT : Temperature change of node
    kox    : Index of restriction in x-direction (0: free, 1: fixed)
    koy    : Index of restriction in x-direction (0: free, 1: fixed)
    ..... (1 to npoin) .....
node  kox  koy  rdis_x  rdis_y
    node    : Node number
    kox     : Index of restriction in x-direction (0: free, 1: fixed)
    koy     : Index of restriction in y-direction (0: free, 1: fixed)
    rdis_x  : Given displacement in x-direction
    tdis_y  : Given displacement in y-direction
    ..... (1 to npfix) .....
elem  i  j  k  l  sec
    elem : Element number
    i    : Node number of start point
    j    : Node number of second point
    k    : Node number of third point
    l    : Node number of end point
    sec  : Material number
    ..... (1 to nele) .....
node  dis-x  dis-y
    node  : Node number
    dis-x : Displacement in x-direction
    dis-y : Displacement in y-direction
    ..... (1 to npoin) .....
elem  sig_x  sig_y  tau_xy  p1  p2  ang
    elem   : Element number
    sig_x  : Normal stress in x-direction
    sig_y  : Normal stress in y-direction
    tau_xy : Shear stress in x-y plane
    p1     : First principal stress
    p2     : Second principal stress
    ang    : Angle of the first principal stress
    ..... (1 to nele) .....
n=(total degrees of freedom)  time=(calculation time)

変数	説明
dis-x, dis-y	x方向変位，y方向変位
sig_x, sig_y, tau_xy	x方向直応力，y方向直応力，せん断応力
p1, p2, ang	第一主応力，第二主応力，第一主応力の方向

プログラミング

1. モジュール読み込み

• numpy：配列演算用
• scipy.sparse.linalg.spsolve：疎行列連立一次方程式解法
• scipy.sparse import csr_matrix：疎行列圧縮格納
• sys：コマンドライン引数使用
• time：計算時間計測用

#===================================
# 2D FEM Stress Analysis (4-nodes)
#===================================
import numpy as np
from scipy.sparse.linalg import spsolve
from scipy.sparse import csr_matrix
import sys
import time

2. 関数：データ入力 (def inpdata_pl4)

(1) 基本数値読み込み（npoin. nele, nsec, npfix, nlod, nstr）
(2) 配列宣言（ae, node, x, deltaT, mpfix, rdis, fp）
(3) 材料特性データ読み込み（配列 ae）
(4) 要素構成節点番号および材料番号読み込み（配列 node）
(5) 節点座標および節点温度変化量読み込み（配列 x, deltaT）
(6) 境界条件読み込み（配列 mpfix, rdis）
(7) 節点荷重データ読み込み（配列 fp）

def inpdata_pl4(fnameR,nod,nfree):
    f=open(fnameR,'r')
    text=f.readline()
    text=text.strip()
    text=text.split()
    npoin =int(text[0]) # Number of nodes
    nele  =int(text[1]) # Number of elements
    nsec  =int(text[2]) # Number of sections
    npfix =int(text[3]) # Number of restricted nodes
    nlod  =int(text[4]) # Number of loaded nodes
    nstr  =int(text[5]) # 0: plane strain, 1: plane stress
    # array declaration
    ae    =np.zeros((7,nsec),dtype=np.float64)      # Section characteristics
    node  =np.zeros((nod+1,nele),dtype=np.int)      # node-element relationship
    x     =np.zeros((2,npoin),dtype=np.float64)     # Coordinates of nodes
    deltaT=np.zeros(npoin,dtype=np.float64)         # Temperature change of node
    mpfix =np.zeros((nfree,npoin),dtype=np.int)     # Ristrict conditions
    rdis  =np.zeros((nfree,npoin),dtype=np.float64) # Ristricted displacement
    fp    =np.zeros(nfree*npoin,dtype=np.float64)   # External force vector
    # section characteristics
    for i in range(0,nsec):
        text=f.readline()
        text=text.strip()
        text=text.split()
        ae[0,i]=float(text[0]) #t    : :Plate thickness
        ae[1,i]=float(text[1]) #E    : Elastic modulus
        ae[2,i]=float(text[2]) #po   : Poisson's ratio
        ae[3,i]=float(text[3]) #alpha: Thermal expansion coefficient
        ae[4,i]=float(text[4]) #gamma: Unit weight of material
        ae[5,i]=float(text[5]) #gkh  : Acceleration in x-direction
        ae[6,i]=float(text[6]) #gkv  : Acceleration in y-direction
    # element-node
    for i in range(0,nele):
        text=f.readline()
        text=text.strip()
        text=text.split()
        node[0,i]=int(text[0]) #node_1
        node[1,i]=int(text[1]) #node_2
        node[2,i]=int(text[2]) #node_3
        node[3,i]=int(text[3]) #node_4
        node[4,i]=int(text[4]) #section characteristic number
    # node coordinates
    for i in range(0,npoin):
        text=f.readline()
        text=text.strip()
        text=text.split()
        x[0,i]   =float(text[0]) # x-coordinate
        x[1,i]   =float(text[1]) # y-coordinate
        deltaT[i]=float(text[2]) # Temperature change of node
    # boundary conditions (0:free, 1:restricted)
    for i in range(0,npfix):
        text=f.readline()
        text=text.strip()
        text=text.split()
        lp=int(text[0])             #fixed node
        mpfix[0,lp-1]=int(text[1])  #fixed in x-direction
        mpfix[1,lp-1]=int(text[2])  #fixed in y-direction
        rdis[0,lp-1]=float(text[3]) #fixed displacement in x-direction
        rdis[1,lp-1]=float(text[4]) #fixed displacement in y-direction
    # load
    if 0<nlod:
        for i in range(0,nlod):
            text=f.readline()
            text=text.strip()
            text=text.split()
            lp=int(text[0])           #loaded node
            fp[2*lp-2]=float(text[1]) #load in x-direction
            fp[2*lp-1]=float(text[2]) #load in y-direction
    f.close()
    return npoin,nele,nsec,npfix,nlod,nstr,ae,node,x,deltaT,mpfix,rdis,fp

3. 関数：入力データ書き出し (def prinp_pl4)

def prinp_pl4(fnameW,npoin,nele,nsec,npfix,nlod,nstr,ae,node,x,deltaT,mpfix,rdis,fp):
    fout=open(fnameW,'w')
    # print out of input data
    print('{0:>5s} {1:>5s} {2:>5s} {3:>5s} {4:>5s} {5:>6s}'
    .format('npoin','nele','nsec','npfix','nlod','nstr'),file=fout)
    print('{0:5d} {1:5d} {2:5d} {3:5d} {4:5d} {5:6d}'
    .format(npoin,nele,nsec,npfix,nlod,nstr),file=fout)
    print('{0:>5s} {1:>15s} {2:>15s} {3:>15s} {4:>15s} {5:>15s} {6:>10s} {7:>10s}'
    .format('sec','t','E','po','alpha','gamma','gkh','gkv'),file=fout)
    for i in range(0,nsec):
        print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e} {4:15.7e} {5:15.7e} {6:10.3f} {7:10.3f}'
        .format(i+1,ae[0,i],ae[1,i],ae[2,i],ae[3,i],ae[4,i],ae[5,i],ae[6,i]),file=fout)
    print('{0:>5s} {1:>15s} {2:>15s} {3:>15s} {4:>15s} {5:>15s} {6:>5s} {7:>5s}'
    .format('node','x','y','fx','fy','deltaT','kox','koy'),file=fout)
    for i in range(0,npoin):
        lp=i+1
        print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e} {4:15.7e} {5:15.7e} {6:5d} {7:5d}'
        .format(lp,x[0,i],x[1,i],fp[2*i],fp[2*i+1],deltaT[i],mpfix[0,i],mpfix[1,i]),file=fout)
    print('{0:>5s} {1:>5s} {2:>5s} {3:>15s} {4:>15s}'
    .format('node','kox','koy','rdis_x','rdis_y'),file=fout)
    for i in range(0,npoin):
        if 0<mpfix[0,i]+mpfix[1,i]:
            lp=i+1
            print('{0:5d} {1:5d} {2:5d} {3:15.7e} {4:15.7e}'
            .format(lp,mpfix[0,i],mpfix[1,i],rdis[0,i],rdis[1,i]),file=fout)
    print('{0:>5s} {1:>5s} {2:>5s} {3:>5s} {4:>5s} {5:>5s}'
    .format('elem','i','j','k','l','sec'),file=fout)
    for ne in range(0,nele):
        print('{0:5d} {1:5d} {2:5d} {3:5d} {4:5d} {5:5d}'
        .format(ne+1,node[0,ne],node[1,ne],node[2,ne],node[3,ne],node[4,ne]),file=fout)
    fout.close()

4. 関数：計算結果書き出し (def prout_pl4)

def prout_pl4(fnameW,npoin,nele,disg,fp,strs):
    fout=open(fnameW,'a')
    # displacement
    print('{0:>5s} {1:>15s} {2:>15s} {3:>15s} {4:>15s}'.format('node','dis-x','dis-y','fr-x','fr-y'),file=fout)
    for i in range(0,npoin):
        lp=i+1
        print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e} {4:15.7e}'.format(lp,disg[2*lp-2],disg[2*lp-1],fp[2*lp-2],fp[2*lp-1]),file=fout)
    # section force
    print('{0:>5s} {1:>15s} {2:>15s} {3:>15s} {4:>15s} {5:>15s} {6:>15s}'
    .format('elem','sig_x','sig_y','tau_xy','p1','p2','ang'),file=fout)
    for ne in range(0,nele):
        sigx =strs[0,ne]
        sigy =strs[1,ne]
        tauxy=strs[2,ne]
        ps1,ps2,ang=pst_pl4(sigx,sigy,tauxy)
        print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e} {4:15.7e} {5:15.7e} {6:15.7e}'
        .format(ne+1,sigx,sigy,tauxy,ps1,ps2,ang),file=fout)
    fout.close()

5. 関数：Dマトリックス作成 (def dmat_pl4)

\begin{align}
&\text{平面ひずみ（nstr=0）：}& &[\boldsymbol{D_e}]=\frac{E}{(1+\nu)(1-2\nu)}\begin{bmatrix} 1-\nu & \nu & 0 \\ \nu & 1-\nu & 0 \\ 0 & 0 & \cfrac{1-2\nu}{2} \end{bmatrix}\\
&\text{平面応力（nstr=1）：}& &[\boldsymbol{D_e}]=\frac{E}{1-\nu^2}\begin{bmatrix} 1 & \nu & 0 \\ \nu & 1 & 0 \\ 0 & 0 & \cfrac{1-\nu}{2} \end{bmatrix} \qquad
\begin{matrix} E & \text{：弾性係数} \\ \nu & \text{：ポアソン比} \end{matrix}
\end{align}

def dmat_pl4(nstr,E,po):
    d=np.zeros((3,3),dtype=np.float64)
    if nstr==0: #plane strain
        d[0,0]=1-po; d[0,1]=po
        d[1,0]=po  ; d[1,1]=1-po
        d[2,2]=0.5*(1-2*po)
        d=E/(1+po)/(1-2*po)*d
    else: # plane stress
        d[0,0]=1  ; d[0,1]=po
        d[1,0]=po ; d[1,1]=1
        d[2,2]=0.5*(1-po)
        d=E/(1-po**2)*d
    return d

6. 関数：Bマトリックスおよびヤコビアン作成 (def bmat_pl4)

\begin{align}
&\cfrac{\partial N_1}{\partial a}=-\cfrac{1}{4}(1-b)  &\cfrac{\partial N_2}{\partial a}=+\cfrac{1}{4}(1-b) & &\cfrac{\partial N_3}{\partial a}=+\cfrac{1}{4}(1+b) & &\cfrac{\partial N_4}{\partial a}=-\cfrac{1}{4}(1+b) \\
&\cfrac{\partial N_1}{\partial b}=-\cfrac{1}{4}(1-a)  &\cfrac{\partial N_2}{\partial b}=-\cfrac{1}{4}(1+a) & &\cfrac{\partial N_3}{\partial b}=+\cfrac{1}{4}(1+a) & &\cfrac{\partial N_4}{\partial b}=+\cfrac{1}{4}(1-a)
\end{align}

\begin{align}
&J_{11}=\frac{\partial N_1}{\partial a} x_i+\frac{\partial N_2}{\partial a} x_j+\frac{\partial N_3}{\partial a} x_k+\frac{\partial N_4}{\partial a} x_l \\
&J_{12}=\frac{\partial N_1}{\partial a} y_i+\frac{\partial N_2}{\partial a} y_j+\frac{\partial N_3}{\partial a} y_k+\frac{\partial N_4}{\partial a} y_l \\
&J_{21}=\frac{\partial N_1}{\partial b} x_i+\frac{\partial N_2}{\partial b} x_j+\frac{\partial N_3}{\partial b} x_k+\frac{\partial N_4}{\partial b} x_l \\
&J_{22}=\frac{\partial N_1}{\partial b} y_i+\frac{\partial N_2}{\partial b} y_j+\frac{\partial N_3}{\partial b} y_k+\frac{\partial N_4}{\partial b} y_l \\
\end{align}

\begin{equation}
\det(J)=J_{11}\cdot J_{22}-J_{12}\cdot J_{21}
\end{equation}

\begin{equation}
[\boldsymbol{B}]
=\begin{bmatrix}
\cfrac{\partial N_1}{\partial x} & 0 & \cfrac{\partial N_2}{\partial x} & 0 & \cfrac{\partial N_3}{\partial x} & 0 & \cfrac{\partial N_4}{\partial x} & 0 \\
0 & \cfrac{\partial N_1}{\partial y} & 0 & \cfrac{\partial N_2}{\partial y} & 0 & \cfrac{\partial N_3}{\partial y} & 0 & \cfrac{\partial N_4}{\partial y} \\
\cfrac{\partial N_1}{\partial y} & \cfrac{\partial N_1}{\partial x} & \cfrac{\partial N_2}{\partial y} & \cfrac{\partial N_2}{\partial x} & \cfrac{\partial N_3}{\partial y} & \cfrac{\partial N_3}{\partial x} & \cfrac{\partial N_4}{\partial y} & \cfrac{\partial N_4}{\partial x}
\end{bmatrix}
\end{equation}

\begin{equation}
\cfrac{\partial N_i}{\partial x}=\cfrac{1}{\det(J)}\left\{J_{22}\cfrac{\partial N_i}{\partial a}-J_{12}\cfrac{\partial N_i}{\partial b}\right\} \qquad
\cfrac{\partial N_i}{\partial y}=\cfrac{1}{\det(J)}\left\{-J_{21}\cfrac{\partial N_i}{\partial a}+J_{11}\cfrac{\partial N_i}{\partial b}\right\}
\end{equation}

def bmat_pl4(a,b,x1,y1,x2,y2,x3,y3,x4,y4):
    bm=np.zeros((3,8),dtype=np.float64)
    #dN/da,dN/db
    dn1a=-0.25*(1.0-b)
    dn2a= 0.25*(1.0-b)
    dn3a= 0.25*(1.0+b)
    dn4a=-0.25*(1.0+b)
    dn1b=-0.25*(1.0-a)
    dn2b=-0.25*(1.0+a)
    dn3b= 0.25*(1.0+a)
    dn4b= 0.25*(1.0-a)
    #Jacobi matrix and det(J)
    J11=dn1a*x1+dn2a*x2+dn3a*x3+dn4a*x4
    J12=dn1a*y1+dn2a*y2+dn3a*y3+dn4a*y4
    J21=dn1b*x1+dn2b*x2+dn3b*x3+dn4b*x4
    J22=dn1b*y1+dn2b*y2+dn3b*y3+dn4b*y4
    detJ=J11*J22-J12*J21
    #[B]=[dN/dx][dN/dy]
    bm[0,0]= J22*dn1a-J12*dn1b
    bm[0,2]= J22*dn2a-J12*dn2b
    bm[0,4]= J22*dn3a-J12*dn3b
    bm[0,6]= J22*dn4a-J12*dn4b
    bm[1,1]=-J21*dn1a+J11*dn1b
    bm[1,3]=-J21*dn2a+J11*dn2b
    bm[1,5]=-J21*dn3a+J11*dn3b
    bm[1,7]=-J21*dn4a+J11*dn4b
    bm[2,0]=-J21*dn1a+J11*dn1b
    bm[2,1]= J22*dn1a-J12*dn1b
    bm[2,2]=-J21*dn2a+J11*dn2b
    bm[2,3]= J22*dn2a-J12*dn2b
    bm[2,4]=-J21*dn3a+J11*dn3b
    bm[2,5]= J22*dn3a-J12*dn3b
    bm[2,6]=-J21*dn4a+J11*dn4b
    bm[2,7]= J22*dn4a-J12*dn4b
    bm=bm/detJ
    return bm,detJ

7. 関数：ガウス積分点指定 (defab_pl4)

プログラム内では，$kk$ の値は，$0\sim 3$ である．

\begin{align}
& i  & j & &  a              & &  b              & & H              & & kk \\
& 1  & 1 & & -0.57735 02692  & & -0.57735 02692  & & 1.00000 00000  & & 1\quad (0) \\
& 1  & 2 & & +0.57735 02692  & & -0.57735 02692  & & 1.00000 00000  & & 2\quad (1) \\
& 2  & 1 & & +0.57735 02692  & & +0.57735 02692  & & 1.00000 00000  & & 3\quad (2) \\
& 2  & 2 & & -0.57735 02692  & & +0.57735 02692  & & 1.00000 00000  & & 4\quad (3)
\end{align}

def ab_pl4(kk):
    if kk==0:
        a=-0.5773502692
        b=-0.5773502692
    if kk==1:
        a= 0.5773502692
        b=-0.5773502692
    if kk==2:
        a= 0.5773502692
        b= 0.5773502692
    if kk==3:
        a=-0.5773502692
        b= 0.5773502692
    return a,b

8. 関数：要素剛性マトリックス作成 (def sm_pl4)

\begin{align}
[\boldsymbol{k}]=&\int_V[\boldsymbol{B}]^T[\boldsymbol{D_e}][\boldsymbol{B}]dV \\
=&t\int_{-1}^{1}\int_{-1}^{1}[\boldsymbol{B}]^{T}[\boldsymbol{D_e}][\boldsymbol{B}]\cdot\det(J)\cdot da\cdot db \\
=&t\cdot \sum_{kk=1}^4 \left\{[\boldsymbol{B}]^{T}[\boldsymbol{D_e}][\boldsymbol{B}]\cdot\det(J)\right\}_{kk}
\end{align}

def sm_pl4(nstr,t,E,po,x1,y1,x2,y2,x3,y3,x4,y4):
    sm=np.zeros((8,8),dtype=np.float64)
    #Stiffness matrix [sm]=[B]T[D][B]*t*det(J)
    d=dmat_pl4(nstr,E,po)
    for kk in range(0,4):
        a,b=ab_pl4(kk)
        bm,detJ=bmat_pl4(a,b,x1,y1,x2,y2,x3,y3,x4,y4)
        sm=sm+np.dot(bm.T,np.dot(d,bm))*t*detJ
    return sm

9. 関数：要素応力計算 (def calst_pl4)

\begin{equation}
\{\boldsymbol{\epsilon}\}=[\boldsymbol{B}]\{\boldsymbol{u_{nd}}\} \qquad \text{（節点変位から求まるひずみ）}
\end{equation}

\begin{align}
N_1=&\frac{1}{4}(1-a)(1-b) \qquad N_2=\frac{1}{4}(1+a)(1-b) \\
N_3=&\frac{1}{4}(1+a)(1+b) \qquad N_4=\frac{1}{4}(1-a)(1+b)
\end{align}

\begin{align}
&T_p=N_1\cdot\phi_i+N_2\cdot\phi_j+N_3\cdot\phi_k+N_4\cdot\phi_l & &\text{（要素内任意点の温度変化量）}\\
&\{\boldsymbol{\epsilon_0}\}=\begin{Bmatrix} (1+\nu)\cdot \alpha\cdot T_p \\ (1+\nu)\cdot \alpha\cdot T_p \\ 0 \end{Bmatrix} & &\text{（平面ひずみ）}\\
&\{\boldsymbol{\epsilon_0}\}=\begin{Bmatrix} \alpha\cdot T_p \\ \alpha\cdot T_p \\ 0 \end{Bmatrix} & &\text{（平面応力）}\\
&\{\boldsymbol{\sigma}\}=[\boldsymbol{D_e}]\{\boldsymbol{\epsilon}-\boldsymbol{\epsilon_0}\} & &\text{（要素応力ベクトル）}
\end{align}

def calst_pl4(nstr,E,po,alpha,dt1,dt2,dt3,dt4,wd,x1,y1,x2,y2,x3,y3,x4,y4):
    eps0  =np.zeros(3,dtype=np.float64)
    stress=np.zeros(3,dtype=np.float64)
    #stress vector {stress}=[D][B]{u}
    d=dmat_pl4(nstr,E,po)
    for kk in range(0,4):
        a,b=ab_pl4(kk)
        bm,detJ=bmat_pl4(a,b,x1,y1,x2,y2,x3,y3,x4,y4)
        eps=np.dot(bm,wd)
        # Strain due to temperature change
        nm1=0.25*(1.0-a)*(1.0-b)
        nm2=0.25*(1.0+a)*(1.0-b)
        nm3=0.25*(1.0+a)*(1.0+b)
        nm4=0.25*(1.0-a)*(1.0+b)
        tem=nm1*dt1+nm2*dt2+nm3*dt3+nm4*dt4
        #Thermal strain
        if nstr==0: # plane strain
            eps0[0]=tem*(1.0+po)*alpha
            eps0[1]=eps0[0]
            eps0[2]=0.0
        else: # plane stress
            eps0[0]=tem*alpha
            eps0[1]=eps0[0]
            eps0[2]=0.0
        stress=stress+np.dot(d,(eps-eps0))
    stress=0.25*stress
    return stress

10. 関数：主応力計算 (def pst_pl4)

主応力は次式で計算する．なお応力の符号は引張が正である．

\begin{align}
&\text{第一主応力} & &\sigma_1=\cfrac{\sigma_x+\sigma_y}{2}+\sqrt{\left(\cfrac{\sigma_x-\sigma_y}{2}\right)^2+\tau_{xy}{}^2} \\
&\text{第二主応力} & &\sigma_2=\cfrac{\sigma_x-\sigma_y}{2}-\sqrt{\left(\cfrac{\sigma_x-\sigma_y}{2}\right)^2+\tau_{xy}{}^2}
\end{align}

第一主応力の方向（x軸から反時計回りの角度）は，一般には次式で計算する．場合分けについてはプログラムを参照．

\begin{equation}
\theta=\cfrac{1}{2}\tan^{-1}\left(\cfrac{2 \tau_{xy}}{\sigma_x-\sigma_y}\right)
\end{equation}

def pst_pl4(sigx,sigy,tauxy):
    ps1=0.5*(sigx+sigy)+np.sqrt(0.25*(sigx-sigy)*(sigx-sigy)+tauxy*tauxy)
    ps2=0.5*(sigx+sigy)-np.sqrt(0.25*(sigx-sigy)*(sigx-sigy)+tauxy*tauxy)
    if sigx==sigy:
        if tauxy >0.0: ang= 45.0
        if tauxy <0.0: ang=135.0
        if tauxy==0.0: ang=  0.0
    else:
        ang=0.5*np.arctan(2.0*tauxy/(sigx-sigy))
        ang=180.0/np.pi*ang
        if sigx>sigy and tauxy>=0.0: ang=ang
        if sigx>sigy and tauxy <0.0: ang=ang+180.0
        if sigx<sigy:                ang=ang+90.0
    return ps1,ps2,ang

11. 関数：節点温度荷重ベクトル作成 (def tfvec_pl4)

\begin{align}
N_1=&\frac{1}{4}(1-a)(1-b) \qquad N_2=\frac{1}{4}(1+a)(1-b) \\
N_3=&\frac{1}{4}(1+a)(1+b) \qquad N_4=\frac{1}{4}(1-a)(1+b)
\end{align}

\begin{equation}
T_p=N_1\cdot\phi_i+N_2\cdot\phi_j+N_3\cdot\phi_k+N_4\cdot\phi_l
\end{equation}

\begin{align}
&\{\boldsymbol{\epsilon_0}\}=\begin{Bmatrix} (1+\nu)\cdot \alpha\cdot T_p \\ (1+\nu)\cdot \alpha\cdot T_p \\ 0 \end{Bmatrix} & &\text{（平面ひずみ）}\\
&\{\boldsymbol{\epsilon_0}\}=\begin{Bmatrix} \alpha\cdot T_p \\ \alpha\cdot T_p \\ 0 \end{Bmatrix} & &\text{（平面応力）}
\end{align}

\begin{align}
\{\boldsymbol{f_t}\}=&\int_V[\boldsymbol{B}]^T[\boldsymbol{D_e}]\{\boldsymbol{\epsilon_0}\}dV \\
=&t\int_{-1}^{1}\int_{-1}^{1}[\boldsymbol{B}]^{T}[\boldsymbol{D_e}]\{\boldsymbol{\epsilon_0}\}\cdot\det(J)\cdot da\cdot db \\
=&t\cdot \sum_{kk=1}^4 \left\{[\boldsymbol{B}]^{T}[\boldsymbol{D_e}]\{\boldsymbol{\epsilon_0}\}\cdot\det(J)\right\}_{kk}
\end{align}

def tfvec_pl4(nstr,t,E,po,alpha,dt1,dt2,dt3,dt4,x1,y1,x2,y2,x3,y3,x4,y4):
    eps0=np.zeros(3,dtype=np.float64)
    tfe =np.zeros(8,dtype=np.float64)
    # {tfe=[B]T[D]{eps0}
    d=dmat_pl4(nstr,E,po)
    for kk in range(0,4):
        a,b=ab_pl4(kk)
        bm,detJ=bmat_pl4(a,b,x1,y1,x2,y2,x3,y3,x4,y4)
        # thermal Strain
        nm1=0.25*(1.0-a)*(1.0-b)
        nm2=0.25*(1.0+a)*(1.0-b)
        nm3=0.25*(1.0+a)*(1.0+b)
        nm4=0.25*(1.0-a)*(1.0+b)
        tem=nm1*dt1+nm2*dt2+nm3*dt3+nm4*dt4
        if nstr==0: # plane strain
            eps0[0]=tem*(1.0+po)*alpha
            eps0[1]=eps0[0]
            eps0[2]=0.0
        else: # plane stress
            eps0[0]=tem*alpha
            eps0[1]=eps0[0]
            eps0[2]=0.0
        tfe=tfe+np.dot(bm.T,np.dot(d,eps0))*t*detJ
    return tfe

12. 関数：節点慣性力ベクトル作成 (def bfvec_pl4)

この関数では，ひずみー変位関係行列 $[\boldsymbol{B}]$ を呼ばないため，$\det(J)$ をこの関数内で計算している．

\begin{align}
&\cfrac{\partial N_1}{\partial a}=-\cfrac{1}{4}(1-b)  &\cfrac{\partial N_2}{\partial a}=+\cfrac{1}{4}(1-b) & &\cfrac{\partial N_3}{\partial a}=+\cfrac{1}{4}(1+b) & &\cfrac{\partial N_4}{\partial a}=-\cfrac{1}{4}(1+b) \\
&\cfrac{\partial N_1}{\partial b}=-\cfrac{1}{4}(1-a)  &\cfrac{\partial N_2}{\partial b}=-\cfrac{1}{4}(1+a) & &\cfrac{\partial N_3}{\partial b}=+\cfrac{1}{4}(1+a) & &\cfrac{\partial N_4}{\partial b}=+\cfrac{1}{4}(1-a)
\end{align}

\begin{align}
&J_{11}=\frac{\partial N_1}{\partial a} x_i+\frac{\partial N_2}{\partial a} x_j+\frac{\partial N_3}{\partial a} x_k+\frac{\partial N_4}{\partial a} x_l \\
&J_{12}=\frac{\partial N_1}{\partial a} y_i+\frac{\partial N_2}{\partial a} y_j+\frac{\partial N_3}{\partial a} y_k+\frac{\partial N_4}{\partial a} y_l \\
&J_{21}=\frac{\partial N_1}{\partial b} x_i+\frac{\partial N_2}{\partial b} x_j+\frac{\partial N_3}{\partial b} x_k+\frac{\partial N_4}{\partial b} x_l \\
&J_{22}=\frac{\partial N_1}{\partial b} y_i+\frac{\partial N_2}{\partial b} y_j+\frac{\partial N_3}{\partial b} y_k+\frac{\partial N_4}{\partial b} y_l
\end{align}

\begin{equation}
\det(J)=J_{11}\cdot J_{22}-J_{12}\cdot J_{21}
\end{equation}

\begin{align}
N_1=&\frac{1}{4}(1-a)(1-b) \qquad N_2=\frac{1}{4}(1+a)(1-b) \\
N_3=&\frac{1}{4}(1+a)(1+b) \qquad N_4=\frac{1}{4}(1-a)(1+b)
\end{align}

\begin{align}
\{\boldsymbol{f_b}\}=&\gamma \cdot \left(\int_V [\boldsymbol{N}]^T [\boldsymbol{N}]dV\right)\{\boldsymbol{w_{nd}}\} \\
=&\gamma\cdot t\cdot\int_{A}[\boldsymbol{N}]^T[\boldsymbol{N}]dA\cdot\{\boldsymbol{w_{nd}}\} \\
=&\gamma\cdot t \cdot\sum_{kk=1}^4 \left\{[\boldsymbol{N}]^T[\boldsymbol{N}]\cdot\det(J)\right\}_ {kk}\cdot\{\boldsymbol{w_{nd}}\}
\end{align}

\begin{equation}
\{\boldsymbol{w_{nd}}\}=\begin{Bmatrix}k_h & k_v & k_h & k_v & k_h & k_v & k_h & k_v \end{Bmatrix}^T
\end{equation}

def bfvec_pl4(t,gamma,gkh,gkv,x1,y1,x2,y2,x3,y3,x4,y4):
    ntn=np.zeros((8,8),dtype=np.float64)
    nm =np.zeros((2,8),dtype=np.float64)
    for kk in range(0,4):
        a,b=ab_pl4(kk)
        #dN/da,dN/db
        dn1a=-0.25*(1.0-b)
        dn2a= 0.25*(1.0-b)
        dn3a= 0.25*(1.0+b)
        dn4a=-0.25*(1.0+b)
        dn1b=-0.25*(1.0-a)
        dn2b=-0.25*(1.0+a)
        dn3b= 0.25*(1.0+a)
        dn4b= 0.25*(1.0-a)
        #Jacobi matrix and det(J)
        J11=dn1a*x1+dn2a*x2+dn3a*x3+dn4a*x4
        J12=dn1a*y1+dn2a*y2+dn3a*y3+dn4a*y4
        J21=dn1b*x1+dn2b*x2+dn3b*x3+dn4b*x4
        J22=dn1b*y1+dn2b*y2+dn3b*y3+dn4b*y4
        detJ=J11*J22-J12*J21
        nm[0,0]=0.25*(1.0-a)*(1.0-b)
        nm[0,2]=0.25*(1.0+a)*(1.0-b)
        nm[0,4]=0.25*(1.0+a)*(1.0+b)
        nm[0,6]=0.25*(1.0-a)*(1.0+b)
        nm[1,1]=nm[0,0]
        nm[1,3]=nm[0,2]
        nm[1,5]=nm[0,4]
        nm[1,7]=nm[0,6]
        ntn=ntn+np.dot(nm.T,nm)*gamma*t*detJ
    w=np.array([gkh,gkv,gkh,gkv,gkh,gkv,gkh,gkv],dtype=np.float64)
    bfe=np.dot(ntn,w)
    return bfe

13. 関数：メインルーチン

メインルーチン処理手順

(1) コマンドラインより入出力ファイル名取得

コマンドライン引数として入出力データファイル名を取得．
入力データファイル名：fnameR，出力データファイル名：fnameW

(2) データ入力

関数 inpdata_pl4 により入力データ取得

(3) 入力データ書き出し

関数 prinp_pl4 により入力データ書き出し

(4) 剛性マトリックス組み立て

各要素毎に，剛性マトリックスsm，温度荷重ベクトル tfe，慣性力ベクトル bfe を計算し，全体剛性方程式に組み込んでいく．ここで剛性マトリックス，温度荷重ベクトル，慣性力ベクトル計算の関数の引数は，全てスカラーとした．要素番号が指定されれば，要素を構成する節点座標や要素物性値は特定できるため，他の要素の値を含んで収納する配列よりは，スカラーを送り込んだほうがプログラム的にスッキリするように思う．
全体剛性マトリックスは変数 gk に，荷重ベクトルは変数 fp に加算されていく．
要素剛性マトリックスを全体剛性マトリックスに組み込むための位置情報は，配列 ir[] に格納されている．

8 x 8 （8 は１節点自由度２と１要素節点数４の積）の次元の要素剛性マトリックスを NｘN （N は１節点自由度２と総節点数の積）の次元の全体剛性マトリックスにどのように配置するかを示したのが下図である．
碁盤目に区切られた大きな正方形は全体剛性マトリックス（N x N）を示し，$i$，$j$，$k$，$l$ は要素を構成する節点番号，$k_{11}\sim k_{88}$ は要素剛性マトリックス（8 x 8）の成分を示す．

(5) 境界条件処理（変位既知自由度の処理）

全体剛性マトリックスの大きさは変えず，変位既知節点の処理を行う．

変数	説明
npoin	総節点数
nfree	１節点自由度（ここでは 2）
mpfix	各節点毎の拘束条件を格納した配列（=1：拘束，=0：自由）
fp	節点荷重ベクトル
rdis	変位拘束する節点の変位量
gk	全体剛性マトリックス
disg	全体剛性方程式の解（変位）

(6) 剛性方程式を解く（変位計算）

csr_matrix により剛性マトリックスを疎行列として圧縮格納した後，spsolve にて解いている．

(7) 境界条件処理（既知変位の組み込み）

連立一次方程式を解いた後，既知変位に相当する自由度の項に，既知変位を入れ直すことを忘れないようにする．

(8) 要素応力計算

関数 calst_pl4 による要素毎の応力計算実施

(9) 計算結果書き出し

関数 prout_pl4 による計算結果出力

(10) 計算時間等情報表示

以下により，計算完了時刻とスタート時刻の差分を処理時間として出力する．

import time

start=time.time()
......
dtime=time.time()-start

メインルーチンプログラム

def main_pl4():
    start=time.time()
    args = sys.argv
    fnameR=args[1] # input data file
    fnameW=args[2] # output data file
    nod=4   # Number of nodes per element
    nfree=2 # Degree of freedom per node
    # data input
    npoin,nele,nsec,npfix,nlod,nstr,ae,node,x,deltaT,mpfix,rdis,fp=inpdata_pl4(fnameR,nod,nfree)
    # print out of input data
    prinp_pl4(fnameW,npoin,nele,nsec,npfix,nlod,nstr,ae,node,x,deltaT,mpfix,rdis,fp)
    # array declaration
    ir=np.zeros(nod*nfree,dtype=np.int)         # Work vector for matrix assembly
    gk=np.zeros((nfree*npoin,nfree*npoin),dtype=np.float64) # Global stiffness matrix
    # assembly of stiffness matrix & load vector
    for ne in range(0,nele):
        i=node[0,ne]-1
        j=node[1,ne]-1
        k=node[2,ne]-1
        l=node[3,ne]-1
        m=node[4,ne]-1
        x1=x[0,i]; y1=x[1,i]
        x2=x[0,j]; y2=x[1,j]
        x3=x[0,k]; y3=x[1,k]
        x4=x[0,l]; y4=x[1,l]
        t    =ae[0,m]
        E    =ae[1,m]
        po   =ae[2,m]
        alpha=ae[3,m]
        gamma=ae[4,m]
        gkh  =ae[5,m]
        gkv  =ae[6,m]
        dt1=deltaT[i]; dt2=deltaT[j]; dt3=deltaT[k]; dt4=deltaT[l] # temperature change
        sm=sm_pl4(nstr,t,E,po,x1,y1,x2,y2,x3,y3,x4,y4)             # element stiffness matrix
        tfe=tfvec_pl4(nstr,t,E,po,alpha,dt1,dt2,dt3,dt4,x1,y1,x2,y2,x3,y3,x4,y4) # thermal load vector
        bfe=bfvec_pl4(t,gamma,gkh,gkv,x1,y1,x2,y2,x3,y3,x4,y4)     # body force
        ir[7]=2*l+1; ir[6]=ir[7]-1
        ir[5]=2*k+1; ir[4]=ir[5]-1
        ir[3]=2*j+1; ir[2]=ir[3]-1
        ir[1]=2*i+1; ir[0]=ir[1]-1
        for i in range(0,nod*nfree):
            it=ir[i]
            fp[it]=fp[it]+tfe[i]+bfe[i]
            for j in range(0,nod*nfree):
                jt=ir[j]
                gk[it,jt]=gk[it,jt]+sm[i,j]
    # treatment of boundary conditions
    for i in range(0,npoin):
        for j in range(0,nfree):
            if mpfix[j,i]==1:
                iz=i*nfree+j
                fp[iz]=0.0
    for i in range(0,npoin):
        for j in range(0,nfree):
            if mpfix[j,i]==1:
                iz=i*nfree+j
                fp=fp-rdis[j,i]*gk[:,iz]
                gk[:,iz]=0.0
                gk[iz,iz]=1.0
    # solution of simultaneous linear equations
    #disg = np.linalg.solve(gk, fp)
    gk = csr_matrix(gk)
    disg = spsolve(gk, fp, use_umfpack=True)
    fp   = np.zeros(nfree*npoin,dtype=np.float64)
    # recovery of restricted displacements
    for i in range(0,npoin):
        for j in range(0,nfree):
            if mpfix[j,i]==1:
                iz=i*nfree+j
                fp[iz]=-disg[iz]
                disg[iz]=rdis[j,i]
    # calculation of element stress
    wd    =np.zeros(8,dtype=np.float64)
    strs  =np.zeros((3,nele),dtype=np.float64)      # Section force vector
    for ne in range(0,nele):
        i=node[0,ne]-1
        j=node[1,ne]-1
        k=node[2,ne]-1
        l=node[3,ne]-1
        m=node[4,ne]-1
        x1=x[0,i]; y1=x[1,i]
        x2=x[0,j]; y2=x[1,j]
        x3=x[0,k]; y3=x[1,k]
        x4=x[0,l]; y4=x[1,l]
        wd[0]=disg[2*i]; wd[1]=disg[2*i+1]
        wd[2]=disg[2*j]; wd[3]=disg[2*j+1]
        wd[4]=disg[2*k]; wd[5]=disg[2*k+1]
        wd[6]=disg[2*l]; wd[7]=disg[2*l+1]
        E    =ae[1,m]
        po   =ae[2,m]
        alpha=ae[3,m]
        dt1=deltaT[i]; dt2=deltaT[j]; dt3=deltaT[k]; dt4=deltaT[l]
        strs[:,ne]=calst_pl4(nstr,E,po,alpha,dt1,dt2,dt3,dt4,wd,x1,y1,x2,y2,x3,y3,x4,y4)
    # print out of result
    prout_pl4(fnameW,npoin,nele,disg,fp,strs)
    # information
    dtime=time.time()-start
    print('n={0}  time={1:.3f}'.format(nfree*npoin,dtime)+' sec')
    fout=open(fnameW,'a')
    print('n={0}  time={1:.3f}'.format(nfree*npoin,dtime)+' sec',file=fout)
    fout.close()

14. メインルーチン実行

#==============
# Execution
#==============
if __name__ == '__main__': main_pl4()

以　上