Pythonによる２次元飽和ー不飽和定常浸透流解析プログラム（四角形要素）の理論を少し詳しく解説する（改訂版）

• 2017.12.31 投稿
• 2018.08.11 プログラム改良

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

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

はじめに

浸透流解析プログラムについて，三角形要素に続き，四角形要素によるものの解説を作ってみた．
三角形要素用のプログラムとの違いは，入力データやマトリックス組み立てを四角形要素用にしていること，透水性マトリックス作成に数値積分を使っていることである．

解析理論は三角形要素用のものと全く同じ，四角形要素の取扱は２次元応力解析プログラムのものと同じなので，ここでは，「プログラミングのために知っているべきこと」として，四角形要素の２次元応力解析で述べた四角形要素の取扱，三角形要素の浸透流解析で述べた有限要素式の結論のみを記載し，「四角形要素による定式化」で四角形要素を用いた透水性マトリックスの作成，要素平均流速計算方法を述べている．

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

• Pythonによる２次元応力解析プログラム（四角形要素）の理論を少し詳しく解説する（改訂版）
• Pythonによる平面骨組弾性解析プログラムの理論を少し詳しく解説する（改訂版）
• Pythonによる２次元飽和ー不飽和定常浸透流解析プログラム（三角形要素）の理論を少し詳しく解説する（改訂版）
• Pythonによる２次元飽和ー不飽和定常浸透流解析プログラム（四角形要素）の理論を少し詳しく解説する（改訂版）

プログラミングのために知っているべきこと

四角形要素の内挿関数関係

この部分は，四角形要素による２次元応力解析の解説で述べた部分と同じである．

以下，$(a,b)$ は，$a=[−1,1]$，$b=[−1,1]$ の範囲を有する正規化パラメータ座標である．
また，$(x,y)$ は，実モデルの $x$ および $y$ 座標である．

内挿関数の数式表現

\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}

ヤコビ行列と行列式

ヤコビ行列 $[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}

\begin{equation}
\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}

ガウス積分

積分点を 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}

浸透流解析のための知識

この部分は，三角形要素の浸透流解析の解説で述べた部分と同じである．

浸透流解析のための有限要素式

\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}}$ は流速ベクトルの境界長さでの積分であり，流量を示すことは容易に理解できる．

要素流速算定式

$x$ 方向および $y$ 方向の流速をそれぞれ $v_x$，$v_y$ とすれば，ダルシーの法則より

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

不飽和透水特性を示す式（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\leq 1\\
&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}

四角形要素の導入

この部分は，四角形要素による浸透流解析を行うために必要な説明である．

要素の有限要素式

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

\begin{equation}
\begin{Bmatrix} q_i \\ q_j \\ q_k \\ q_l \end{Bmatrix}
=\begin{bmatrix}
k_{11} & k_{12} & k_{13} & k_{14}  \\
k_{21} & k_{22} & k_{23} & k_{24}  \\
k_{31} & k_{32} & k_{33} & k_{34}  \\
k_{41} & k_{42} & k_{43} & k_{44} 
\end{bmatrix}
\begin{Bmatrix} \phi_i \\ \phi_j \\ \phi_k \\ \phi_l \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$の全水頭}\\
&q_l & &\text{：節点$l$の流量} & & \phi_l & &\text{：節点$l$の全水頭}
\end{align}

透水性マトリックス

面積積分の積分点を4点，線積分の積分点を2点とすれば，重みは1となるため，それぞれの行列・ベクトルは積分点での計算値の総和として次のように示される．

\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 \\
=&\int_{-1}^{1}\int_{-1}^{1} \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)\cdot det(J)\cdot da\cdot db \\
=&\sum_{kk=1}^{4} \left[\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)\cdot det(J) \right]_ {kk}
\end{align}

等方均質材料に対しては，$k_x=k_y=k$ として，

\begin{equation}
[\boldsymbol{k}]=
k\cdot \sum_{kk=1}^{4} \left[\left(\frac{\partial [\boldsymbol{N}]^T}{\partial x}\frac{\partial [\boldsymbol{N}]}{\partial x}+\frac{\partial [\boldsymbol{N}]^T}{\partial y}\frac{\partial [\boldsymbol{N}]}{\partial y}\right)\cdot det(J) \right]_ {kk}
\end{equation}

要素内平均流速

要素内任意点の全水頭 $\phi$ は内挿関数 $[\boldsymbol{N}]$ と節点全水頭ベクトル ${\boldsymbol{h}}$ を用いて
ダルシーの法則より要素内任意点の流速 $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}\} \\
&v_y=-k_y\frac{\partial \phi}{\partial y} =-k_y\frac{\partial [\boldsymbol{N}]}{\partial y}\{\boldsymbol{\phi}\}
\end{align}

要素内平均流速は，要素内の4つのGauss積分点の平均とすればよいため，以下の通り求められる．

\begin{align}
&v_x=-\frac{k_x}{4}\sum_{kk=1}^{4}\left(\frac{\partial [\boldsymbol{N}]}{\partial x}\{\boldsymbol{\phi}\}\right)_ {kk} \\
&v_y=-\frac{k_y}{4}\sum_{kk=1}^{4}\left(\frac{\partial [\boldsymbol{N}]}{\partial y}\{\boldsymbol{\phi}\}\right)_ {kk}
\end{align}

節点流量ベクトル

節点流量ベクトル ${\boldsymbol{q}}$ については，等価節点ベクトルとして直接入力すればよい．

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

プログラムの概要

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

扱える境界条件の説明

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

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

三角形一次要素に対し，数値積分をおこなうため，5.および 6. が追加されている．
メインルーチンの内容は，三角形要素用を四角形要素用に書き換えただけで，同じである．

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

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

解析実行コマンド

python3 py_fem_seep4x.py inp.txt out.txt

変数	説明
py_fem_seep4.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  node-4  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, node-4, 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  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  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：計算時間計測用

#=======================================
# 2DFEM seepage Flow 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_s4)

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

# data input
def inpdata_s4(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]) #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
        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_s4)

# print out of input data
def prinp_s4(fnameW,npoin,nele,nsec,koh,koq,kou,idan,ae,x,node,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} {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_s4)

# print out of calculation result
def prout_s4(fnameW,npoin,nele,nsec,koh,koq,kou,idan,ae,x,node,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_s4)

まず，要素は飽和していると仮定し，比透水係数を１にセット（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_s4(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 dmat_s4)

\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}
\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}

# differential of shape function [n]
def dnmat_s4(a,b,x1,y1,x2,y2,x3,y3,x4,y4):
    dnx=np.zeros(4,dtype=np.float64)
    dnz=np.zeros(4,dtype=np.float64)
    #[dn/dadn/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
    #[dn/dx][dn/dy]
    dnx[0]=( J22*dn1a-J12*dn1b)/detJ
    dnx[1]=( J22*dn2a-J12*dn2b)/detJ
    dnx[2]=( J22*dn3a-J12*dn3b)/detJ
    dnx[3]=( J22*dn4a-J12*dn4b)/detJ
    dnz[0]=(-J21*dn1a+J11*dn1b)/detJ
    dnz[1]=(-J21*dn2a+J11*dn2b)/detJ
    dnz[2]=(-J21*dn3a+J11*dn3b)/detJ
    dnz[3]=(-J21*dn4a+J11*dn4b)/detJ
    return dnx,dnz,detJ

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

プログラム上での $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}

# Gauss integration point
def cal_ab_s4(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_s4)

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

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

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

\begin{equation}
[\boldsymbol{k}]=
k\cdot \sum_{kk=1}^{4} \left[\left(\frac{\partial [\boldsymbol{N}]^T}{\partial x}\frac{\partial [\boldsymbol{N}]}{\partial x}+\frac{\partial [\boldsymbol{N}]^T}{\partial y}\frac{\partial [\boldsymbol{N}]}{\partial y}\right)\cdot det(J) \right]_ {kk}
\end{equation}

# element permeability matrix
def sm_s4(ak0,alpha,em,h1,h2,h3,h4,x1,y1,x2,y2,x3,y3,x4,y4,idan):
    sm=np.zeros([4,4],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
    p4=h4-y4 # pressure head at node_l
    for kk in range(0,4):
        a,b=cal_ab_s4(kk)
        dnx,dnz,detJ=dnmat_s4(a,b,x1,y1,x2,y2,x3,y3,x4,y4)
        if idan==0:
            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)
            pp=nm1*p1+nm2*p2+nm3*p3+nm4*p4
            akre=cal_kre_s4(pp,alpha,em)
        if idan==1:
            akre=1.0
        wxz=np.zeros([4,4],dtype=np.float64)
        for i in range(0,4):
            for j in range(0,4):
                wxz[i,j]=wxz[i,j]+dnx[i]*dnx[j]+dnz[i]*dnz[j]
        sm=sm+ak0*akre*wxz*detJ
        del wxz
    return sm

9. 関数：要素流速計算 (def calv_s4)

\begin{align}
&v_x=-\frac{k}{4}\sum_{kk=1}^{4}\left(\frac{\partial [\boldsymbol{N}]}{\partial x}\{\boldsymbol{\phi}\}\right)_ {kk} \\
&v_y=-\frac{k}{4}\sum_{kk=1}^{4}\left(\frac{\partial [\boldsymbol{N}]}{\partial y}\{\boldsymbol{\phi}\}\right)_ {kk}
\end{align}

# velocity calculation
def calv_s4(ak0,alpha,em,h1,h2,h3,h4,x1,y1,x2,y2,x3,y3,x4,y4,idan):
    vx=0.0
    vz=0.0
    akr=0.0
    # 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
    p4=h4-y4 # pressure head at node_l
    for kk in range(0,4):
        a,b=cal_ab_s4(kk)
        dnx,dnz,detJ=dnmat_s4(a,b,x1,y1,x2,y2,x3,y3,x4,y4)
        if idan==0:
            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)
            pp=nm1*p1+nm2*p2+nm3*p3+nm4*p4
            akre=cal_kre_s4(pp,alpha,em)
        if idan==1:
            akre=1.0
        akr=akr+akre
        wvx=dnx[0]*h1+dnx[1]*h2+dnx[2]*h3+dnx[3]*h4
        wvz=dnz[0]*h1+dnz[1]*h2+dnz[2]*h3+dnz[3]*h4
        vx=vx+wvx
        vz=vz+wvz
    akr=0.25*akr
    vx=-0.25*ak0*akr*vx
    vz=-0.25*ak0*akr*vz
    vm=np.sqrt(vx**2+vz**2)
    return akr,vx,vz,vm

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

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

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

(2) 基本数値定義

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

    itmax=100
    nod=4               # 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_s4(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
    prinp_s4(fnameW,npoin,nele,nsec,koh,koq,kou,idan,ae,x,node,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
            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]
            ak0  =ae[0,m]
            alpha=ae[1,m]
            em   =ae[2,m]
            h1=hvec[i]
            h2=hvec[j]
            h3=hvec[k]
            h4=hvec[l]
            sm=sm_s4(ak0,alpha,em,h1,h2,h3,h4,x1,y1,x2,y2,x3,y3,x4,y4,idan)
            ir[3]=l; 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 nodaal 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
        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]
        ak0  =ae[0,m]
        alpha=ae[1,m]
        em   =ae[2,m]
        h1=hvec[i]
        h2=hvec[j]
        h3=hvec[k]
        h4=hvec[l]
        akr,vx,vz,vm=calv_s4(ak0,alpha,em,h1,h2,h3,h4,x1,y1,x2,y2,x3,y3,x4,y4,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]+qvec[l]):
            if vimax<vm:
                kimax=ne+1
                vimax=vm
        if (qvec[i]+qvec[j]+qvec[k]+qvec[l])<0.0:
            if vomax<vm:
                komax=ne+1
                vomax=vm

(7) 計算結果書き出し

    # print out of result
    prout_s4(fnameW,npoin,nele,nsec,koh,koq,kou,idan,ae,x,node,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()

11. メインルーチン実行

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

以　上