Modern FortranでFEMのプログラムを作る（２）

Modern FortranでFEMのプログラムを作る（１）の続きです。
前回解説していなかった部分を追加で説明していきます。

##有限要素方程式に関するオブジェクト
有限要素法では各要素で有限要素方程式を構築します。

K_{ij} u_j = f_{i}

$K_{ij}$および$f_i$は要素の節点数×節点自由度の行列サイズになります。三次元の四面体二次要素で変位を解く場合には10×3=30です。左辺の要素剛性行列$K_{ij}$はBマトリックスとDマトリックスを積分することで計算します。一般に、ガウス積分が用いられます。$f_i$は熱応力など計算する際に、各要素で計算されますが今回のestif_elas中では単純な弾性解析問題なので0として計算します。

\begin{equation*}
{\bf K =\int{B^tDB~}}
\end{equation*}dV

本プログラムでは、有限要素方程式に関するオブジェクト：t_femElemを用意しており基本的な機能は以下の通りです。

• 計算の次元を設定し、要素形状を選択する。
• 要素節点の座標を入力する。
• 積分点の数・位置を決定する。
  • 標準的なガウス積分点の数と位置を選択可能（データベースとして内部的に記述している）。マニュアルインプットも可能（重みもマニュアル設定可能）。
• 各節点や積分点で、スカラー量、ベクトル量、テンソル量を設定する。
• 形状関数やBマトリックスを計算する。
• 要素内でガウス積分を行う。

こういった使い方をします。

• 各積分点にDマトリックスを設定して$\int{B^tDB}dV$を計算する。
• 各節点に変位ベクトルを設定して、積分点上でのひずみテンソルを計算する。
• 各積分点に応力テンソルを設定して$\int{B^t\sigma dV}$を計算する(節点力の計算) 。

**有限要素方程式が理解できれば、ガウス積分など知らなくても使えるオブジェクトの作成を目指しました。**要素内積分に関するサブルーチンをいくつか用意し、熱応力の算出、熱容量行列の計算など、有限要素方程式で記述される多数の計算などが可能です。
また、前回の記事と重複する部分もありますが、プログラムの思想は以下の通りです。OOPだと当たり前だろ、という点も記載しますが、ご容赦ください。

• 変数は全てprivateにする。内部変数へのアクセス・参照は、サブルーチンもしくはファンクションによって実装する（ゲッター・セッター）。
  • 例えば、ガウス点の数はfunction ngausとして実装。プログラム中では、ig=feme%ngaus()として使用。
  • これほどまでに変数への直接アクセスを嫌っているのは、ユーザーによう予期せぬ変数変更を防ぐためです。
• 使用するサブルーチンのみをpublicとして使用できるようにする。
• 内部パラメータは外部から分からないようにする。例えば、各要素の種類で番号付けをしているが、要素の選択などはサブルーチンを介して実施する（confElemType_iTETR4）。
  • こうしておくと、内部的な変数名の変更など、プログラムの修正した場合にも、ユーザー側のプログラムの変更が最小限で済むことが多いです。後でプログラムを修正できると考えると、かなり気楽に作成できます。
• オブジェクト内にデバッグルーチンを用意する。
  • 特にobj_femelem中では、１要素有限要素解析ができるデバッグサブルーチンを用意しており（Paraviewへの出力ファイルも作成）、verificationができるようにしています。
  • **この際に、他の出力に関するオブジェクトを呼ばず、obj_femElem中のコードですべて完結させています。**重複するプログラムが多くなる、というデメリットはありますが、このプログラム単体で動くという点は大きなメリットと考えています。

  type :: t_femElem
    !
    private
    character(50) :: filename = "femelem"
    integer :: ndof,ndim,nmat
    integer :: enon,nump
    integer :: elemtype, intgorder
    !
    integer :: vgtid(6) = [ 3,1,2,3,1,2 ] 
    !
    real    :: vole
    !
    logical :: flag_calcShapf = .true.
    logical :: flag_calcBmat  = .true.
    !
    type(t_fe_point),allocatable :: fepoint(:)
    type(t_fe_node) ,allocatable :: fenode(:)
    !
    contains
    private
    !
    procedure,public :: init => init_femElem
    procedure,public :: fnlz => fnlz_femElem
    !
    ! setter
    procedure,public :: confElemType_iLin2 
    procedure,public :: confElemType_iQUAD4
    procedure,public :: confElemType_iTETR4
    procedure,public :: confElemType_iBRIC8
    procedure,public :: confElemType_iTET10
    procedure,public :: confElemType_iBRI20
    !
    procedure,public :: set_calcShapf
    procedure,public :: set_calcBmat
    !
    procedure,public :: setInfoElem
    procedure,public :: confOrderOfQuadrature
    !
    procedure,public :: addPoint_gauss  ! add gauss points
    procedure,public :: addPoint_xyz    ! add new point, after setInfoElem
    procedure,public :: addPoint_xietzt !
           procedure :: addPoint
           procedure :: gradsh
    !
    ! setter at Point
    procedure,public :: setPoint_scalar !
    procedure,public :: setPoint_voigt  !
    procedure,public :: setPoint_cmat   !
    !
    ! setter at Node
    procedure,public :: setNode_scalar  !
    procedure,public :: setNode_vector  !
    !
    ! getter integrate value
    procedure,public :: integrateScalar
    procedure,public :: integrateVector
    procedure,public :: integrateNtN      !
    procedure,public :: integrateNtScalar => integrateNtVector !
    procedure,public :: integrateNtVector !
    procedure,public :: integrateNtVoigtB !
    procedure,public :: integrateBtCB     !
    procedure,public :: integrateBtVoigt  !
    procedure,public :: integrateBtCVoigt  !
    !
    ! getter interporate value at point
    procedure,public :: NtNode        !
    procedure,public :: NtNode_scalar !
    procedure,public :: BtNode        !
    !
    ! getter
    procedure,public :: ngaus => getNumpoint
    procedure,public :: getNumpoint
    procedure,public :: getPoint_xietzt
    procedure,public :: xyzToXietzt
    !
    procedure,public :: deF => get_deF
    procedure,public :: J => get_Jacobian
    !
    procedure,public :: get_be_xyz
    !
    procedure :: Nm => get_Nmatrix
    procedure :: Bm => get_Bmatrix
    !
    procedure :: init_femElem
    procedure :: fnlz_femElem
    !
    ! for debug
    !
    procedure,public :: debug_LIN2
    procedure,public :: debug_QUAD4
    procedure,public :: debug_BRIC8 !
    !
    procedure,public :: solve_CG
    procedure,public :: setFileName
    procedure,public :: output_vtk
    !
  endtype t_femElem

##obj_femElemを使った有限要素解析
前回紹介したコード中のestif_elas：要素剛性行列作成のコードを簡単な追加解説付きで再掲します。

  !
  ! 要素番号ieにおける要素剛性行列:ske_ij*u_j=sfe_iのサブルーチン
  !
  subroutine estif_elas( ie, ske, sfe )

    use obj_femElem

    implicit none

    type(t_femElem)  :: feme

    integer,intent(in)  :: ie
    real,   intent(out) :: ske(:,:)
    real,   intent(out) :: sfe(:)

    integer  :: ig, mat
    real,allocatable :: enod(:,:)

    mat = mesh%elem(ie)%mat !材料番号に対応

    ske = 0e0
    sfe = 0e0

    allocate( enod(3,mesh%elem(ie)%non ) )
    !
    enod(:,:) = mesh%elem(ie)%ncod(:,:) !要素の節点座標
    !
    call feme%init( ndimela,dofela ) !オブジェクトの初期化、次元、節点自由度の追加
    call feme%confElemType_iTET10() !要素形状の選択。ここでは四面体二次要素
    call feme%setInfoElem( enod ) !節点座標の節点
    call feme%addPoint_gauss() !積分点として、標準のガウス点を追加
    !
    ! Kmatrix
    do ig = 1, feme%ngaus() !ガウス積分点のループ
      call feme%setPoint_cmat( matElas%dmatrix(:,:,mat),ig ) !各積分点にDマトリックスを設定する。
    enddo
    ske = feme%integrateBtCB() !ガウス積分を実施。
    !
    ! vector
    sfe = 0e0
    !
    call feme%fnlz()

    deallocate( enod )

  endsubroutine estif_elas

##レビュー
obj_femElemですが、目標はOOPにより使用しやすいコードを作ることだったのでその点では成功したと考えています。ただし、スピードに関して最適化できているかは怪しいところ、というところが率直な気持ちです。体感としてそれほど遅くなった気がしていないので、ヨシとしました。この辺に詳しい人の情報は随時お待ちしています。

##おわりに
今回は、有限要素解析における主要オブジェクトの概要・その使い方を紹介しました。引き続き、内容については追加していきます。熱伝導解析や熱応力解析まで説明できれば、と思っています。

##おまけ
なお、t_femElemは対称テンソルを前提とした微小変形に関するオブジェクトですが、その拡張として有限変形を考慮可能なオブジェクトを作成する構想もありましたが、結局完成はしていません。。

 type, extends( t_femElem ) :: t_NLfemElem
    contains
    private
    procedure,public :: init => init_NLfemElem
    procedure,public :: fnlz => fnlz_NLfemElem
    !　この辺に有限変形用に拡張したprocedure
 endtype

ですが、extendsの使用方法としては、こういったものが適切かなと考えているので参考までに。