module GrapeClass use LeafClass use StemClass use RootClass implicit none type :: Grape_ ! 節-節点データ構造 type(Mesh_) :: struct integer(int32), allocatable :: leaf2stem(:, :) integer(int32), allocatable :: stem2stem(:, :) integer(int32), allocatable :: root2stem(:, :) integer(int32), allocatable :: root2root(:, :) real(real64) :: mainstem_length real(real64) :: mainstem_width integer(int32) :: mainstem_node real(real64) :: mainroot_length real(real64) :: mainroot_width integer(int32) :: mainroot_node integer(int32) :: num_branch integer(int32) :: num_branch_node integer(int32) :: num_branch_root integer(int32) :: num_branch_root_node real(real64) :: ms_angle_ave = 0.0d0 real(real64) :: ms_angle_sig = 0.0d0 real(real64), allocatable :: br_angle_ave_x(:) real(real64), allocatable :: br_angle_sig_x(:) real(real64), allocatable :: br_angle_ave_z(:) real(real64), allocatable :: br_angle_sig_z(:) real(real64), allocatable :: peti_size_ave(:) real(real64), allocatable :: peti_size_sig(:) real(real64), allocatable :: peti_width_ave(:) real(real64), allocatable :: peti_width_sig(:) real(real64), allocatable :: peti_angle_ave(:) real(real64), allocatable :: peti_angle_sig(:) real(real64), allocatable :: leaf_thickness_ave(:) real(real64), allocatable :: leaf_thickness_sig(:) real(real64), allocatable :: leaf_angle_ave(:) real(real64), allocatable :: leaf_angle_sig(:) real(real64), allocatable :: leaf_length_ave(:) real(real64), allocatable :: leaf_length_sig(:) real(real64), allocatable :: leaf_width_ave(:) real(real64), allocatable :: leaf_width_sig(:) integer(int32), allocatable :: br_node(:) integer(int32), allocatable :: br_from(:) real(real64), allocatable :: br_length(:) real(real64), allocatable :: br_width(:) integer(int32) :: num_leaf integer(int32) :: num_stem integer(int32) :: num_root ! 器官オブジェクト配列 type(FEMDomain_), allocatable :: leaf_list(:) type(FEMDomain_), allocatable :: stem_list(:) type(FEMDomain_), allocatable :: root_list(:) character(:), allocatable :: LeafSurfaceData type(Leaf_), allocatable :: Leaf(:) type(Stem_), allocatable :: Stem(:) type(Root_), allocatable :: Root(:) contains procedure :: create => createGrape procedure, public :: msh => mshGrape procedure, public :: vtk => vtkGrape procedure, public :: stl => stlGrape procedure, public :: json => jsonGrape procedure, public :: move => moveGrape end type contains ! ############################################################# subroutine createGrape(obj, config, LeafSurfaceData) class(Grape_), intent(inout) :: obj character(*), intent(in) :: config character(*), optional, intent(in) ::LeafSurfaceData character(:), allocatable :: line type(IO_) :: grapeconfig type(Random_) :: random integer(int32)::i, n, j, k, num_leaf, num_stem_node, num_branch_branch if (present(LeafSurfaceData)) then obj%LeafSurfaceData = LeafSurfaceData else obj%LeafSurfaceData = grapeconfig%parse(config, key1="LeafSurfaceData") end if obj%mainstem_length = freal(grapeconfig%parse(config, key1="Mainstem", key2="Length")) obj%mainstem_width = freal(grapeconfig%parse(config, key1="Mainstem", key2="Width")) obj%mainstem_node = fint(grapeconfig%parse(config, key1="Mainstem", key2="Node")) obj%ms_angle_ave = freal(grapeconfig%parse(config, key1="Mainstem", key2="ms_angle_ave")) obj%ms_angle_sig = freal(grapeconfig%parse(config, key1="Mainstem", key2="ms_angle_sig")) ! get number of branch && number of node obj%num_branch = 1 obj%num_branch_node = 0 do line = grapeconfig%parse(config, key1="Branch#"//str(obj%num_branch), key2="Node") if (len(trim(line)) == 0) then obj%num_branch = obj%num_branch - 1 exit else obj%num_branch = obj%num_branch + 1 obj%num_branch_node = obj%num_branch_node + fint(line) ! Further branch cycle end if end do allocate (obj%br_node(obj%num_branch)) allocate (obj%br_from(obj%num_branch)) allocate (obj%br_length(obj%num_branch)) allocate (obj%br_width(obj%num_branch)) allocate (obj%br_angle_ave_x(obj%num_branch)) allocate (obj%br_angle_sig_x(obj%num_branch)) allocate (obj%br_angle_ave_z(obj%num_branch)) allocate (obj%br_angle_sig_z(obj%num_branch)) allocate (obj%peti_size_ave(obj%num_branch)) allocate (obj%peti_size_sig(obj%num_branch)) allocate (obj%peti_width_ave(obj%num_branch)) allocate (obj%peti_width_sig(obj%num_branch)) allocate (obj%peti_angle_ave(obj%num_branch)) allocate (obj%peti_angle_sig(obj%num_branch)) allocate (obj%leaf_thickness_ave(obj%num_branch)) allocate (obj%leaf_thickness_sig(obj%num_branch)) allocate (obj%leaf_angle_ave(obj%num_branch)) allocate (obj%leaf_angle_sig(obj%num_branch)) allocate (obj%leaf_length_ave(obj%num_branch)) allocate (obj%leaf_length_sig(obj%num_branch)) allocate (obj%leaf_width_ave(obj%num_branch)) allocate (obj%leaf_width_sig(obj%num_branch)) do i = 1, obj%num_branch line = grapeconfig%parse(config, key1="Branch#"//str(i), key2="Node") obj%br_node(i) = fint(line) line = grapeconfig%parse(config, key1="Branch#"//str(i), key2="From") obj%br_from(i) = fint(line) line = grapeconfig%parse(config, key1="Branch#"//str(i), key2="Length") obj%br_length(i) = freal(line) line = grapeconfig%parse(config, key1="Branch#"//str(i), key2="Width") obj%br_width(i) = freal(line) obj%br_angle_ave_x(i) = freal(grapeconfig%parse( & config, key1="Branch#"//str(i), key2="br_angle_ave_x")) obj%br_angle_ave_z(i) = freal(grapeconfig%parse( & config, key1="Branch#"//str(i), key2="br_angle_ave_z")) obj%br_angle_sig_x(i) = freal(grapeconfig%parse( & config, key1="Branch#"//str(i), key2="br_angle_sig_x")) obj%br_angle_sig_z(i) = freal(grapeconfig%parse( & config, key1="Branch#"//str(i), key2="br_angle_sig_z")) obj%peti_size_ave(i) = freal(grapeconfig%parse( & config, key1="Branch#"//str(i), key2="peti_size_ave")) obj%peti_size_sig(i) = freal(grapeconfig%parse( & config, key1="Branch#"//str(i), key2="peti_size_sig")) obj%peti_width_ave(i) = freal(grapeconfig%parse( & config, key1="Branch#"//str(i), key2="peti_width_ave")) obj%peti_width_sig(i) = freal(grapeconfig%parse( & config, key1="Branch#"//str(i), key2="peti_width_sig")) obj%peti_angle_ave(i) = freal(grapeconfig%parse( & config, key1="Branch#"//str(i), key2="peti_angle_ave")) obj%peti_angle_sig(i) = freal(grapeconfig%parse( & config, key1="Branch#"//str(i), key2="peti_angle_sig")) obj%leaf_thickness_ave(i) = freal(grapeconfig%parse( & config, key1="Branch#"//str(i), key2="leaf_thickness_ave")) obj%leaf_thickness_sig(i) = freal(grapeconfig%parse( & config, key1="Branch#"//str(i), key2="leaf_thickness_sig")) obj%leaf_angle_ave(i) = freal(grapeconfig%parse( & config, key1="Branch#"//str(i), key2="leaf_angle_ave")) obj%leaf_angle_sig(i) = freal(grapeconfig%parse( & config, key1="Branch#"//str(i), key2="leaf_angle_sig")) obj%leaf_length_ave(i) = freal(grapeconfig%parse( & config, key1="Branch#"//str(i), key2="leaf_length_ave")) obj%leaf_length_sig(i) = freal(grapeconfig%parse( & config, key1="Branch#"//str(i), key2="leaf_length_sig")) obj%leaf_width_ave(i) = freal(grapeconfig%parse( & config, key1="Branch#"//str(i), key2="leaf_width_ave")) obj%leaf_width_sig(i) = freal(grapeconfig%parse( & config, key1="Branch#"//str(i), key2="leaf_width_sig")) end do obj%mainroot_length = freal(grapeconfig%parse(config, key1="Mainroot", key2="Length")) obj%mainroot_width = freal(grapeconfig%parse(config, key1="Mainroot", key2="Width")) obj%mainroot_node = fint(grapeconfig%parse(config, key1="Mainroot", key2="Node")) ! get number of branch && number of node obj%num_branch_root = 1 obj%num_branch_root_node = 0 do line = grapeconfig%parse(config, key1="Branchroot#"//str(obj%num_branch_root), key2="Node") if (len(trim(line)) == 0) then obj%num_branch_root = obj%num_branch_root - 1 exit else obj%num_branch_root = obj%num_branch_root + 1 obj%num_branch_root_node = obj%num_branch_root_node + fint(line) cycle end if end do obj%num_leaf = obj%num_branch_node + obj%mainstem_node obj%num_stem = obj%num_branch_node + obj%mainstem_node obj%num_root = obj%num_branch_root_node + obj%mainroot_node allocate (obj%leaf_list(obj%num_leaf)) allocate (obj%stem_list(obj%num_stem)) allocate (obj%root_list(obj%num_root)) allocate (obj%leaf(obj%num_leaf*2)) allocate (obj%stem(obj%num_stem*2)) allocate (obj%root(obj%num_root*2)) obj%leaf2stem = zeros(obj%num_leaf, obj%num_stem*2) obj%stem2stem = zeros(obj%num_stem*2, obj%num_stem*2) obj%root2stem = zeros(obj%num_root, obj%num_stem*2) obj%root2root = zeros(obj%num_root, obj%num_root) ! set mainstem do i = 1, obj%mainstem_node call obj%stem(i)%init() call obj%stem(i)%resize( & x=obj%mainstem_width, & y=obj%mainstem_width, & z=obj%mainstem_length/dble(obj%mainstem_node) & ) call obj%stem(i)%rotate( & x=radian(random%gauss(mu=obj%ms_angle_ave, sigma=obj%ms_angle_sig)), & y=radian(random%gauss(mu=obj%ms_angle_ave, sigma=obj%ms_angle_sig)), & z=radian(random%gauss(mu=obj%ms_angle_ave, sigma=obj%ms_angle_sig)) & ) end do do i = 1, obj%mainstem_node - 1 call obj%stem(i + 1)%connect("=>", obj%stem(i)) obj%stem2stem(i + 1, i) = 1 end do ! set branches k = obj%mainstem_node do i = 1, size(obj%br_node) do j = 1, obj%br_node(i) k = k + 1 call obj%stem(k)%init() call obj%stem(k)%resize( & x=obj%br_width(i), & y=obj%br_width(i), & z=obj%br_length(i)/dble(obj%br_node(i)) & ) call obj%stem(k)%rotate( & x=radian(random%gauss(mu=obj%br_angle_ave_x(i), sigma=obj%br_angle_sig_x(i))), & y=0.0d0, & z=radian(random%gauss(mu=obj%br_angle_ave_z(i), sigma=obj%br_angle_sig_z(i))) & ) if (j == 1) then call obj%stem(k)%connect("=>", obj%stem(obj%br_from(i))) obj%stem2stem(k, obj%br_from(i)) = 1 else call obj%stem(k)%connect("=>", obj%stem(k - 1)) obj%stem2stem(k, k - 1) = 1 end if end do end do ! ! peti and leaf num_stem_node = k num_leaf = 0 do i = 1, size(obj%br_node) do j = 1, obj%br_node(i) ! ３複葉 ! add peti num_stem_node = num_stem_node + 1 call obj%stem(num_stem_node)%init() call obj%stem(num_stem_node)%resize( & x=random%gauss(mu=obj%peti_width_ave(i), sigma=obj%peti_width_sig(i)), & y=random%gauss(mu=obj%peti_width_ave(i), sigma=obj%peti_width_sig(i)), & z=random%gauss(mu=obj%peti_size_ave(i), sigma=obj%peti_size_sig(i)) & ) call obj%stem(num_stem_node)%rotate( & x=radian(random%gauss(mu=obj%peti_angle_ave(i), sigma=obj%peti_angle_sig(i))), & y=0.0d0, & z=radian(360.0d0*random%random()) & ) !obj%leaf2stem(num_stem_node,i) = 1 if (i == 1) then call obj%stem(num_stem_node)%connect("=>", obj%stem(j + obj%mainroot_node)) obj%stem2stem(num_stem_node, j + obj%mainroot_node) = 1 else call obj%stem(num_stem_node)%connect("=>", obj%stem(j + sum(obj%br_node(1:i - 1)) + obj%mainroot_node)) obj%stem2stem(num_stem_node, j + sum(obj%br_node(1:i - 1)) + obj%mainroot_node) = 1 end if ! add leaves num_leaf = num_leaf + 1 !call obj%leaf(num_leaf)%create(filename=obj%LeafSurfaceData) call obj%leaf(num_leaf)%create(filename=obj%LeafSurfaceData) call obj%leaf(num_leaf)%femdomain%resize( & z=random%gauss(mu=obj%leaf_thickness_ave(i), sigma=obj%leaf_thickness_sig(i)), & x=random%gauss(mu=obj%leaf_length_ave(i), sigma=obj%leaf_length_sig(i)), & y=random%gauss(mu=obj%leaf_width_ave(i), sigma=obj%leaf_width_sig(i)) & ) call obj%leaf(num_leaf)%femdomain%rotate( & x=radian(random%gauss(mu=obj%leaf_angle_ave(i), sigma=obj%leaf_angle_sig(i))), & y=0.0d0, & z=radian(random%random()*360.0d0) & ) if (.not. allocated(obj%leaf(num_leaf)%i_planenodeid)) then obj%leaf(num_leaf)%i_planenodeid = [1] end if call obj%leaf(num_leaf)%connect("=>", obj%stem(num_stem_node)) obj%leaf2stem(num_leaf, num_stem_node) = 1 end do end do end subroutine ! ############################################################# ! ######################################## subroutine mshGrape(obj, name, num_threads) class(Grape_), intent(inout) :: obj character(*), intent(in) :: name integer(int32), optional, intent(in) :: num_threads integer(int32) :: i, n n = input(default=1, option=num_threads) !$OMP parallel num_threads(n) private(i) !$OMP do do i = 1, size(obj%stem) !if(obj%stem(i)%femdomain%mesh%empty() .eqv. .false. )then call obj%stem(i)%msh(name=name//"_stem"//str(i)) !endif end do !$OMP end do !$OMP end parallel !$OMP parallel num_threads(n) private(i) !$OMP do do i = 1, size(obj%root) !if(obj%root(i)%femdomain%mesh%empty() .eqv. .false. )then call obj%root(i)%msh(name=name//"_root"//str(i)) !endif end do !$OMP end do !$OMP end parallel !$OMP parallel num_threads(n) private(i) !$OMP do do i = 1, size(obj%leaf) !if(obj%leaf(i)%femdomain%mesh%empty() .eqv. .false. )then call obj%leaf(i)%msh(name=name//"_leaf"//str(i)) !endif end do !$OMP end do !$OMP end parallel end subroutine ! ######################################## ! ######################################## subroutine vtkGrape(obj, name, num_threads) class(Grape_), intent(inout) :: obj character(*), intent(in) :: name integer(int32), optional, intent(in) :: num_threads integer(int32) :: i, n n = input(default=1, option=num_threads) if (allocated(obj%stem)) then !$OMP parallel num_threads(n) private(i) !$OMP do do i = 1, size(obj%stem) !if(obj%stem(i)%femdomain%mesh%empty() .eqv. .false. )then call obj%stem(i)%vtk(name=name//"_stem"//str(i)) !endif end do !$OMP end do !$OMP end parallel end if if (allocated(obj%root)) then !$OMP parallel num_threads(n) private(i) !$OMP do do i = 1, size(obj%root) !if(obj%root(i)%femdomain%mesh%empty() .eqv. .false. )then call obj%root(i)%vtk(name=name//"_root"//str(i)) !endif end do !$OMP end do !$OMP end parallel end if if (allocated(obj%leaf)) then !$OMP parallel num_threads(n) private(i) !$OMP do do i = 1, size(obj%leaf) !if(obj%leaf(i)%femdomain%mesh%empty() .eqv. .false. )then call obj%leaf(i)%vtk(name=name//"_leaf"//str(i)) !endif end do !$OMP end do !$OMP end parallel end if end subroutine ! ######################################## ! ######################################## subroutine jsonGrape(obj, name) class(Grape_), intent(inout) :: obj character(*), intent(in) :: name integer(int32) :: i, countnum type(IO_) :: f call f%open(name//".json") call f%write("{") countnum = 0 do i = 1, size(obj%stem) if (obj%stem(i)%femdomain%mesh%empty() .eqv. .false.) then countnum = countnum + 1 call f%write('"'//"stem"//str(i)//'":') call obj%stem(i)%femdomain%json(name=name//"_stem"//str(i), fh=f%fh, endl=.false.) end if end do call f%write('"num_stem":'//str(countnum)//',') countnum = 0 do i = 1, size(obj%root) if (obj%root(i)%femdomain%mesh%empty() .eqv. .false.) then countnum = countnum + 1 call f%write('"'//"root"//str(i)//'":') call obj%root(i)%femdomain%json(name=name//"_root"//str(i), fh=f%fh, endl=.false.) end if end do call f%write('"num_root":'//str(countnum)//',') countnum = 0 do i = 1, size(obj%leaf) if (obj%leaf(i)%femdomain%mesh%empty() .eqv. .false.) then countnum = countnum + 1 call f%write('"'//"leaf"//str(i)//'":') call obj%leaf(i)%femdomain%json(name=name//"_leaf"//str(i), fh=f%fh, endl=.false.) end if end do call f%write('"num_leaf":'//str(countnum)//',') call f%write('"return_Grape":0') call f%write("}") call f%close() end subroutine ! ######################################## ! ######################################## subroutine stlGrape(obj, name, num_threads) class(Grape_), intent(inout) :: obj character(*), intent(in) :: name integer(int32), optional, intent(in) :: num_threads integer(int32) :: i, n n = input(default=1, option=num_threads) !call execute_command_line("echo ' ' > "//name//".stl") !$OMP parallel num_threads(n) private(i) !$OMP do do i = 1, size(obj%stem) if (obj%stem(i)%femdomain%mesh%empty() .eqv. .false.) then call obj%stem(i)%stl(name=name//"_stem"//str(i)) !call execute_command_line("cat "//name//"_stem"//str(i)//"_000001.stl >> "//name//".stl") end if end do !$OMP end do !$OMP end parallel !$OMP parallel num_threads(n) private(i) !$OMP do do i = 1, size(obj%root) if (obj%root(i)%femdomain%mesh%empty() .eqv. .false.) then call obj%root(i)%stl(name=name//"_root"//str(i)) !call execute_command_line("cat "//name//"_root"//str(i)//"_000001.stl >> "//name//".stl") end if end do !$OMP end do !$OMP end parallel !$OMP parallel num_threads(n) private(i) !$OMP do do i = 1, size(obj%leaf) if (obj%leaf(i)%femdomain%mesh%empty() .eqv. .false.) then call obj%leaf(i)%stl(name=name//"_leaf"//str(i)) !call execute_command_line("cat "//name//"_leaf"//str(i)//"_000001.stl >> "//name//".stl") end if end do !$OMP end do !$OMP end parallel end subroutine ! ######################################## ! ######################################## subroutine moveGrape(obj, x, y, z) class(Grape_), intent(inout) :: obj real(real64), optional, intent(in) :: x, y, z integer(int32) :: i do i = 1, size(obj%stem) if (obj%stem(i)%femdomain%mesh%empty() .eqv. .false.) then call obj%stem(i)%move(x=x, y=y, z=z) end if end do do i = 1, size(obj%root) if (obj%root(i)%femdomain%mesh%empty() .eqv. .false.) then call obj%root(i)%move(x=x, y=y, z=z) end if end do do i = 1, size(obj%leaf) if (obj%leaf(i)%femdomain%mesh%empty() .eqv. .false.) then call obj%leaf(i)%move(x=x, y=y, z=z) end if end do end subroutine ! ######################################## end module GrapeClass