StemClass.f90 Source File


Contents

Source Code


Source Code

module StemClass
    use, intrinsic :: iso_fortran_env
    use KinematicClass
    use FEMDomainClass
    implicit none
    
    type :: Stem_
        type(FEMDomain_)    ::  FEMDomain
        real(real64)             ::  Thickness,length,width
        real(real64)             ::  MaxThickness,Maxlength,Maxwidth
        real(real64)             ::  center_bottom(3),center_top(3)
        real(real64)             ::  radius_bottom(3),radius_top(3)
        real(real64)             ::  outer_normal_bottom(3),outer_normal_top(3)
        real(real64)             ::  rot_x = 0.0d0
        real(real64)             ::  rot_y = 0.0d0
        real(real64)             ::  rot_z = 0.0d0
        real(real64)             ::  disp_x = 0.0d0
        real(real64)             ::  disp_y = 0.0d0
        real(real64)             ::  disp_z = 0.0d0
        integer(int32)           ::  EdgeNodeID(4)
        integer(int32)           ::  EdgeElemID(4)
        real(real64)             ::  maxdiameter,mindiameter,minlength
        integer(int32),allocatable  :: I_planeNodeID(:)
        integer(int32),allocatable  :: I_planeElementID(:)
        integer(int32),allocatable  :: II_planeNodeID(:)
        integer(int32),allocatable  :: II_planeElementID(:)
        integer(int32)  :: A_PointNodeID
        integer(int32)  :: B_PointNodeID
        integer(int32)  :: A_PointElementID
        integer(int32)  :: B_PointElementID
        integer(int32)  :: xnum = 10
        integer(int32)  :: ynum = 10
        integer(int32)  :: znum = 10

        ! physical parameter
        real(real64),allocatable :: DryDensity(:)
        real(real64),allocatable :: WaterContent(:)

        integer(int32)             ::  Division
        type(Stem_),pointer ::  pStem
    contains
        procedure, public :: Init => initStem
        procedure, public :: rotate => rotateStem
        procedure, public :: grow => growStem
        procedure, public :: resize => resizeStem
        procedure, public :: move => moveStem
        procedure, public :: connect => connectStem
        procedure, public :: getCoordinate => getCoordinateStem
        procedure, public :: getLength => getLengthStem
        procedure, public :: gmsh => gmshStem
        procedure, public :: msh => mshStem
        procedure, public :: vtk => vtkStem
        procedure, public :: stl => stlStem
        procedure, public :: export => exportStem
        procedure, public :: getVolume => getVolumeStem
        procedure, public :: getBiomass => getBiomassStem
    end type
contains



! ########################################
subroutine initStem(obj,config,regacy,Thickness,length,width,MaxThickness,Maxlength,Maxwidth,rotx,roty,rotz,location)
    class(Stem_),intent(inout) :: obj
    real(real64),optional,intent(in)::  Thickness,length,width
    real(real64),optional,intent(in)::  MaxThickness,Maxlength,MaxWidth
    real(real64),optional,intent(in)::  rotx,roty,rotz,location(3)
    logical, optional,intent(in) :: regacy
    character(*),optional,intent(in) :: config
    type(IO_) :: stemconf,f
    character(200) :: fn,conf,line
    integer(int32),allocatable :: buf(:)
    integer(int32) :: id,rmc,n,node_id,node_id2,elemid,blcount,i,j
    real(real64) :: loc(3)
    logical :: debug=.false.

    ! 節を生成するためのスクリプトを開く
    if(.not.present(config)  .or. index(config,"json")==0 )then
        ! デフォルトの設定を生成
        if(debug) print *, "New stem-configuration >> stemconfig.json"
        call stemconf%open("stemconfig.json")
        write(stemconf%fh,*) '{'
        write(stemconf%fh,*) '   "type": "stem",'
        write(stemconf%fh,*) '   "minlength": 0.001,'
        write(stemconf%fh,*) '   "mindiameter": 0.001,'
        write(stemconf%fh,*) '   "maxlength": 0.07,'
        write(stemconf%fh,*) '   "maxdiameter": 0.01,'
        write(stemconf%fh,*) '   "drydensity": 0.0,'
        write(stemconf%fh,*) '   "watercontent": 0.0,'
        write(stemconf%fh,*) '   "xnum": 10,'
        write(stemconf%fh,*) '   "ynum": 10,'
        write(stemconf%fh,*) '   "znum": 10'
        write(stemconf%fh,*) '}'
        conf="stemconfig.json"
        call stemconf%close()
    else
        conf = trim(config)
    endif

    
    call stemconf%open(trim(conf),"r")
    blcount=0
    do
        read(stemconf%fh,'(a)') line
        if(debug) print *, trim(line)
        if( adjustl(trim(line))=="{" )then
            blcount=1
            cycle
        endif
        if( adjustl(trim(line))=="}" )then
            exit
        endif
        
        if(blcount==1)then
            
            if(index(line,"type")/=0 .and. index(line,"stem")==0 )then
                print *, "ERROR: This config-file is not for stem"
                return
            endif

            if(index(line,"maxlength")/=0 )then
                ! 生育ステージ
                rmc=index(line,",")
                ! カンマがあれば除く
                if(rmc /= 0)then
                    line(rmc:rmc)=" "
                endif
                id = index(line,":")
                read(line(id+1:),*) obj%maxlength
            endif


            if(index(line,"maxdiameter")/=0 )then
                ! 種子の長さ
                rmc=index(line,",")
                ! カンマがあれば除く
                if(rmc /= 0)then
                    line(rmc:rmc)=" "
                endif
                id = index(line,":")
                read(line(id+1:),*) obj%maxdiameter
            endif

            if(index(line,"minlength")/=0 )then
                ! 生育ステージ
                rmc=index(line,",")
                ! カンマがあれば除く
                if(rmc /= 0)then
                    line(rmc:rmc)=" "
                endif
                id = index(line,":")
                read(line(id+1:),*) obj%minlength
            endif


            if(index(line,"mindiameter")/=0 )then
                ! 種子の長さ
                rmc=index(line,",")
                ! カンマがあれば除く
                if(rmc /= 0)then
                    line(rmc:rmc)=" "
                endif
                id = index(line,":")
                read(line(id+1:),*) obj%mindiameter
            endif



            if(index(line,"xnum")/=0 )then
                ! 種子の長さ
                rmc=index(line,",")
                ! カンマがあれば除く
                if(rmc /= 0)then
                    line(rmc:rmc)=" "
                endif
                id = index(line,":")
                read(line(id+1:),*) obj%xnum
            endif



            if(index(line,"ynum")/=0 )then
                ! 種子の長さ
                rmc=index(line,",")
                ! カンマがあれば除く
                if(rmc /= 0)then
                    line(rmc:rmc)=" "
                endif
                id = index(line,":")
                read(line(id+1:),*) obj%ynum
            endif



            if(index(line,"znum")/=0 )then
                ! 種子の長さ
                rmc=index(line,",")
                ! カンマがあれば除く
                if(rmc /= 0)then
                    line(rmc:rmc)=" "
                endif
                id = index(line,":")
                read(line(id+1:),*) obj%znum
            endif
            cycle

        endif

    enddo
    call stemconf%close()
    
    
    ! グラフ構造とメッシュ構造を生成する。



    !
    !                                    <II>        
    !                     # # #           
    !                 #   # B       # #           
    !               #     #       #   #           
    !              # #     #           
    !              #      #    #      #           
    !              #      #    #      #           
    !              #      #    #      #           
    !              #      #    #      #           
    !              #      #           #           
    !              #      # D #           
    !              #    #            #             
    !              #  C     A  #   #   <I>           
    !              # #         # #                 
    !              # # # # # # #                        
    !                                             

    ! メッシュを生成
    
    call obj%FEMdomain%create(meshtype="Cube",x_num=obj%xnum,y_num=obj%ynum,z_num=obj%znum,&
    x_len=obj%mindiameter/2.0d0,y_len=obj%mindiameter/2.0d0,z_len=obj%minlength )

    ! initialize physical parameters
    obj%DryDensity = zeros( obj%FEMDomain%ne() )
    obj%watercontent = zeros( obj%FEMDomain%ne() )
    
    obj%DryDensity(:) = freal(stemconf%parse(conf,key1="drydensity"))
    obj%watercontent(:) = freal(stemconf%parse(conf,key1="watercontent"))

    ! <I>面に属する要素番号、節点番号、要素座標、節点座標のリストを生成
    obj%I_planeNodeID = obj%FEMdomain%mesh%getNodeList(zmax=0.0d0)
    obj%I_planeElementID = obj%FEMdomain%mesh%getElementList(zmax=0.0d0)
    
    ! <I>面に属する要素番号、節点番号、要素座標、節点座標のリストを生成
    obj%II_planeNodeID = obj%FEMdomain%mesh%getNodeList(zmin=obj%minlength)
    obj%II_planeElementID = obj%FEMdomain%mesh%getElementList(zmin=obj%minlength)
    
    buf   = obj%FEMDomain%mesh%getNodeList(&
        xmin=obj%mindiameter/2.0d0 - obj%mindiameter/dble(obj%xnum)/2.0d0 ,&
        xmax=obj%mindiameter/2.0d0 + obj%mindiameter/dble(obj%xnum)/2.0d0 ,&
        ymin=obj%mindiameter/2.0d0 - obj%mindiameter/dble(obj%ynum)/2.0d0 ,&
        ymax=obj%mindiameter/2.0d0 + obj%mindiameter/dble(obj%ynum)/2.0d0 ,&
        zmax=0.0d0)
    obj%A_PointNodeID = buf(1)

    buf   = obj%FEMDomain%mesh%getNodeList(&
        xmin=obj%mindiameter/2.0d0 - obj%mindiameter/dble(obj%xnum)/2.0d0 ,&
        xmax=obj%mindiameter/2.0d0 + obj%mindiameter/dble(obj%xnum)/2.0d0 ,&
        ymin=obj%mindiameter/2.0d0 - obj%mindiameter/dble(obj%ynum)/2.0d0 ,&
        ymax=obj%mindiameter/2.0d0 + obj%mindiameter/dble(obj%ynum)/2.0d0 ,&
        zmin=obj%minlength)
    obj%B_PointNodeID = buf(1)
    buf    = obj%FEMDomain%mesh%getElementList(&
        xmin=obj%mindiameter/2.0d0 - obj%mindiameter/dble(obj%xnum)/2.0d0 ,&
        xmax=obj%mindiameter/2.0d0 + obj%mindiameter/dble(obj%xnum)/2.0d0 ,&
        ymin=obj%mindiameter/2.0d0 - obj%mindiameter/dble(obj%ynum)/2.0d0 ,&
        ymax=obj%mindiameter/2.0d0 + obj%mindiameter/dble(obj%ynum)/2.0d0 ,&
        zmax=0.0d0)
    obj%A_PointElementID = buf(1)

    buf    = obj%FEMDomain%mesh%getElementList(&
        xmin=obj%mindiameter/2.0d0 - obj%mindiameter/dble(obj%xnum)/2.0d0 ,&
        xmax=obj%mindiameter/2.0d0 + obj%mindiameter/dble(obj%xnum)/2.0d0 ,&
        ymin=obj%mindiameter/2.0d0 - obj%mindiameter/dble(obj%ynum)/2.0d0 ,&
        ymax=obj%mindiameter/2.0d0 + obj%mindiameter/dble(obj%ynum)/2.0d0 ,&
        zmin=obj%minlength)
    obj%B_PointElementID = buf(1)

    if(debug) print *, obj%A_PointNodeID
    if(debug) print *, obj%B_PointNodeID
    if(debug) print *, obj%A_PointElementID
    if(debug) print *, obj%B_PointElementID
! デバッグ用
!    call f%open("I_phaseNodeID.txt")
!    do i=1,size(obj%I_planeNodeID)
!        write(f%fh,*) obj%femdomain%mesh%NodCoord( obj%I_planeNodeID(i) ,:)
!    enddo
!    call f%close()
!
!    call f%open("II_phaseNodeID.txt")
!    do i=1,size(obj%II_planeNodeID)
!        write(f%fh,*) obj%femdomain%mesh%NodCoord( obj%II_planeNodeID(i) ,:)
!    enddo
!    call f%close()
!
!    call f%open("I_phaseElementID.txt")
!    do i=1,size(obj%I_planeElementID)
!        do j=1,size(obj%femdomain%mesh%elemnod,2)
!            write(f%fh,*) obj%femdomain%mesh%NodCoord( &
!            obj%femdomain%mesh%elemnod(obj%I_planeElementID(i),j),:)
!        enddo
!    enddo
!    call f%close()
!
!    call f%open("II_phaseElementID.txt")
!    do i=1,size(obj%II_planeElementID)
!        do j=1,size(obj%femdomain%mesh%elemnod,2)
!            write(f%fh,*) obj%femdomain%mesh%NodCoord( &
!            obj%femdomain%mesh%elemnod(obj%II_planeElementID(i),j),:)
!        enddo
!    enddo
!    call f%close()
!    return

    ! Aについて、要素番号、節点番号、要素座標、節点座標のリストを生成




    !
    !           %%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
    !         %%                        %   %
    !        %%                    %      %%  
    !      %%                 %          %%    
    !     %%            %              %%      
    !     %%      %                  %%        
    !     %%                       %%          
    !       %%                  %%            
    !         %%%%%%%%%%%%%%%%%                               














    if( present(regacy))then
        if(regacy .eqv. .true.)then
            loc(:)=0.0d0
            if(present(location) )then
                loc(:)=location(:)
            endif
        
            obj%Thickness   = input(default=0.010d0,option= Thickness     )
            obj%length      = input(default=0.050d0,option= length      )
            obj%width       = input(default=0.010d0,option= width)
        
            obj%MaxThickness   = input(default=0.50d0  ,option=MaxThickness      )
            obj%Maxlength      = input(default=10.0d0  ,option=Maxlength       )
            obj%Maxwidth       = input(default=0.50d0  ,option=Maxwidth )
        
            obj%outer_normal_bottom(:)=0.0d0
            obj%outer_normal_bottom(1)=1.0d0
            obj%outer_normal_top(:)=0.0d0
            obj%outer_normal_top(1)=1.0d0
            ! rotate
            obj%outer_normal_Bottom(:) = Rotation3D(vector=obj%outer_normal_bottom,rotx=rotx,roty=roty,rotz=rotz)
            obj%outer_normal_top(:) = Rotation3D(vector=obj%outer_normal_top,rotx=rotx,roty=roty,rotz=rotz)
        
            obj%center_bottom(:)=loc(:)
            obj%center_top(:) = obj%center_bottom(:) + obj%length*obj%outer_normal_bottom(:)
        endif
    endif


end subroutine 
! ########################################

subroutine resize(obj,x,y,z)
    class(Stem_),intent(inout) :: obj
    real(real64),optional,intent(in) :: x,y,z

    call obj%femdomain%resize(x,y,z)

end subroutine

! ########################################
subroutine exportStem(obj,FileName,StemID)
    class(Stem_),intent(in)::obj
    character(*),intent(in) :: FileName
    integer(int32),optional,intent(inout) :: StemID
    real(real64) :: radius

    
    radius=0.50d0*obj%width+0.50d0*obj%Thickness
    
    open(13,file=FileName)
    write(13,'(A)') "//+"
    write(13,'(A)') 'SetFactory("OpenCASCADE");'
    write(13,*) "Cylinder(",input(default=1,option=StemID),") = {",&
    obj%center_bottom(1),",", obj%center_bottom(2),",", obj%center_bottom(3),",",&
    obj%center_top(1)-obj%center_bottom(1),",", obj%center_top(2)-obj%center_bottom(2),",",&
     obj%center_top(3)-obj%center_bottom(3),",",&
    radius,", 2*Pi};"
    close(13)
    StemID=StemID+1

end subroutine
! ########################################


! ########################################
recursive subroutine rotateStem(obj,x,y,z,reset)
    class(Stem_),intent(inout) :: obj
    real(real64),optional,intent(in) :: x,y,z
    logical,optional,intent(in) :: reset
    real(real64),allocatable :: origin1(:),origin2(:),disp(:)

    if(present(reset) )then
        if(reset .eqv. .true.)then
            call obj%femdomain%rotate(-obj%rot_x,-obj%rot_y,-obj%rot_z)
            obj%rot_x = 0.0d0
            obj%rot_y = 0.0d0
            obj%rot_z = 0.0d0
        endif
    endif

    origin1 = obj%getCoordinate("A")
    call obj%femdomain%rotate(x,y,z)
    obj%rot_x = obj%rot_x + input(default=0.0d0, option=x)
    obj%rot_y = obj%rot_y + input(default=0.0d0, option=y)
    obj%rot_z = obj%rot_z + input(default=0.0d0, option=z)
    origin2 = obj%getCoordinate("A")
    disp = origin1
    disp(:) = origin1(:) - origin2(:)
    call obj%femdomain%move(x=disp(1),y=disp(2),z=disp(3) )

    
end subroutine
! ########################################


! ########################################
recursive subroutine moveStem(obj,x,y,z,reset)
    class(Stem_),intent(inout) :: obj
    real(real64),optional,intent(in) :: x,y,z
    logical,optional,intent(in) :: reset
    real(real64),allocatable :: origin1(:),origin2(:),disp(:)

    if(present(reset) )then
        if(reset .eqv. .true.)then
            call obj%femdomain%move(-obj%disp_x,-obj%disp_y,-obj%disp_z)
            obj%disp_x = 0.0d0
            obj%disp_y = 0.0d0
            obj%disp_z = 0.0d0
        endif
    endif

    call obj%femdomain%move(x,y,z)
    obj%disp_x = obj%disp_x + input(default=0.0d0, option=x)
    obj%disp_y = obj%disp_y + input(default=0.0d0, option=y)
    obj%disp_z = obj%disp_z + input(default=0.0d0, option=z)
    
end subroutine
! ########################################

subroutine connectStem(obj,direct,stem)
    class(Stem_),intent(inout) :: obj,stem
    character(2),intent(in) :: direct
    real(real64),allocatable :: x1(:),x2(:),disp(:)

    if(direct=="->" .or. direct=="=>")then
        ! move obj to connect stem (stem is not moved.)
        x1 =  obj%getCoordinate("A")
        x2 = stem%getCoordinate("B")
        disp = x2 - x1
        call obj%move(x=disp(1),y=disp(2),z=disp(3) )
    endif


    if(direct=="<-" .or. direct=="<=")then
        ! move obj to connect stem (stem is not moved.)
        x1 = stem%getCoordinate("A")
        x2 =  obj%getCoordinate("B")
        disp = x2 - x1
        call stem%move(x=disp(1),y=disp(2),z=disp(3) )
    endif
end subroutine


! ########################################
function getCoordinateStem(obj,nodetype) result(ret)
    class(Stem_),intent(inout) :: obj
    character(*),intent(in) :: nodetype
    real(real64),allocatable :: ret(:)
    integer(int32) :: dimnum

    dimnum = size(obj%femdomain%mesh%nodcoord,2)
    allocate(ret(dimnum) )
    if( trim(nodetype)=="A" .or. trim(nodetype)=="a")then
        ret = obj%femdomain%mesh%nodcoord(obj%A_PointNodeID,:)
    endif
    if( trim(nodetype)=="B" .or. trim(nodetype)=="B")then
        ret = obj%femdomain%mesh%nodcoord(obj%B_PointNodeID,:)
    endif

end function
! ########################################

subroutine gmshStem(obj,name)
    class(Stem_),intent(inout) :: obj
    character(*),intent(in) ::name
    if(obj%femdomain%mesh%empty() )then
        return
    endif
    
    call obj%femdomain%gmsh(Name=name)
end subroutine

subroutine mshStem(obj,name)
    class(Stem_),intent(inout) :: obj
    character(*),intent(in) ::name
    if(obj%femdomain%mesh%empty() )then
        return
    endif
    
    call obj%femdomain%msh(Name=name)
end subroutine

subroutine vtkStem(obj,name)
    class(Stem_),intent(inout) :: obj
    character(*),intent(in) ::name
    
    if(obj%femdomain%mesh%empty() )then
        return
    endif
    
    call obj%femdomain%vtk(Name=name)
end subroutine

subroutine stlStem(obj,name)
    class(Stem_),intent(inout) :: obj
    character(*),intent(in) ::name
    if(obj%femdomain%mesh%empty() )then
        return
    endif
    
    call obj%femdomain%stl(Name=name)
end subroutine
! ########################################
subroutine resizeStem(obj,x,y,z)
    class(Stem_),optional,intent(inout) :: obj
    real(real64),optional,intent(in) :: x,y,z
    real(real64),allocatable :: origin1(:), origin2(:),disp(:)

    origin1 = obj%getCoordinate("A")
    call obj%femdomain%resize(x_len=x,y_len=y,z_len=z)
    origin2 = obj%getCoordinate("A")
    disp = origin1 - origin2
    call obj%move(x=disp(1),y=disp(2),z=disp(3) )
end subroutine
! ########################################

! ########################################
subroutine growStem(obj,length,length_rate,width_rate)
    class(Stem_),intent(inout) :: obj 
    real(real64),optional,intent(in) :: length,length_rate,width_rate
    real(real64) :: length_r,width_r,l_0,w_0,clength
    real(real64),allocatable :: origin(:),top(:),n1(:),coord(:),center(:),vert(:)
    integer(int32) :: i

    origin = obj%getCoordinate("A")
    top    = obj%getCoordinate("B")
    l_0    = sqrt(dot_product(top-origin, top-origin) ) 
    n1     = origin
    n1     = top - origin
    n1     = 1.0d0/norm(n1)*n1
    coord  = origin
    
    ! length-ratio = new length / old length
    if(present(length) )then
        length_r = length/l_0
    elseif(present(length_rate) )then
        length_r = length_rate
    else
        length_r = 1.0d0
    endif

    width_r = input(default=1.0d0, option=width_rate)

    ! enlong & fatten
    do i=1,obj%femdomain%nn()
        coord(:) = obj%femdomain%mesh%nodcoord(i,:) - origin(:)
        center   = coord
        clength   = dot_product(coord, n1)
        center(:)  = clength*n1(:)
        vert = coord - center
        ! origin -> center -> current coordinate
        coord(:) = length_r*center(:) + width_r*vert(:)
        obj%femdomain%mesh%nodcoord(i,:) = origin(:) + coord(:)
    enddo

    

end subroutine
! ########################################

! ########################################
subroutine rescaleStem(obj,x,y,z)
    class(Stem_),optional,intent(inout) :: obj
    real(real64),optional,intent(in) :: x,y,z
    real(real64),allocatable :: origin1(:), origin2(:),disp(:)

    origin1 = obj%getCoordinate("A")
    call obj%femdomain%resize(x_rate=x,y_rate=y,z_rate=z)
    origin2 = obj%getCoordinate("A")
    disp = origin1 - origin2
    call obj%move(x=disp(1),y=disp(2),z=disp(3) )
end subroutine
! ########################################


function getLengthStem(obj) result(ret)
        class(Stem_),intent(in) :: obj
    real(real64) :: ret

    if(obj%femdomain%mesh%empty() )then
        ret = 0.0d0
    else
        ret = norm(&
            obj%femdomain%mesh%nodcoord(obj%A_PointNodeID,: ) &
            - obj%femdomain%mesh%nodcoord(obj%B_PointNodeID,:) )
    endif

end function


function getVolumeStem(obj) result(ret)
    class(Stem_),intent(in) :: obj
    real(real64) :: ret
    integer(int32) :: i,j
    
    ret =0.0d0
    if(obj%femdomain%mesh%empty() ) then
        return
    endif

    do i=1,obj%femdomain%ne()
        ret = ret + obj%femdomain%getVolume(elem=i)
    enddo

end function


function getBiomassStem(obj) result(ret)
    class(Stem_),intent(in) :: obj
    real(real64) :: ret
    integer(int32) :: i,j

    ret =0.0d0
    if(obj%femdomain%mesh%empty() ) then
        return
    endif
    
    do i=1,obj%femdomain%ne()
        ret = ret + obj%femdomain%getVolume(elem=i)*obj%drydensity(i)
    enddo

end function


end module