2017.12.30 投稿
2018.08.11 プログラム改良

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

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

はじめに

力学系の有限要素法の次に，２次元浸透流解析のプログラムの解説を作ってみた．本来私は土木構造物の構造設計が本職であるが，水力発電所土木設計関係の仕事をする関係上，浸透流解析や熱伝導解析を用いた仕事も行ってきた．ここに示す解説とプログラムは，それらの経験に基づき作成したものである．対象とする問題は，ダムや矢板の周りの浸透流，トンネルや地下空洞周辺の圧力水頭分布であり，浸潤面境界への流入（浸潤面境界への雨等の流入）は考慮していない．ただし，トンネル周辺領域への地下水の流入は，解析領域周囲の圧力水頭指定や流量指定により考慮できる．

これは「私見」であるが，有限要素法のプログラムは，線形問題であれば，力学問題にしろ浸透流のようなポテンシャル問題にしろ，基礎微分方程式をルールに従って連立一次方程式を解く問題に帰着させることなので，それほど違いがあるものではない．

この解説では，まず２次元飽和定常浸透流解析の理論を解説し，そのあと飽和ー不飽和解析への拡張という説明の進め方をしている．
解析プログラム上，不飽和領域を含む計算は，透水係数を水頭値の関数として収束計算を行うものであり，飽和解析と比較して，それほど大々的なプログラム変更が必要なものではないという印象である．

また，ここでは三角形一次要素による定式化を行っているが，１節点１自由度なので，２次元応力解析プログラムにおける扱いよりも定式化自体は単純である．とはいえ，三角形一次要素では，２次元問題でも浸透流問題でも，剛性マトリックスあるいは透水性マトリックスは解析的に積分され定数の集まりに集約されるので，プログラム上は扱いが単純でうれしい．

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

２次元定常飽和浸透流解析の理論

Darcyの法則が成立する2次元直交異方性材料の直交座標系 $x$-$y$ における定常飽和浸透流問題の支配方程式は以下の通りである．なお奥行き方向の厚さは考慮しない．

\begin{equation}
k_x\frac{\partial^2 \phi}{\partial x^2}+k_y\frac{\partial^2 \phi}{\partial y^2}=0
\end{equation}

\begin{align}
&\phi & &\text{：全水頭}\\
&k_x  & &\text{：$x$方向透水係数}\\\
&k_y  & &\text{：$y$方向透水係数}
\end{align}

ここで，Galerkin法により，任意の $\delta \phi$ に対し，上式の弱形式である次式を満足する未知数 $\phi$ を決定することを考える．

\begin{equation}
\int_A\delta \phi\left(k_x\frac{\partial^2 \phi}{\partial x^2}+k_y\frac{\partial^2 \phi}{\partial y^2}\right)dA=0
\end{equation}

有限要素において，要素内の任意箇所の全水頭 $\phi$ が，内挿関数 $[\boldsymbol{N}]$ と節点全水頭ベクトル ${\boldsymbol{\phi}}$ を用いて，

\begin{equation}
\phi(x,y)=[\boldsymbol{N}(x,y)]\{\boldsymbol{\phi}\}
\end{equation}

と表せれば，弱形式は，

\begin{equation}
\{\boldsymbol{\delta \phi}\}^T\int_A [\boldsymbol{N}]^T \left(k_x\frac{\partial^2 \phi}{\partial x^2}+k_y\frac{\partial^2 \phi}{\partial y^2}\right)dA=0
\end{equation}

となるが，${\boldsymbol{\delta \phi}}$ は任意の節点水頭ベクトルであるため，面積積分値が $0$ となる必要がある．

次に部分積分公式

\begin{equation*}
f\cdot g'=(f\cdot g)'-f'\cdot g
\end{equation*}

およびGreenの公式

\begin{equation}
\iint_D\left(\frac{\partial P}{\partial x}+\frac{\partial Q}{\partial y}\right)dA=\int_C \left(P\frac{\partial x}{\partial n}+Q\frac{\partial y}{\partial n}\right)ds
\end{equation}

を用いることにより，

\begin{align}
0=&\int_A [\boldsymbol{N}]^T \left(k_x\frac{\partial^2 \phi}{\partial x^2}+k_y\frac{\partial^2 \phi}{\partial y^2}\right)dA \\
=&\int_A 
\left\{
k_x\left(
\frac{\partial \left([\boldsymbol{N}]^T\frac{\partial \phi}{\partial x}\right)}{\partial x}-\frac{\partial [\boldsymbol{N}]^T}{\partial x}\frac{\partial \phi}{\partial x}
\right)
+k_y\left(
\frac{\partial \left([\boldsymbol{N}]^T\frac{\partial \phi}{\partial y}\right)}{\partial y}-\frac{\partial [\boldsymbol{N}]^T}{\partial y}\frac{\partial \phi}{\partial y}
\right)
\right\}
dA \\
=&\int_s [\boldsymbol{N}]^T \left(k_x\frac{\partial \phi}{\partial x}\frac{\partial x}{\partial n}+k_y\frac{\partial \phi}{\partial y}\frac{\partial y}{\partial n}\right)ds
-\int_A\left(k_x\frac{\partial [\boldsymbol{N}]^T}{\partial x}\frac{\partial [\boldsymbol{N}]}{\partial x}+k_y\frac{\partial [\boldsymbol{N}]^T}{\partial y}\frac{\partial [\boldsymbol{N}]}{\partial y}\right)dA\cdot\{\boldsymbol{\phi}\}
\end{align}

また，$x$ 方向および $y$ 方向の流速をそれぞれ $v_x$，$v_y$ とすれば，Darcyの法則より

\begin{equation}
v_x=-k_x\frac{\partial \phi}{\partial x} \qquad v_y=-k_x\frac{\partial \phi}{\partial y}
\end{equation}

以上により，要素有限要素式は以下の通り表現できる．

\begin{align}
&[\boldsymbol{k}]\{\boldsymbol{\phi}\}=\{\boldsymbol{q}\} \\
&[\boldsymbol{k}]=\int_A \left(k_x\frac{\partial [\boldsymbol{N}]^T}{\partial x}\frac{\partial [\boldsymbol{N}]}{\partial x}+k_y\frac{\partial [\boldsymbol{N}]^T}{\partial y}\frac{\partial [\boldsymbol{N}]}{\partial y}\right)dA \\
&\{\boldsymbol{q}\}=-\int_s [\boldsymbol{N}]^T \left(v_x \frac{\partial x}{\partial n}+v_y \frac{\partial y}{\partial n}\right)ds
\end{align}

\begin{align}
&[\boldsymbol{k}]   & &\text{：透水性マトリックス}\\
&\{\boldsymbol{\phi}\} & &\text{：節点全水頭ベクトル}\\
&\{\boldsymbol{q}\} & &\text{：節点流量ベクトル}
\end{align}

なお ${\boldsymbol{q}}$ は流速ベクトルの境界長さでの積分であり，流量を示すことは容易に理解できる．

三角形一次要素の導入

要素の有限要素式

２次元浸透流解析用の三角形要素は，１要素４節点であり，１節点１自由度（節点全水頭もしくは節点流量）を有する．
よって要素の有限要素式の次元は３（１節点１自由度ｘ１要素３節点）であり，具体的には以下の形をとる．
ここに，三角形要素を構成する節点番号を反時計回りに，$i$，$j$，$k$ とする．

\begin{equation}
\begin{Bmatrix} q_i \\ q_j \\ q_k  \end{Bmatrix}
=\begin{bmatrix}
k_{11} & k_{12} & k_{13}  \\
k_{21} & k_{22} & k_{23}  \\
k_{31} & k_{32} & k_{33} 
\end{bmatrix}
\begin{Bmatrix} \phi_i \\ \phi_j \\ \phi_k \end{Bmatrix}
\end{equation}

\begin{align}
&q_i & &\text{：節点$i$の流量} & & \phi_i & &\text{：節点$i$の全水頭}\\
&q_j & &\text{：節点$j$の流量} & & \phi_j & &\text{：節点$j$の全水頭}\\
&q_k & &\text{：節点$k$の流量} & & \phi_k & &\text{：節点$k$の全水頭}
\end{align}

内挿関数

２次元浸透流解析に用いる三角形要素では，１節点１自由度３節点であるため，１要素あたり３自由度を有する．
このため，３個の未定係数 $\alpha_{1\sim3}$ を用い，三角形要素内任意点の全水頭 $\phi$ を要素内座標 ($x,y$) により以下の通り仮定する．

\begin{align}
\phi=&\alpha_1+\alpha_2 x+\alpha_3 y
=\begin{bmatrix}1 & x & y\end{bmatrix}
\begin{Bmatrix}
\alpha_1 \\
\alpha_2 \\
\alpha_3
\end{Bmatrix}
\end{align}

上式より，各節点の全水頭は，各節点座標を($x_i, y_i$)，($x_j, y_j$)，($x_k, y_k$)とすれば以下の通りとなる．

\begin{equation*}
\begin{Bmatrix}
\phi_i \\
\phi_j \\
\phi_k
\end{Bmatrix}=
\begin{bmatrix}
1 & x_i & y_i \\
1 & x_j & y_j \\
1 & x_k & y_k
\end{bmatrix}
\begin{Bmatrix}
\alpha_1 \\
\alpha_2 \\
\alpha_3
\end{Bmatrix}
\end{equation*}

これを ${\alpha}$ について解くと

\begin{equation}
\begin{Bmatrix}
\alpha_1 \\
\alpha_2 \\
\alpha_3
\end{Bmatrix}=
\cfrac{1}{2\Delta}
\begin{bmatrix}
a_i & a_j & a_k \\
b_i & b_j & b_k \\
c_i & c_j & c_k
\end{bmatrix}
\begin{Bmatrix}
\phi_i \\
\phi_j \\
\phi_k
\end{Bmatrix}
\end{equation}

\begin{equation}
\begin{matrix}
a_i=x_j y_k-x_k y_j & \qquad & b_i=y_j-y_k & \qquad & c_i=x_k-x_j \\
a_j=x_k y_i-x_i y_k & \qquad & b_j=y_k-y_i & \qquad & c_j=x_i-x_k \\
a_k=x_i y_j-x_j y_i & \qquad & b_k=y_i-y_j & \qquad & c_k=x_j-x_i
\end{matrix}
\end{equation}

\begin{equation}
2\Delta=(x_k-x_j)y_i+(x_i-x_k)y_j+(x_j-x_i)y_k
\end{equation}

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

\begin{align}
[\boldsymbol{N}]=&\begin{bmatrix}1 & x & y\end{bmatrix}\cdot
\cfrac{1}{2\Delta}
\begin{bmatrix}
a_i & a_j & a_k \\
b_i & b_j & b_k \\
c_i & c_j & c_k
\end{bmatrix}\\
=&
\begin{bmatrix}
\cfrac{a_i+b_i x+c_i y}{2\Delta}
&\cfrac{a_j+b_j x+c_j y}{2\Delta}
&\cfrac{a_k+b_k x+c_k y}{2\Delta}
\end{bmatrix}
\end{align}

さらに透水性マトリックス作成のため，内挿関数の微分を求めておく．

\begin{align}
&\frac{\partial[\boldsymbol{N}]}{\partial x}=\left[\frac{b_i}{2\Delta} \quad \frac{b_j}{2\Delta} \quad\frac{b_k}{2\Delta} \right] \
&\frac{\partial[\boldsymbol{N}]}{\partial y}=\left[\frac{c_i}{2\Delta} \quad \frac{c_j}{2\Delta} \quad\frac{c_k}{2\Delta} \right]
\end{align}

透水性マトリックス

\begin{align}
[\boldsymbol{k}]=&\int_A \left(k_x\frac{\partial [\boldsymbol{N}]^T}{\partial x}\frac{\partial [\boldsymbol{N}]}{\partial x}+k_y\frac{\partial [\boldsymbol{N}]^T}{\partial y}\frac{\partial [\boldsymbol{N}]}{\partial y}\right)dA \\
=&\frac{k_x}{4\Delta}
\begin{bmatrix}
b_i b_i & b_j b_i & b_k b_i \\
b_i b_j & b_j b_j & b_k b_j \\
b_i b_k & b_j b_k & b_k b_k
\end{bmatrix}
+\frac{k_y}{4\Delta}
\begin{bmatrix}
c_i c_i & c_j c_i & c_k c_i \\
c_i c_j & c_j c_j & c_k c_j \\
c_i c_k & c_j c_k & c_k c_k
\end{bmatrix}
\end{align}

等方均質透水体の場合

等方均質透水体の場合，$k_x=k_y=k$ として，透水マトリックスは以下の通りとなる．

\begin{equation}
[\boldsymbol{k}]=
\frac{k}{4\Delta}
\begin{bmatrix}
b_i b_i + c_i c_i & b_j b_i + c_j c_i & b_k b_i + c_k c_i \\
b_i b_j + c_i c_j & b_j b_j + c_j c_j & b_k b_j + c_k c_j \\
b_i b_k + c_i c_k & b_j b_k + c_j c_k & b_k b_k + c_k c_k
\end{bmatrix}
\end{equation}

理論では，異方性透水体としたが，プログラムでは等方均質体としている．

要素内平均流速

要素内任意点の全水頭 $\phi$ は内挿関数 $[\boldsymbol{N}]$ と節点全水頭ベクトル${\boldsymbol{\phi}}$を用いて

\begin{equation}
\phi=[\boldsymbol{N}]\{\boldsymbol{\phi}\}
\end{equation}

であることから，Darcyの法則より要素内任意点の流速 $v_x$ および $v_y$ は

\begin{align}
&v_x=-k_x\frac{\partial \phi}{\partial x} =-k_x\frac{\partial [\boldsymbol{N}]}{\partial x}\{\boldsymbol{\phi}\}=-\frac{k_x}{2\Delta}(b_i \phi_i+b_j \phi_j+b_k \phi_k) \\
&v_y=-k_y\frac{\partial \phi}{\partial y} =-k_y\frac{\partial [\boldsymbol{N}]}{\partial y}\{\boldsymbol{\phi}\}=-\frac{k_y}{2\Delta}(c_i \phi_i+c_j \phi_j+c_k \phi_k)
\end{align}

不飽和浸透流解析への拡張

飽和ー不飽和浸透流解析の特徴

飽和浸透流解析では，下表に示すように，透水性マトリックスは定数（材料特性のみに依存）であり，また全水頭が既知の節点は流量が未知，流量が既知の節点は全水頭が未知という関係があるため，これを考慮した有限要素式の修正（境界条件処理）を行うことにより，連立方程式を１回解けばよかった．

飽和解析	節点全水頭	節点流量
水頭指定境界節点	既知	未知
流量指定境界節点	未知	既知
無指定節点	未知	既知（ゼロ）
透水性マトリックス	定数
入力する材料特性	飽和透水係数

飽和ー不飽和浸透流解析では，解析領域の中に飽和領域と不飽和領域が混在し，下表に示すように，透水性マトリックスは全水頭（正確にはサクション：負の圧力水頭）の関数であり，また浸潤面境界では節点全水頭・節点流量の双方が条件付きの未知数（事前に確定していない値）であるため，反復計算が必要となる．

飽和ー不飽和解析	節点全水頭	節点流量
水頭指定境界節点	既知	未知
流量指定境界節点	未知	既知
無指定節点	未知	ゼロ
浸潤面境界節点	未知（条件付き）	未知（条件付き）
透水性マトリックス	全水頭の関数
入力する材料特性	飽和透水係数および不飽和透水特性

地盤の不飽和透水特性

地盤の不飽和透水特性として，Van Genuchten モデルがよく知られている．
このモデルは，不飽和透水特性が連続関数で与えられており，解析上プログラムを単純化でき便利である．

\begin{align}
&S_e=\left\{\cfrac{1}{1+(\alpha\cdot H_s)^n}\right\}^m \qquad\qquad n=\cfrac{1}{1-m} \qquad (0<m<1)\\
&K_r=(S_e)^{0.5}\cdot\left\{1-\left(1-S_e{}^{1/m}\right)^m\right\}^2 \qquad (0<S_e, K_R\leq1)\\
&K=K_r \cdot K_0
\end{align}

\begin{align}
&S_e    & &\text{：有効飽和度}\\
&K_r    & &\text{：比透水係数}\\
&H_S    & &\text{：サクション（負圧が正）}\\
&\alpha & &\text{：不飽和特性を表すパラメータ（正の実数）}\\
&m      & &\text{：不飽和特性を表すパラメータ（０を超え１未満の実数）}\\
&K      & &\text{：不飽和透水係数}\\
&K_0    & &\text{：飽和透水係数}
\end{align}

上記のうち，プログラムの入力データとして必要なのは，飽和透水係数 $K_0$，不飽和特性を表すパラメータ $\alpha$ および $m$ である．
サクション $H_s$ については，計算過程において，全水頭と位置水頭の差として算出できる．ただし上式では負圧を正として扱うことに注意する．

反復計算

反復計算においては，浸潤面境界上節点の条件を補正しながら逐次全水頭ベクトルと透水性マトリックスを修正し，全水頭値が収束するまで反復を繰り返している．
浸潤面境界上節点の拘束条件は以下の通り．

浸潤面が出現する境界では，流出のみが許容される（ここでは雨等の影響は考えない）．
浸潤面が出現する境界は，大気に接する境界であるため，圧力水頭がゼロ以下でなければならない．

反復計算の方法としては，全節点に全水頭の初期値を与え，得られた解を次の反復計算の入力値とする，逐次代入法を用いている．
この方法は，全節点の水頭を仮定するため，全要素の透水マトリックスを，仮定した全水頭から計算でき，計算フローを単純化できる．
全水頭の初期値は，全水頭既知境界以外は，飽和領域の最小全水頭とするのが経験的に有効なようである．

反復計算のおおまかな過程は以下の通りである．

全節点の全水頭を仮定（入力データ）
仮定した全水頭を用いた透水マトリックスの作成と組み立て
流量ベクトル・全水頭ベクトルを確からしい値に補正（全水頭および流量指定値の反映，浸潤面境界の拘束条件反映）
補正した条件を用いた連立方程式の拘束条件処理
連立方程式を解く（全水頭の計算）
全水頭による収束判定
収束しなければ全水頭ベクトルを更新して 2. に戻る．

三角形一次要素を用いた２次元飽和ー不飽和浸透流解析プログラム

プログラムの概要

プログラミング言語はPythonとする
等方均質透水地盤の２次元飽和ー不飽和定常浸透流解析を行う．等方均質地盤であることに注意．
座標システムは，右方向を x 軸正方向，上方向あるいは標高方向を z 軸正方向とする．（理論式では x-y 座標としているため，プログラムでは理論式の y を z と読み替えることに注意）
要素は，１要素３節点の三角形一次要素とする．地盤の奥行きは入力値の単位に応じた単位奥行きとする．
境界条件として，不透水境界，全水頭指定境界，流量指定境界，浸潤面境界が扱える．
浸潤面境界への流入はないものとする（雨等の影響は考えない）
連立一次方程式は，その大きさ（１節点自由度ｘ総節点数）を変えずに，全水頭既知節点処理を行い解いている．
連立一次方程式は，大きな自由度での計算効率を考慮し，疎行列計算ライブラリの scipy.sparse.linalg.spsolve(A,b) で解く
出力は，節点全水頭，節点圧力水頭，節点流量および要素平均流速とする
入力データファイル，出力データファイルの文字区切りは空白とする
入力データファイル名，出力データファイル名は，コマンドライン引数として入力する

扱える境界条件の説明

項目	説明
不透水境界	何も指定しなければ不透水境界となる
全水頭指定境界	既知全水頭を指定する節点数・節点番号・全水頭値を入力する．正の流量が節点への流入，負の流量が節点からの流出を示す．
流量指定境界	既知流量を指定する節点数・節点番号・流量値を入力する
浸潤面境界	浸潤面が現れる可能性のある領域の節点数・節点番号を入力する

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

1. モジュール読み込み
2. 入力データ読み込み (def inpdata_s3)  
    (1) 基本数値読み込み  
    (2) 配列宣言  
    (3) 材料特性データ読み込み  
    (4) 要素構成節点番号および材料番号読み込み  
    (5) 節点座標および節点全水頭初期値読み込み  
    (6) 境界条件読み込み  
3. 関数：入力データ書き出し (def prinp_s3)
4. 関数：計算結果書き出し (def prout_s3)
5. 関数：不飽和透水特性計算 (def cal_kre_s3)
6. 関数：要素透水マトリックス作成 (def sm_s3)
7. 関数：要素流速計算 (def calv_s3)
8. 関数：メインルーチン  
    (1) 入出力ファイル名取得
    (2) 基本数値定義
    (3) 入力データ読み込み
    (4) 入力データ書き出し    
    (5) 反復計算過程  
        a. 透水性マトリックス組み立て  
        b. 節点流量計算  
        c. 作業用流量ベクトル更新  
        d. 境界条件処理（全水頭既知自由度の処理）  
        e. 連立一次方程式を解く（全水頭計算）  
        f. 境界条件処理（既知全水頭の組み込み）  
        g. 収束判定  
    (6) 要素流速ベクトル計算  
    (7) 計算結果書き出し  
9. メインルーチン実行

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

解析実行コマンド

python3 py_fem_seep3x.py inp.txt out.txt

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

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

npoin  nele  nsec  koh  koq  kou  idan  # Basic values for analysis
Ak0  alpha  em                          # Material properties
    ..... (1 to nsec) .....             #
node-1  node-2  node-3  isec            # Element connectivity
    ..... (1 to nele) .....             #   (counter-clockwise order of node number)
x  z  hvec0                             # Coordinate, initial value of total head
    ..... (1 to npoin) .....            #
nokh  Hinp                              # Node with given total head, given total head
    ..... (1 to koh) .....              #     (omit data input if koh=0)
nokq  Qinp                              # Node with given duscharge, given discharge
    ..... (1 to koq) .....              #     (omit data input if koq=0)
noku                                    # Node with saturation point
    ..... (1 to kou) .....              #     (omit data input if kou=0)

変数	説明
npoin, nele, nsec	節点数，要素数，材料特性数
koh, koq, kou	全水頭指定節点数，流量指定節点数，浸潤面境界節点数
idan	解析断面（0：鉛直断面，1：水平断面）
Ak0	要素の飽和透水係数
alpha, em	要素の不飽和透水特性
node1, node2, node3, isec	節点1，節点2，節点3，材料特性番号
x, y, hvec0	節点x座標，節点y座標，全水頭初期値
nokh, Hinp	全水頭指定節点番号，節点全水頭値
nokq, Qinp	流量指定節点番号，節点流量値
noku	浸潤面境界節点番号

（コメント）

変数 idan は，圧力水頭計算時に，鉛直断面では（全水頭マイナス位置水頭），水平断面では（イコール全水頭）で計算する必要があるため，この識別を行うために導入している．
要素の不飽和透水特性 alpha および em は，飽和解析を行いたい場合は，負値（例えば alpha=-1, em=-1）を入力する．これにより，全要素が飽和状態であることを明確にしている．
節点全水頭初期値 hvec0 は，全節点同一の値を入れて構わない．節点全水頭指定境界では，プログラムの中で指定値に入れ替えている．

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

npoin  nele  nsec   koh   koq   kou  idan
    (Each value of above)
sec  Ak0  alpha  em
    Ak0   : Permeability coefficient (saturated)
    alpha : Un-saturated parameter
    em    : Un-saturated parameter
    ..... (1 to nsec) .....
node  x  z  hvec  qvec  koh  koq  kou
    node : Node number
    x    : x-coordinate
    z    : z-coordinate
    hvec : Initial total head
    qvec : Initial discharge
    koh  : Index of node with given total head (0: free, 1: fixed)
    koq  : Index of node with given discharge  (0: free, 1: fixed)
    kou  : Index of node with saturation point (0: free, 1: fixed)
    ..... (1 to npoin) .....
node Hinp
    node : Node number
    Hinp : Given total head
    ..... (1 to koh) .....
node  Qinp
    node : Node number
    Qinp : Given discharge
    ..... (1 to koq) .....
elem  i  j  k  sec
    elem : Element number
    i    : Node number of start point
    j    : Node number of second point
    k    : Node number of third point
    sec  : Material number
    ..... (1 to nele) .....
node  hvec  pvec  qvec  koh  koq  kou
    node : Node number
    hvec : Total head
    pvec : Pressure head (Total head minus elevation of the node)
    qvec : Discharge (+: inflow, -: outflow)
    koh  : Index of node with given total head (0: free, 1: fixed)
    koq  : Index of node with given discharge  (0: free, 1: fixed)
    kou  : Index of node with saturation point (0: free, 1: fixed)
    ..... (1 to npoin) .....
elem  vx  vz  vm  kr
    elem : Element number
    vx   : Flow velocity in x-direction
    vz   : Flow velocity in z-direction
    vm   : Composed flow velocity
    kr   : Parameter of saturation (1: saturated, less than 1: unsaturated)
    ..... (1 to nele) .....
Total inflow = (total inflow discharge: positive sign)
Total outflow= (total outflow discharge: negative sign)
Max.velocity in all area     =  (maximum velocity in all area)
Max.velocity in inflow area  =  (maximum velocity in inflow area)
Max.velocity in outflow area =  (maximum velocity in outflow area)
iii=(number of itaration)
icount=(converged number of degrees of freedom)
kop=(number of pressured nodes)
n=(total degrees of freedom)  time=(calculation time)

変数	説明
hvec, pvec, qvec	全水頭，圧力水頭，節点流量
vx, vzm vm	要素 x 方向流速，要素 y 方向流速，合成流速
kr	飽和度パラメータ（1：飽和，1未満：不飽和）

プログラミング

プログラム本体を，「プログラムの構成」に沿って以下に示します．

1. モジュール読み込み

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

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

2. 入力データ読み込み (def inpdata_s3)

(1) 基本数値読み込み
(2) 配列宣言
(3) 材料特性データ読み込み
(4) 要素構成節点番号および材料番号読み込み
(5) 節点座標および節点全水頭初期値読み込み
(6) 境界条件読み込み

# data input
def inpdata_s3(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
    koh   =int(text[3]) # Number of nodes with given total head
    koq   =int(text[4]) # Number of nodes with given discharge
    kou   =int(text[5]) # Number of nodes with saturation point
    idan  =int(text[6]) # section (0: vertical, 1: horizontal)
    # array declaration
    ae    =np.zeros((3,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
    hvec  =np.zeros(nfree*npoin,dtype=np.float64)   # Total head vector
    nokh  =np.zeros(koh,dtype=np.int)               # nodes with given total head
    hinp  =np.zeros(koh,dtype=np.float64)           # Total head at specified nodes
    nokq  =np.zeros(koq,dtype=np.int)               # nodes with given discharge
    qinp  =np.zeros(koq,dtype=np.float64)           # Discharge at specified nodes (+: inflow, -: outflow)
    noku  =np.zeros(kou,dtype=np.int)               # nodes with saturation point
    # section characteristics
    for i in range(0,nsec):
        text=f.readline()
        text=text.strip()
        text=text.split()
        ae[0,i]=float(text[0]) #ak0   : Permeability coefficient (saturated)
        ae[1,i]=float(text[1]) #alpha : Unsaturated property
        ae[2,i]=float(text[2]) #em    : Unsaturated property
    # 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]) #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
        hvec[i]=float(text[2]) # Initial total head
    # boundary conditions
    if 0<koh:
        for i in range(0,koh):
            text=f.readline()
            text=text.strip()
            text=text.split()
            nokh[i]=int(text[0])
            hinp[i]=float(text[1])
    if 0<koq:
        for i in range(0,koq):
            text=f.readline()
            text=text.strip()
            text=text.split()
            nokq[i]=int(text[0])
            qinp[i]=float(text[1])
    if 0<kou:
        for i in range(0,kou):
            text=f.readline()
            text=text.strip()
            text=text.split()
            noku[i]=int(text[0])
    f.close()
    return npoin,nele,nsec,koh,koq,kou,idan,ae,node,x,hvec,nokh,hinp,nokq,qinp,noku

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

# print out of input data
def prinp_s3(fnameW,npoin,nele,nsec,koh,koq,kou,idan,ae,node,x,nokh,nokq,noku,hinp,qinp,hvec,qvec):
    fout=open(fnameW,'w')
    # print out of input data
    print('{0:>5s} {1:>5s} {2:>5s} {3:>5s} {4:>5s} {5:>5s} {6:>5s}'
    .format('npoin','nele','nsec','koh','koq','kou','idan'),file=fout)
    print('{0:5d} {1:5d} {2:5d} {3:5d} {4:5d} {5:5d} {6:5d}'
    .format(npoin,nele,nsec,koh,koq,kou,idan),file=fout)
    print('{0:>5s} {1:>15s} {2:>15s} {3:>15s}'
    .format('sec','ak0','alpha','em'),file=fout)
    for i in range(0,nsec):
        print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e}'
        .format(i+1,ae[0,i],ae[1,i],ae[2,i]),file=fout)
    print('{0:>5s} {1:>15s} {2:>15s} {3:>15s} {4:>15s} {5:>5s} {6:>5s} {7:>5s}'
    .format('node','x','z','hvec','qvec','koh','koq','kou'),file=fout)
    for i in range(0,npoin):
        lp=i+1
        kh=0; kq=0; ku=0
        if 0<koh:
            for k in range(0,koh):
                if i==nokh[k]-1: kh=1
        if 0<koq:
            for k in range(0,koq):
                if i==nokq[k]-1: kq=1
        if 0<kou:
            for k in range(0,kou):
                if i==noku[k]-1: ku=1
        print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e} {4:15.7e} {5:5d} {6:5d} {7:5d}'
        .format(lp,x[0,i],x[1,i],hvec[i],qvec[i],kh,kq,ku),file=fout)
    if 0<koh:
        print('{0:>5s} {1:>5s}'.format('node','Hinp'),file=fout)
        for k in range(0,koh):
            print('{0:5d} {1:15.7e}'
            .format(nokh[k],hinp[k]),file=fout)
    if 0<koq:
        print('{0:>5s} {1:>5s}'.format('node','qinp'),file=fout)
        for k in range(0,koq):
            print('{0:5d} {1:15.7e}'
            .format(nokq[k],qinp[k]),file=fout)
    print('{0:>5s} {1:>5s} {2:>5s} {3:>5s} {4:>5s}'
    .format('elem','i','j','k','sec'),file=fout)
    for ne in range(0,nele):
        print('{0:5d} {1:5d} {2:5d} {3:5d} {4:5d}'
        .format(ne+1,node[0,ne],node[1,ne],node[2,ne],node[3,ne]),file=fout)
    fout.close()

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

# print out of calculation result
def prout_s3(fnameW,npoin,nele,nsec,koh,koq,kou,idan,ae,node,x,nokh,nokq,noku,hvec,qvec,evx,evz,evm,ekr):
    fout=open(fnameW,'a')
    # total head, pressure head, discharge
    print('{0:>5s} {1:>15s} {2:>15s} {3:>15s} {4:>5s} {5:>5s} {6:>5s}'.format('node','hvec','pvec','qvec','koh','koq','kou'),file=fout)
    for i in range(0,npoin):
        lp=i+1
        kh=0; kq=0; ku=0
        if 0<koh:
            for k in range(0,koh):
                if i==nokh[k]-1: kh=1
        if 0<koq:
            for k in range(0,koq):
                if i==nokq[k]-1: kq=1
        if 0<kou:
            for k in range(0,kou):
                if i==noku[k]-1: ku=1
        if idan==0:
            print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e} {4:5d} {5:5d} {6:5d}'.format(lp,hvec[i],hvec[i]-x[1,i],qvec[i],kh,kq,ku),file=fout)
        if idan==1:
            print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e} {4:5d} {5:5d} {6:5d}'.format(lp,hvec[i],hvec[i],qvec[i],kh,kq,ku),file=fout)
    # velocity
    print('{0:>5s} {1:>15s} {2:>15s} {3:>15s} {4:>15s}'
    .format('elem','vx','vz','vm','kr'),file=fout)
    for ne in range(0,nele):
        print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e} {4:15.7e}'
        .format(ne+1,evx[ne],evz[ne],evm[ne],ekr[ne]),file=fout)
    fout.close()

5. 関数：不飽和透水特性計算 (def cal_kre_s3)

まず，要素は飽和していると仮定し，比透水係数を１にセット（akre=1.0）する．
次に，サクション（負圧）が生じていれば要素は不飽和であると判断し，飽和度 $S_e$（変数：se）および比透水係数 $K_r$（変数：akre）を算定している．
戻り値は，比透水係数 $K_r$（変数：akre）のみである．

\begin{align}
&S_e=\left\{\cfrac{1}{1+(\alpha\cdot H_s)^n}\right\}^m \qquad\qquad n=\cfrac{1}{1-m} \qquad (0<m<1)\\
&K_r=(S_e)^{0.5}\cdot\left\{1-\left(1-S_e{}^{1/m}\right)^m\right\}^2 \qquad (0<S_e, K_R\leq1)
\end{align}

# calculation of degree of saturation
def cal_kre_s3(pp,alpha,em):
    sucp=-pp
    akre=1.0
    if 0.0<sucp:
        en=1.0/(1.0-em)
        se=(1.0/(1.0+(alpha*sucp)**en))**em
        akre=np.sqrt(se)*(1.0-(1.0-se**(1.0/em))**em)**2
    return akre

6. 関数：要素透水マトリックス作成 (def sm_s3)

要素の飽和ー不飽和透水係数 $k$ は，飽和透水係数 $k_0$ と不飽和時の比透水係数 $k_r$ の積となる．

\begin{equation}
k=k_0\cdot k_r
\end{equation}

要素の比透水係数は，各節点での圧力水頭から計算した比透水係数の平均値としている．

\begin{equation}
[\boldsymbol{k}]=
\frac{k}{4\Delta}
\begin{bmatrix}
b_i b_i + c_i c_i & b_j b_i + c_j c_i & b_k b_i + c_k c_i \\
b_i b_j + c_i c_j & b_j b_j + c_j c_j & b_k b_j + c_k c_j \\
b_i b_k + c_i c_k & b_j b_k + c_j c_k & b_k b_k + c_k c_k
\end{bmatrix}
\end{equation}

# element permeability matrix
def sm_s3(ak0,alpha,em,h1,h2,h3,x1,y1,x2,y2,x3,y3,idan):
    sm=np.zeros([3,3],dtype=np.float64)
    # calculation of pressure head (total head - position head)
    p1=h1-y1 # pressure head at node_i
    p2=h2-y2 # pressure head at node_j
    p3=h3-y3 # pressure head at node_k
    delta=0.5*((x3-x2)*y1+(x1-x3)*y2+(x2-x1)*y3)
    b1=y2-y3; c1=x3-x2
    b2=y3-y1; c2=x1-x3
    b3=y1-y2; c3=x2-x1
    sm[0,0]=b1*b1+c1*c1; sm[0,1]=b2*b1+c2*c1; sm[0,2]=b3*b1+c3*c1
    sm[1,0]=b1*b2+c1*c2; sm[1,1]=b2*b2+c2*c2; sm[1,2]=b3*b2+c3*c2
    sm[2,0]=b1*b3+c1*c3; sm[2,1]=b2*b3+c2*c3; sm[2,2]=b3*b3+c3*c3
    if idan==0:
        akr1=cal_kre_s3(p1,alpha,em)
        akr2=cal_kre_s3(p2,alpha,em)
        akr3=cal_kre_s3(p3,alpha,em)
        akr=(akr1+akr2+akr3)/3
    if idan==1:
        akr=1.0
    sm=ak0*akr/4/delta*sm
    return sm

7. 関数：要素流速計算 (def calv_s3)

\begin{align}
&v_x=-k\frac{\partial \phi}{\partial x} =-k\frac{\partial [\boldsymbol{N}]}{\partial x}\{\boldsymbol{\phi}\}=-\frac{k}{2\Delta}(b_i \phi_i+b_j \phi_j+b_k \phi_k) \\
&v_y=-k\frac{\partial \phi}{\partial y} =-k\frac{\partial [\boldsymbol{N}]}{\partial y}\{\boldsymbol{\phi}\}=-\frac{k}{2\Delta}(c_i \phi_i+c_j \phi_j+c_k \phi_k)
\end{align}

# velocity calculation
def calv_s3(ak0,alpha,em,h1,h2,h3,x1,y1,x2,y2,x3,y3,idan):
    # calculation of pressure head (total head - position head)
    p1=h1-y1 # pressure head at node_i
    p2=h2-y2 # pressure head at node_j
    p3=h3-y3 # pressure head at node_k
    delta=0.5*((x3-x2)*y1+(x1-x3)*y2+(x2-x1)*y3)
    b1=y2-y3; c1=x3-x2
    b2=y3-y1; c2=x1-x3
    b3=y1-y2; c3=x2-x1
    vx=0.0
    vz=0.0
    akr=0.0
    if idan==0:
        akr1=cal_kre_s3(p1,alpha,em)
        akr2=cal_kre_s3(p2,alpha,em)
        akr3=cal_kre_s3(p3,alpha,em)
        akr=(akr1+akr2+akr3)/3
    if idan==1:
        akr=1.0
    vx=b1*h1+b2*h2+b3*h3
    vz=c1*h1+c2*h2+c3*h3
    vx=-ak0*akr*vx/2/delta
    vz=-ak0*akr*vz/2/delta
    vm=np.sqrt(vx**2+vz**2)
    return akr,vx,vz,vm

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

(1) 入出力ファイル名取得

# Main routine
def main_seep3():
    start=time.time()
    args = sys.argv
    fnameR=args[1] # input data file
    fnameW=args[2] # output data file

(2) 基本数値定義

収束計算の上限回数 itmax，１要素節点数 nod，１節点自由度数 nfree を指定する．

    itmax=100      # maximum number of iteration
    nod=3          # number of nodes per element
    nfree=1        # degree of freedom per node

(3) 入力データ読み込み

関数 inpdata_s3 により入力データを読み込む．これにより，総節点数，要素数などの基本数値が定まるため，必要な配列の初期化を行う．

    # data input
    npoin,nele,nsec,koh,koq,kou,idan,ae,node,x,hvec,nokh,hinp,nokq,qinp,noku=inpdata_s3(fnameR,nod,nfree)
    # array declaration
    qvec  =np.zeros(nfree*npoin,dtype=np.float64)   # Discharge vector
    nokp  =np.zeros(kou,dtype=np.int)               # nodes with saturation point
    ir    =np.zeros(nod*nfree,dtype=np.int)         # Work vector for matrix assembly

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

全水頭ベクトル hvec および流量ベクトル qvec は初期化により，全要素ゼロとなっているが，入力データ書き出し前に，全水頭指定節点および流量指定節点の全水頭および流量値を指定値（hinp, qinp）に更新している．

    # setting of given head and discharge
    if 0<koh:
        for k in range(0,koh):
            hvec[nokh[k]-1]=hinp[k] # set given total head
    if 0<koq:
        for k in range(0,koq):
            qvec[nokq[k]-1]=qinp[k] # set given discharge
    # print out of input data
    prinp_s3(fnameW,npoin,nele,nsec,koh,koq,kou,idan,ae,node,x,nokh,nokq,noku,hinp,qinp,hvec,qvec)

(5) 反復計算過程

反復計算に入る前に，収束判定に用いる全水頭ベクトル h1vec の初期設定を行っている．

    # iterative calculation
    h1vec=hvec
    for iii in range(0,itmax):

a. 透水性マトリックス組み立て

仮定した全水頭ベクトル h1vec を用い，要素透水性マトリックスを作成すし，全体透水性マトリックス gk に組み込んでいく．

        # Assembly of permeability matrix, [k] is calculated using {h1}
        gk=np.zeros((nfree*npoin,nfree*npoin),dtype=np.float64)  # Global stiffness matrix
        for ne in range(0,nele):
            i=node[0,ne]-1
            j=node[1,ne]-1
            k=node[2,ne]-1
            m=node[3,ne]-1
            x1=x[0,i]; x2=x[0,j]; x3=x[0,k]
            y1=x[1,i]; y2=x[1,j]; y3=x[1,k]
            ak0  =ae[0,m]
            alpha=ae[1,m]
            em   =ae[2,m]
            h1=hvec[i]
            h2=hvec[j]
            h3=hvec[k]
            sm=sm_s3(ak0,alpha,em,h1,h2,h3,x1,y1,x2,y2,x3,y3,idan)
            ir[2]=k; ir[1]=j; ir[0]=i
            for i in range(0,nod*nfree):
                it=ir[i]
                for j in range(0,nod*nfree):
                    jt=ir[j]
                    gk[it,jt]=gk[it,jt]+sm[i,j]

b. 節点流量計算

組み立てた全体透水性マトリックス gk と仮定された全水頭ベクトル h1vec から，節点流量ベクトル qvec を算定する．

        # calculation of nodal discharge {q}=[k]{h1}
        qvec=np.dot(gk,h1vec)

c. 作業用流量ベクトル更新

連立一次方程式作成に先立ち，流量ベクトルを確からしい値に修正するため，以下の操作を行う．

作業用流量ベクトル wvec を作成・初期化（全要素ゼロセット）する．
全水頭指定境界での流量は，先に計算した節点流量ベクトルの値に更新する
流量指定境界での流量は，入力データで指定した指定値に更新する
浸潤面境界での流量は，先に計算した節点流量ベクトルの値に更新する．ただし値が正（流入）の場合はゼロとする．

        #Update of nodal discharge {q}->{q0}
        #(In-flow is not allowed at seepage point. In-flow means positive discharge.)
        #(If discharge has positive value at seepage point, the discharge is set to zero)
        wvec=np.zeros(nfree*npoin,dtype=np.float64)   # Work vector for discharge
        if 0<koh:
            for k in range(0,koh):
                wvec[nokh[k]-1]=qvec[nokh[k]-1]
        if 0<koq:
            for k in range(0,koq):
                wvec[nokq[k]-1]=qinp[k]
        if 0<kou:
            for k in range(0,kou):
                wvec[noku[k]-1]=qvec[noku[k]-1]
                if 0.0<wvec[noku[k]-1]: wvec[noku[k]-1]=0.0
        # Memory of node numbers with positive pressure head in seepage points
        kop=0
        if 0<kou:
            for k in range(0,kou):
                i=noku[k]-1
                if x[1,i]<h1vec[i]: # if h1[i] is greater than z[i] (exclude equal to)
                    kop=kop+1
                    nokp[kop-1]=i+1
        # Amendment of nodal values {q}={q0}-[k]{h(known)}
        # To set total head at the node with given total head
        # To set zero pressure head at the node with positive presure head in seepage points
        #(Total head = position head)

d. 境界条件処理（全水頭既知自由度の処理）

浸潤面境界節点は，大気に接しているため，圧力水頭は少なくともゼロ以下である必要がある．
これを境界条件処理の条件に取り入れ，以下の操作を行う．

浸潤面境界節点のうち，圧力水頭が正，すなわち全水頭が位置水頭より大きい（x[1,i] < h1vec[i]，x[1,i] は節点 i の z 座標）節点の個数（kop）と番号（nokp）を記憶する．
境界条件処理において，全水頭指定節点の流量ベクトルの補正を指定全水頭を用いて行う．（ wvec=wvec-hinp[k]*gk[:,iz] ）
先に記憶した，浸潤面境界節点のうち圧力水頭が正となっている節点において，流量ベクトルの補正を，全水頭を位置水頭に等しいとして行う．（ wvec=wvec-x[1,iz]*gk[:,iz] ）

        # treatment of boundary conditions
        if 0<koh:
            for k in range(0,koh):
                iz=nokh[k]-1
                wvec[iz]=0.0
        if 0<kop:
            for k in range(0,kop):
                iz=nokp[k]-1
                wvec[iz]=0.0
        if 0<koh:
            for k in range(0,koh):
                iz=nokh[k]-1
                wvec=wvec-hinp[k]*gk[:,iz]
                gk[:,iz]=0.0
                gk[iz,iz]=1.0
        if 0<kop:
            for k in range(0,kop):
                iz=nokp[k]-1
                wvec=wvec-x[1,iz]*gk[:,iz]
                gk[:,iz]=0.0
                gk[iz,iz]=1.0

e. 連立一次方程式を解く（全水頭計算）

全体透水性マトリックス gk を，疎行列として圧縮保存した後，spsplve で連立一次方程式を解いている．

        # solution of simultaneous linear equations
        # hvec = np.linalg.solve(gk, wvec)
        gk = csr_matrix(gk)
        hvec = spsolve(gk, wvec, use_umfpack=True)
        del gk,wvec

f. 境界条件処理（既知全水頭の組み込み）

連立一次方程式の解として得られた全水頭に対し，既知全水頭の組み込みとして以下の操作を行う．

全水頭指定節点の全水頭指定値を全水頭ベクトルに組み込む．
先に記憶した，浸潤面境界節点のうち圧力水頭が正となっている節点において，位置水頭を全水頭ベクトルに組み込む．
補正した全水頭ベクトルを h2vec として定義する．

        if 0<koh:
            for k in range(0,koh):
                hvec[nokh[k]-1]=hinp[k]
        if 0<kop:
            for k in range(0,kop):
                hvec[nokp[k]-1]=x[1,nokp[k]-1]
        h2vec=hvec

g. 収束判定

前回の反復過程で得られた全水頭ベクトル h1vec と今回の反復過程で得られた全水頭ベクトル h2vec の差分により収束判定を行う．
前回と今回の全水頭ベクトルの差分が，全節点で $1 \times 10^{-6}$ 未満となった時，収束とみなし，反復計算のループを抜ける．

収束しなければ，全水頭ベクトルを更新（ h1vec = h2vec ）して，反復計算を繰り返す．

        # Judgement of convergence |h2-h1|<1e-6
        icount=0
        for i in range(0,nfree*npoin):
            if np.abs(h2vec[i]-h1vec[i])<1.0e-6: icount=icount+1
        print('nnn={0}  icount={1}  kop={2}'.format(iii,icount,kop))
        if icount==nfree*npoin: break
        # Update of total head {h1}={h2}
        h1vec=h2vec

(6) 要素流速ベクトル計算

このプログラムの前提では，解析領域の全流入量と全流出量が等しくなる必要があるため，これを確認するため，全流入量 qti と全流出量 qto を求めている．
その後の for ne in range(0,nele): のループ内が，要素流速ベクトルの計算となっている．出力の一部とするため，流入領域と流筒領域での流速の最大値検索も同時に行っている．

    # Calculation of total in-flow and out-flow
    qti=0.0  #Total in-flow
    qto=0.0  #Total out-flow
    for i in range(0,nfree*npoin):
        if 0.0<qvec[i]: qti=qti+qvec[i]
        if qvec[i]<0.0: qto=qto+qvec[i]
    # Calculation of element velocity & to search max.velocity
    vmmax=0.0  #Max.velocity
    vimax=0.0  #Max.velocity on in-flow area
    vomax=0.0  #max.velocity on out-flow area
    kmmax=0    #Element number with max.velocity
    kimax=0    #Element number with max.velocity on in-flow area
    komax=0    #Element number with max.velocity on out-flow area
    evx  =np.zeros(nele,dtype=np.float64)
    evz  =np.zeros(nele,dtype=np.float64)
    evm  =np.zeros(nele,dtype=np.float64)
    ekr  =np.zeros(nele,dtype=np.float64)
    for ne in range(0,nele):
        i=node[0,ne]-1
        j=node[1,ne]-1
        k=node[2,ne]-1
        m=node[3,ne]-1
        x1=x[0,i]; x2=x[0,j]; x3=x[0,k]
        y1=x[1,i]; y2=x[1,j]; y3=x[1,k]
        ak0  =ae[0,m]
        alpha=ae[1,m]
        em   =ae[2,m]
        h1=hvec[i]
        h2=hvec[j]
        h3=hvec[k]
        akr,vx,vz,vm=calv_s3(ak0,alpha,em,h1,h2,h3,x1,y1,x2,y2,x3,y3,idan)
        evx[ne]=vx
        evz[ne]=vz
        evm[ne]=vm
        ekr[ne]=akr
        if vmmax<vm:
            kmmax=ne+1
            vmmax=vm
        if 0.0<(qvec[i]+qvec[j]+qvec[k]):
            if vimax<vm:
                kimax=ne+1
                vimax=vm
        if (qvec[i]+qvec[j]+qvec[k])<0.0:
            if vomax<vm:
                komax=ne+1
                vomax=vm

(7) 計算結果書き出し

    # print out of result
    prout_s3(fnameW,npoin,nele,nsec,koh,koq,kou,idan,ae,node,x,nokh,nokq,noku,hvec,qvec,evx,evz,evm,ekr)
    # information
    dtime=time.time()-start
    print('n={0}  time={1:.3f}'.format(nfree*npoin,dtime)+' sec')
    fout=open(fnameW,'a')
    print('Total inflow ={0:15.7e}'.format(qti),file=fout)
    print('Total outflow={0:15.7e}'.format(qto),file=fout)
    print('Max.velocity in all area     ={0:15.7e} (ne={1})'.format(vmmax,kmmax),file=fout)
    print('Max.velocity in inflow area  ={0:15.7e} (ne={1})'.format(vimax,kimax),file=fout)
    print('Max.velocity in outflow area ={0:15.7e} (ne={1})'.format(vomax,komax),file=fout)
    print('iii={0}  icount={1}  kop={2}'.format(iii+1,icount,kop),file=fout)
    print('n={0}  time={1:.3f}'.format(nfree*npoin,dtime)+' sec',file=fout)
    fout.close()

9. メインルーチン実行

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

以　上