初めに

初投稿です．拙い部分しかないと思いますが，まだまだ初学者ですゆえどうかご了承ください．

この記事はfortranでスパコンなどの大規模並列計算機を使うシミュレーションプログラムを作成するためのヒントとなればいいなあという記事です．ですが，まず初めにこれからシミュレーションプログラムを作成しようと思われる方はこの記事を読む前にこれとこれを買って読んでください．読むべきです．話はこれからです．シミュレーションにおける，まさにバイブルとも呼べる本です（個人の感想です）．私は本当にこの本と同じ時代に生まれたことに感謝しています．シミュレーションでなにをすべきか，全てがここに書いてあります．

ではなぜいまさら私ごときがこんな記事を書くかというと，確かに本の中には「これはC++でもfortranでも同じである」とは書いてあるのですが，中身はC++で書かれてるんですよねー．しかしオブジェクト指向に関するfortranの記事は少ないという…．

というわけで，まずは内容で一番最初に躓いた，構造体とポインタについて書いてみようと思いました．いまのところ以下のやり方で特に不具合はありませんが，もしかしたら問題があるかもしれません．詳しい方は是非コメント欄に指摘ください．全力で勉強させて頂きます．

データモデルの表記法

本の中では有限要素法を例に，メモリレイアウトについて言及されています．以下の図のようなデータ構造を保持したいとします．

米印はポインタであることを表しています．私は有限要素法に明るくないのですが，それでも全てのノードでループを回したいとき，エレメントごとにループを回したいときがそれぞれあると思うのでこのメモリ配置がベストなんだと思います．

fortranのコード

あとは上図を実現するコードを書きましょう．

Node.f90

! クラス名にアンダーバーをつけてモジュール名を定義． 
module Node_
! モジュール名＝クラス名の場合，インデントを下げない方が見やすそう．（エディタが許すかは分かりませんが・・・） 
implicit none
public :: Node

! ここからいわゆるヘッダーに相当する部分
type Node
	private
	real(8) :: x(2), u(2)
	
contains
	procedure :: init => Node_init
end type

! ここからいわゆるソースに相当する部分
contains

subroutine Node_init(this, x, u)
	class(Node) :: this
	real(8), intent(in) :: x(2), u(2)
	
	this%x = x
	this%u = u
end subroutine

end module

他言語で言われているように，１モジュール１ソースコードが良さそうです．また，fortranは言語仕様上１モジュールに複数クラスを定義できますが，管理が破綻しそうなので１モジュール１クラスが良いと思われます．

Element.f90

! 同様にモジュール（クラス）を定義していきます． 
module Element_
use Node_
implicit none
public :: Element
	
type Element
	private
	real(8) :: p, m(4, 4)
	type(Node), pointer :: node1, node2, node3, node4
	
contains
	procedure :: init => Element_init
	procedure :: link => Element_link
end type

contains

subroutine Element_init(this, p, m)
	class(Element) :: this
	real(8),intent(in) :: p, m(4, 4)
	
	this%p = p
	this%m = m
end subroutine

subroutine Element_link(this, node1, node2, node3, node4)
	class(Element) :: this
	type(Node),target  :: node1, node2, node3, node4
	
	this%node1 => node1
	this%node2 => node2
	this%node3 => node3
	this%node4 => node4
end subroutine

end module

値の初期化（この例ではsampleという意味のない値を用いていますが，実際にはちゃんと意味のある値を入れましょう）と，構造体の互いの関係をそれぞれサブルーチンinitとlinkで定義します．

Area.f90

! 同様に(ry
module Area_
use Node_
use Element_
implicit none
	
type Area
	integer :: num_nodes, num_elems
	type(Node),    allocatable :: nodes(:)
	type(Element), allocatable :: elems(:)
contains
	procedure :: init => Area_init
	procedure :: link => Area_link
end type

contains

subroutine Area_init(this, num_nodes, num_elems)
	class(Area) :: this
	integer, intent(in) :: num_nodes, num_elems
	integer :: i
	real(8) :: sample1 = 1d0, sample2(2) = 2d0, sample44(4, 4) = 4d0
	
	this%num_nodes = num_nodes
	this%num_elems = num_elems
	
	allocate(this%nodes(num_nodes))
	do i = 1, this%num_nodes
		call this%nodes(i)%init(sample2, sample2)
	end do
	
	allocate(this%elems(num_elems))
	do i = 1, this%num_elems
		call this%elems(i)%init(sample1, sample44)
	end do
end subroutine
	
subroutine Area_link(this)
	class(Area) :: this
	integer     :: i, j

	do i = 1, this%num_elems
		j = 3 * (i - 1)
		call this%elems(i)%link(this%nodes(j+1), this%nodes(j+2), this%nodes(j+3), this%nodes(j+4))
	end do
end subroutine

end module

これで図に示す構造体の関係はできました．実際の使い方は以下の様になると思います．

main.f90

program main
	use Area_
	implicit none
	
	integer :: num_nodes = 7
	integer :: num_elems = 2
	
	type(Area) :: this_area
	
	call this_area%init(num_nodes, num_elems)
	call this_area%link()
		
end program

これで出来たと思います．どうでしょうか．

更新

2020/02/29　クラスごとにファイルを分けるような記載に変更．