plantFEMで1次元定常拡散方程式ソルバーを自作して動かす。

plantFEMは、植物のFEM解析に特化したライブラリです。
主に3次元拡散、3次元非線形弾性変形、およびそれらのカップリング計算に対する解析を行えます。
追記：拡散係数マトリクスに誤りがありましたので修正しました。@2021/05/13

plantFEMには4つのレベルがあります。
- Fortranの純粋な拡張としてのstd
- 有限要素法解析に必要なデータ構造と簡単なプリポストを備えたfem
- 個別のソルバーとポストプロセッシング処理を備えたsim
- 種子、葉、花、土といった植物と周辺環境のオブジェクトとそれらへの操作を備えたobj

今回は、plantFEMの機能のうち、stdとfemを使い、1次元定常拡散方程式のソルバーを自作して動かしてみます。
$$
\begin{equation}
- k \frac{d^2 h(x) }{d x^2} = f
\end{equation}
$$

最初に、plantFEMをインストールします。Ubuntuのターミナル上で、

git clone https://github.com/kazulagi/plantFEM.git && cd plantFEM && ./plantfem setup

と打ってエンターキーを押してしばらく待つと、インストールが完了します。

次に、今いるディレクトリでdiffusion1D.f90というファイル（ファイル名は何でも良いです）を開き、編集していきます。
最初に、plantFEMのライブラリをuseして、全機能をアンロックしましょう。
なお、ここでは、plantFEMのライブラリの内、stdとfemしか使いませんので、それらのみをuseしても大丈夫です。
ここで、useはFortranの文法で、pythonのimport等に対応します。
implicit noneはFortranを使うときに唱えねばならない呪文です。

program main
    use plantFEM 
    implicit none

end program main

次に、ライブラリの中から、今回使う2つのアプリを呼び出してきます。
一つは、有限要素法シミュレーションのためのデータベースアプリであるFEMDomain_です。ここでは、domainという名前で実体化します。
もう一つは、連立方程式を解くアプリです。ここでは、LinearSolverという名前で実体化します。
type(アプリ名)::実体化したアプリの名前
のようにして、アプリを実体化し起動できます。
（正確には、FEMDomain_型のインスタンスと、LinearSolver_型のインスタンスを生成しています。）

program main
    use plantFEM 
    implicit none

    type(FEMDomain_)::domain ! 有限要素法シミュレーションのためのデータベースアプリ
    type(LinearSolver_)::LinearSolver ! 連立方程式を解くアプリ

end program main

まず最初に、節点情報、要素コネクティビティ情報、材料情報および境界情報（ここではディリクレ条件のみ）を読み込みます。
FEMDomain_アプリに、読み込み（import）命令を行い、情報の種類（node, element, materialinfo, dirichlet）とファイル名を指定します。

program main
    use plantFEM 
    implicit none

    type(FEMDomain_)::domain ! 有限要素法シミュレーションのためのデータベースアプリ
    type(LinearSolver_)::LinearSolver ! 連立方程式を解くアプリ

    ! FEMデータベースアプリを起動して、データを読み込む

    ! 呼び出し方 >> 領域(domain)　を(%)　読み込む(import)（ファイル名）
    call domain%import(node=.true.,file="./Node.txt")
    call domain%import(element=.true.,file="./Element.txt")
    call domain%import(materialinfo=.true.,file="./Materialinfo.txt")
    call domain%import(Dirichlet=.true.,file="./DBoundary.txt")


end program main

入力ファイルの形式と例を以下に挙げます。

Node.txt

4   1 \\節点数　要素の次元
0.0　　\\節点1番の座標
1.0　　\\節点2番の座標
2.0　　\\節点3番の座標
3.0　　\\節点4番の座標

Element.txt

3 2 \\要素数　1要素あたりの節点数
1 2 \\要素1番を構成する節点番号を列挙
2 3 \\要素2番を構成する節点番号を列挙
3 4 \\要素3番を構成する節点番号を列挙

Materialinfo.txt

3 \\要素数
1 \\要素1番の材料番号
1 \\要素2番の材料番号
1 \\要素3番の材料番号
1 2  \\材料数　材料パラメータの数
1.0 0.0 \\透水係数　領域内の湧き出し/消失速度

DBoundary.txt

2 \\ディリクレ条件を課す節点の数
1 \\1つ目の境界条件の節点ID
4 \\2つ目の境界条件の節点ID
0.5 \\1つ目の境界条件の値
0.2 \\2つ目の境界条件の値

次に、一つ飛ばして連立方程式ソルバーの使用と、結果の表示を記述します。
ここでは、アルゴリズムにBiCGSTABを、データ格納形式にCRS方式を採用します。

program main
    use plantFEM 
    implicit none

    type(FEMDomain_)::domain ! 有限要素法シミュレーションのためのデータベースアプリ
    type(LinearSolver_)::LinearSolver ! 連立方程式を解くアプリ

    ! FEMデータベースアプリを起動して、データを読み込む

    ! 呼び出し方 >> 領域(domain)　を(%)　読み込む(import)（ファイル名）
    call domain%import(node=.true.,file="./Node.txt")
    call domain%import(element=.true.,file="./Element.txt")
    call domain%import(materialinfo=.true.,file="./Materialinfo.txt")
    call domain%import(Dirichlet=.true.,file="./DBoundary.txt")

    ! 一次元定常拡散方程式について、有限要素法で離散化した連立方程式を作り、ソルバーにセットする。
    ! （以下では、このアプリを自作します。）
    call getDiffusionEquation1D(domain,LinearSolver)

    ! ソルバーを使って、連立方程式を解く。
    ! 連立方程式を解くアプリ(LinearSolver) に(%) 連立方程式を解け(solve)、と命令する。（CRSBiCGSTABで）
    call LinearSolver%solve(Solver="BiCGSTAB",CRS=.true.)

    ! 結果を表示する。
    ! 連立方程式を解くアプリ(LinearSolver) の(%)　解(x)を表示する。
    print *, LinearSolver%x(:)

end program main

あとは、domainから情報を呼び出して連立方程式を作り、LinearSolverへ格納するサブルーチンgetDiffusionEquation1Dを定義します。

subroutine getDiffusionEquation1D(domain,LinearSolver)
    use plantFEM
    implicit none

    type(FEMDomain_),intent(inout)::domain
    type(LinearSolver_),intent(inout)::LinearSolver
    integer(int32) :: i,j,n,elemnum,&
        nodeid_1,nodeid_2,elemid,numdirichlet,node_id,dboundid,materialid
    real(real64)   :: k,f,h_1,h_2,Le,x_1,x_2,val
    real(real64)   :: Kmat(2,2),fvec(2)

    elemnum = size(domain%mesh%elemnod,1)
    do i=1, elemnum
        elemid = i
        nodeid_1 = domain%mesh%elemnod(  elemid,1)
        nodeid_2 = domain%mesh%elemnod(  elemid,2)
        materialid = domain%mesh%elemmat(elemid)

        x_1 = domain%mesh%nodcoord(nodeid_1,1)
        x_2 = domain%mesh%nodcoord(nodeid_2,1)

        Le = abs(x_2 - x_1)
        k  = domain%materialprop%matpara(materialid,1)
        f  = domain%materialprop%matpara(materialid,2)

        Kmat(1,1)=  k/Le
        Kmat(1,2)= -k/Le
        Kmat(2,1)= -k/Le
        Kmat(2,2)=  k/Le

        fvec(1) = 0.50d0*f*Le
        fvec(2) = 0.50d0*f*Le

        call LinearSolver%set(nodeid_1,nodeid_1, entryvalue= Kmat(1,1) )
        call LinearSolver%set(nodeid_1,nodeid_2, entryvalue= Kmat(1,2) )
        call LinearSolver%set(nodeid_2,nodeid_1, entryvalue= Kmat(2,1) )
        call LinearSolver%set(nodeid_2,nodeid_2, entryvalue= Kmat(2,2) )

        call LinearSolver%set(nodeid_1, entryvalue=fvec(1) )
        call LinearSolver%set(nodeid_2, entryvalue=fvec(2) )

    enddo

    ! ディリクレ境界条件を導入する。
    numdirichlet = size(domain%Boundary%DboundNodID,1)

    do i=1,numdirichlet
        dboundid = i
        node_id = domain%Boundary%DboundNodID(dboundid,1)
        val     = domain%Boundary%DboundVal(dboundid,1)
        call LinearSolver%fix(node_id, entryvalue=val)
    enddo
end subroutine

以下に、FEMDomainのデータ構造と利用方法を示します。

まず、domain%meshは、メッシュ情報を格納しています。
- domain%mesh%elemnod(要素番号、要素内節点番号)で節点番号を返します。
- domain%mesh%elemmat(要素番号)で要素番号を返します。
- domain%mesh%nodcoord(節点番号,次元)で、節点座標を返します。

Leでは、要素長を得ています。
domain%materialprop%matpara(材料番号,材料パラメータ番号(ここでは1:透水係数、2:湧き出し/消失)　)で、材料パラメータを返します

要素マトリクスを2×2行列のKmatとして、その(1,1)成分から(2,2)成分まで値を入れます。要素は2節点1次要素です。
右辺ベクトルを2列のベクトルfvecとして、その(1)成分と(2)成分にそれぞれ値を入れます。
そして、連立方程式ソルバに格納していきます。

全体マトリクス（全体の連立方程式の係数行列）について、
call LinearSolver%set(行番号に対応する節点番号,列番号に対応する節点番号, entryvalue= 成分(連立方程式の係数) )として、ソルバに全体マトリクスの情報を格納していきます。（中ではCRS形式で保存されています。）

全体ベクトル（全体の連立方程式の右辺定数項）について、
call LinearSolver%set(行番号に対応する節点番号, entryvalue= 成分(右辺定数項の値) )として、ソルバに全体ベクトルの情報を格納していきます。

これを、全要素に渡って行ったあと、ディリクレ境界条件を導入します。
境界条件については、domain%Boundaryに格納されています。

• size(domain%Boundary%DboundNodID,1)は、ディリクレ境界条件を導入する節点の数を返します。
• domain%Boundary%DboundNodID(ディリクレ境界条件番号,次元)は、該当する節点番号を返します。
• domain%Boundary%DboundVal(ディリクレ境界条件番号,次元)は、該当する境界値を返します。

そして、それぞれの境界値を固定します。
call LinearSolver%fix(節点番号, entryvalue=ディリクレ条件として固定する値)

以上になります。

これらをまとめて、一つのスクリプトが出来上がります。

diffusion1dD.f90

subroutine getDiffusionEquation1D(domain,LinearSolver)
    use plantFEM
    implicit none

    type(FEMDomain_),intent(inout)::domain
    type(LinearSolver_)::LinearSolver
    integer(int32) :: i,j,n,elemnum,&
        nodeid_1,nodeid_2,elemid,numdirichlet,node_id,dboundid,materialid
    real(real64)   :: k,f,h_1,h_2,Le,x_1,x_2,val
    real(real64)   :: Kmat(2,2),fvec(2)

    elemnum = size(domain%mesh%elemnod,1)
    do i=1, elemnum
        elemid = i
        nodeid_1 = domain%mesh%elemnod(  elemid,1)
        nodeid_2 = domain%mesh%elemnod(  elemid,2)
        materialid = domain%mesh%elemmat(elemid)

        x_1 = domain%mesh%nodcoord(nodeid_1,1)
        x_2 = domain%mesh%nodcoord(nodeid_2,1)

        Le = abs(x_2 - x_1)
        k  = domain%materialprop%matpara(materialid,1)
        f  = domain%materialprop%matpara(materialid,2)

        Kmat(1,1)=  k/Le
        Kmat(1,2)= -k/Le
        Kmat(2,1)= -k/Le
        Kmat(2,2)=  k/Le

        fvec(1) = 0.50d0*f*Le
        fvec(2) = 0.50d0*f*Le

        call LinearSolver%set(nodeid_1,nodeid_1, entryvalue= Kmat(1,1) )
        call LinearSolver%set(nodeid_1,nodeid_2, entryvalue= Kmat(1,2) )
        call LinearSolver%set(nodeid_2,nodeid_1, entryvalue= Kmat(2,1) )
        call LinearSolver%set(nodeid_2,nodeid_2, entryvalue= Kmat(2,2) )

        call LinearSolver%set(nodeid_1, entryvalue=fvec(1) )
        call LinearSolver%set(nodeid_2, entryvalue=fvec(2) )

    enddo

    ! ディリクレ境界条件を導入する。
    numdirichlet = size(domain%Boundary%DboundNodID,1)

    do i=1,numdirichlet
        dboundid = i
        node_id = domain%Boundary%DboundNodID(dboundid,1)
        val     = domain%Boundary%DboundVal(dboundid,1)
        call LinearSolver%fix(node_id, entryvalue=val)
    enddo
end subroutine


program main
    use plantFEM 
    implicit none

    type(FEMDomain_)::domain ! 有限要素法シミュレーションのためのデータベースアプリ
    type(LinearSolver_)::LinearSolver ! 連立方程式を解くアプリ

    ! FEMデータベースアプリを起動して、データを読み込む

    ! 呼び出し方 >> 領域(domain)　を(%)　読み込む(import)（ファイル名）
    call domain%import(node=.true.,file="./Node.txt")
    call domain%import(element=.true.,file="./Element.txt")
    call domain%import(materialinfo=.true.,file="./Materialinfo.txt")
    call domain%import(Dirichlet=.true.,file="./DBoundary.txt")

    ! 一次元定常拡散方程式について、有限要素法で離散化した連立方程式を作り、ソルバーにセットする。
    ! （以下では、このアプリを自作します。）
    call getDiffusionEquation1D(domain,LinearSolver)

    ! ソルバーを使って、連立方程式を解く。
    ! 連立方程式を解くアプリ(LinearSolver) に(%) 連立方程式を解け(solve)、と命令する。（CRSBiCGSTABで）
    call LinearSolver%solve(Solver="BiCGSTAB",CRS=.true.)

    ! 結果を表示する。
    ! 連立方程式を解くアプリ(LinearSolver) の(%)　解(x)を表示する。
    print *, LinearSolver%x(:)

end program main

先述の入力ファイル群を作成してから、スクリプトをplantFEMで実行します。

./plantfem diffusion1d.f90

すると、結果が得られます。

>>> diffusion1D.f90
  0.50000000000000000       0.40000000000000002       0.29999999999999999       0.20000000000000001     

今回は、ノイマン条件を省きましたが、同様にして、import時のオプションで、neumann=.true.としてノイマン条件を読み込み、上記サブルーチンに少し加筆することで導入できます。