Modern Fortranで書いたFEMのプログラムの紹介
かつて研究室で書いていたFortranのFEMコードをごく簡単に紹介します。
なお、main文のみの紹介ですのでご留意ください（これを引っ張っても、なにも動きません）。
また、一部説明用に省略・修正などしているので、実際のものとは多少異なります。
-
ソースコードのポイントと設計思想
- Fortran2003以降のModern Fortranを利用し、オブジェクト指向プログラミング（OOP）にて作成した。
- メッシュ、境界条件、行列ソルバ、出力などの各オブジェクトをmoduleで用意した。main文中にて組み合わせることで、一連のFEM解析を実行する。
- 各モジュールは互いに独立するように作っている。また、各クラス事にデバッグ・検証問題用のサブルーチンを入れている。
- 設定ファイルは極力使用しない。条件設定は、main文中にべた書き or 設定用のサブルーチンをcallすることで実行。
- グローバル変数は一切使用しない。すべてサブルーチンの引数で引き渡す。
- interface文を用いて、要素剛性行列のsubroutineを引数とするアセンブリングのsubroutineを作成（なお、これはFortran 2003以降の機能ではないですが、非常に便利だったので紹介）。
-
何がやりたかったことか？
- 力学、熱伝導、反応解析やそれらの連成解析など諸々の解析をできるだけ同じサブルーチン使いまわして実現する。
- 行列ソルバのオブジェクトはAx=bを解くだけの機能として設計。変位だろうと、温度だろうと化学種だろうと全部コレで解く。
- プログラムのユーザである研究室の後輩には、自分のやりたい解析を行わせるためにmain文だけをいじらせる。moduleは触らせない（コピペすらさせない）。
- バグや機能追加は管理者の自分が作成（または監修）して、ユーザにモジュールを配布。オブジェクトの使いかたを教える。
- 一方で、自分は自分で書きたいmain文を好き勝手作成する。
- 力学、熱伝導、反応解析やそれらの連成解析など諸々の解析をできるだけ同じサブルーチン使いまわして実現する。
-
やってよかったと感じたこと
- 気軽に・機動的にプログラム（main文）を書くことができて、開発スピードが上がった。
- 気を遣ってOOPで書くとオブジェクト同士の依存性が無くなり、既存のコードの修正・追加が非常にやりやすかった。
- 各クラスごとにデバッグが可能となると問題の切り分けが非常に用意になり、プログラムの追加・修正が格段にやりやすくなった。
- setterを作って条件設定の値を直接引数とせず、subroutineとすることで修正がやりやすかった。
- 例えば、ソルバ選択の部分。各ソルバにパラメータとして番号を割り当ててもよいが、内部の変数が変わったりするとややこしくなるので、、
- 後輩のコードを確認するときに、後輩が作ったmain文の問題か、既存のモジュールの問題か、その切り分けがやりやすかった。
- interface文を利用したことで、コードの汎用性が飛躍的に向上した。
- 材料構成則ごとにアセンブリングのプログラムを作る、なんてことは無くなった。
-
結論（所感）
- OOPによるFEMプログラムは研究室で管理する上で非常に使いやすかった。
- C++などで書いている人は慣れていて当たり前なのかもしれないが、FortranユーザーからするとOOPに触れることがほとんどないので、ここでのステップアップは劇的に感じた。
- アセンブリングのsubrtoutineなど一部はオブジェクト化しておらず、どこまでオブジェクト化していくか、そこの設計思想はまだまだ検討の余地ありだと感じている。
- OOPによるFEMプログラムは研究室で管理する上で非常に使いやすかった。
##線形弾性解析のソースコード
program main_elastic3d
use obj_mesh
use obj_boundary
use obj_matElas
use obj_solver
use mod_commfem, only : setBoundary,countmain
use mod_commfem, only : assemElem, assemNodeFromNeumBoundary
use mod_commPost,only : initMeshdataForPost
use obj_postproc
use obj_femElem
implicit none
type(t_mesh) :: mesh !メッシュオブジェクト
type(t_boundary) :: boundElas !境界条件オブジェクト
type(t_matElas) :: matElas !材料オブジェクト（線形弾性）
type(t_solver) :: solveElas !ソルバオブジェクト
type(t_postproc) :: post !結果出力ファイル作成オブジェクト
integer :: ndimela,dofela, neqela !自由度
real,allocatable :: Uxyz(:) !節点変位
real,allocatable :: vms(:) !応力
real,allocatable :: tns(:,:) !応力テンソル
ndimEla = 3
dofela = 3
mesh%fread = "mesh_elastic.grd"
call mesh%read_mesh_mentat()
call mesh%convertTo2ndOrder()
!
! setMaterial Elast
call matElas%init( mesh%NOM )
call matElas%read()
!
boundElas%fread = "BCtable.dat"
call setBoundary( dofela,mesh,boundElas )
neqela = dofela*mesh%NON
call solveElas%init( neqEla,1e-7,1000 ) !ソルバの収束残差、最大反復数を設定
call solveElas%confBiCG() !利用するソルバの選択
call alloc( .true. )
!
call setPost() !出力条件の設定
!
call countmain( mesh,solveElas,dofela ) !メッシュデータをもとに、非ゼロの行列位置をソルバオブジェクトに認識させる。
!
call assemNodeFromNeumBoundary( dofela,mesh,boundElas,solveElas ) !ノイマン境界条件の設定
call assemElem( dofela,mesh,boundElas,estif_elas,solveElas ) !アセンブリング
!
call solveElas%solve() !行列方程式Ax=bの求解
Uxyz(:) = solveElas%phi()
!
call calcStress() !変位から応力の算出 Uxyz→tns,vms
!
call postmain() !結果の出力
!
contains
!
! 要素番号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 )
enddo
ske = feme%integrateBtCB()
!
! vector
sfe = 0e0
!
call feme%fnlz()
deallocate( enod )
endsubroutine estif_elas
!
! subroutine alloc, setPost, calcStress, postmainはmain文中に記載。
! 詳細は省略。
!
endprogram main_elastic3d
##おわりに
今回の記事では解説しきれてない部分やまとまっていない内容もあるので、内容は追加していければと思います。
Fortranによる理想的なFEMコードとは何か？ の議論はどこかで深めていきたいです。
今度、モダンFortran勉強会とかでネタにできればいいなと思います。
ModernFortran勉強会
以上です
##追記