module SeedClass
   use, intrinsic :: iso_fortran_env
   use MathClass
   use RandomClass
   use LsystemClass
   use FEMDomainClass
   !use PreprocessingClass
   use sim
   implicit none

   type :: Seed_
      type(FEMDomain_) :: FEMDomain ! regacy

      type(WaterAbsorption_) :: seedDomain
      type(FEMDomain_) :: water
      type(FEMDomain_) :: tissue

      type(MaterialProp_):: Permiability, YoungsModulus, PoissonRatio
      type(MaterialProp_):: density, Cohesion, phi, psi
      type(MaterialProp_):: a_Psi
      type(MaterialProp_):: a_P
      type(MaterialProp_):: theta_eq
      type(MaterialProp_):: Psi_eq
      type(MaterialProp_):: a_E
      type(MaterialProp_):: a_v
      type(MaterialProp_):: E_eq
      type(MaterialProp_):: v_eq

      type(Boundary_) :: disp_x
      type(Boundary_) :: disp_y
      type(Boundary_) :: disp_z
      type(Boundary_) :: traction_x
      type(Boundary_) :: traction_y
      type(Boundary_) :: traction_z
      type(Boundary_) :: flux
      type(Boundary_) :: const

      integer(int32) :: num_of_seed = 1 ! num of seed
      real(real64) :: mass = 0.0d0 ! seed mass g/cm^3
      real(real64) :: water_content = 0.0d0 ! seed water_content %
      real(real64) :: radius = 0.0d0 ! seed radius (cm)
      real(real64) :: width1 = 0.00d0
      real(real64) :: width2 = 0.00d0
      real(real64) :: width3 = 0.0d0
      real(real64) :: width1_origin = 0.00d0
      real(real64) :: width2_origin = 0.00d0
      real(real64) :: width3_origin = 0.0d0
      real(real64) :: location(3) = 0.0d0 ! seed location (x,y,z)

      character(200) :: path = " "
      character(200) :: name = " "

      ! not exported

      integer(int32) :: step = 0
   contains
      procedure :: create => createSeed
      !procedure :: create_seed_shape  => create_seed_shapeSeed
      procedure :: env => envSeed
      procedure :: material => materialSeed
      procedure :: grow => growSeed
      procedure :: move => moveSeed
      procedure :: rotate => rotateSeed
      procedure :: gmsh => gmshSeed

      procedure :: init => initSeed
      procedure :: import => importSeed
      procedure :: createMesh => createMeshSeed
      procedure :: export => exportSeed
      procedure :: save => saveSeed
      procedure :: open => openSeed
      procedure :: flush => flushSeed
      procedure :: show => showSeed
      procedure :: load => loadSeed
      procedure :: remove => removeSeed

      ! measure
      procedure :: length => lengthSeed

      procedure :: convertMeshType => convertMeshTypeSeed
      procedure :: result => resultSeed

      !MPI

   end type
contains

!#####################################
   function lengthSeed(obj) result(length)
      class(Seed_), intent(inout) :: obj
      real(real64)::length(3)

      length(:) = obj%tissue%length()

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

!#####################################
   subroutine resultSeed(obj, path, name, step)
      class(Seed_), intent(inout) :: obj
      character(*), intent(in) :: path, name
      integer(int32), intent(in) :: step

      call obj%seedDomain%result(path, name, step)
   end subroutine
!#####################################

!#####################################
   subroutine removeSeed(obj)
      class(Seed_), intent(inout)::obj

      call obj%FEMDomain%remove() ! regacy
      call obj%seedDomain%remove()
      call obj%water%remove()
      call obj%tissue%remove()
      call obj%Permiability%remove()
      call obj%YoungsModulus%remove()
      call obj%PoissonRatio%remove()
      call obj%density%remove()
      call obj%Cohesion%remove()
      call obj%phi%remove()
      call obj%psi%remove()
      call obj%a_Psi%remove()
      call obj%a_P%remove()
      call obj%theta_eq%remove()
      call obj%Psi_eq%remove()
      call obj%a_E%remove()
      call obj%a_v%remove()
      call obj%E_eq%remove()
      call obj%v_eq%remove()
      call obj%disp_x%remove()
      call obj%disp_y%remove()
      call obj%disp_z%remove()
      call obj%traction_x%remove()
      call obj%traction_y%remove()
      call obj%traction_z%remove()
      call obj%flux%remove()
      call obj%const%remove()

      obj%num_of_seed = 1 ! num of seed
      obj%mass = 0.0d0 ! seed mass g/cm^3
      obj%water_content = 0.0d0 ! seed water_content %
      obj%radius = 0.0d0 ! seed radius (cm)
      obj%width1 = 0.00d0
      obj%width2 = 0.00d0
      obj%width3 = 0.0d0
      obj%width1_origin = 0.00d0
      obj%width2_origin = 0.00d0
      obj%width3_origin = 0.0d0
      obj%location(3) = 0.0d0 ! seed location (x,y,z)

      obj%path = " "
      obj%name = " "

      ! not exported

      obj%step = 0

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

!########################################################
   subroutine loadSeed(obj, path, name)
      class(Seed_), intent(inout)::obj
      character(*), optional, intent(in) :: path
      character(*), optional, intent(in) :: name

      call obj%remove()
      call obj%open(path=path, name=name)
   end subroutine
!########################################################

!########################################################
   subroutine openSeed(obj, path, name)
      class(Seed_), intent(inout)::obj
      character(*), optional, intent(in) :: path
      character(*), optional, intent(in) :: name
      character(200) :: data_type
      type(IO_) :: f

      if (present(path)) then
         obj%path = path
      else
         obj%path = "home/"
      end if

      if (present(name)) then
         obj%name = name
      else
         obj%name = "untitled"
      end if

      ! check data-type
      call f%open(obj%path//obj%name//"/type")
      call f%read(data_type)
      call f%close()
      if (data_type /= "seed") then
         print *, "openSeed >> ERROR :: "//obj%path//obj%name//"/type /= seed"
         stop
      end if

      call f%open(obj%path//obj%name//"/input/seed/seedDomain/", "seed", ".prop")
      read (f%fh, *) obj%num_of_seed ! num of seed
      read (f%fh, *) obj%mass ! seed mass g/cm^3
      read (f%fh, *) obj%water_content ! seed water_content %
      read (f%fh, *) obj%radius ! seed radius (cm)
      read (f%fh, *) obj%width1
      read (f%fh, *) obj%width2
      read (f%fh, *) obj%width3
      read (f%fh, *) obj%width1_origin
      read (f%fh, *) obj%width2_origin
      read (f%fh, *) obj%width3_origin
      read (f%fh, *) obj%location(1) ! seed location (x,y,z)
      read (f%fh, *) obj%location(2) ! seed location (x,y,z)
      read (f%fh, *) obj%location(3) ! seed location (x,y,z)
      read (f%fh, '(A)') obj%path
      obj%path = obj%path
      read (f%fh, '(A)') obj%name
      obj%name = obj%name
      call f%close()

      ! only for SeedClass
      ! FEMDomain
      call obj%water%open(path=obj%path//obj%name//"/input/seed", name="water")

      call obj%tissue%open(path=obj%path//obj%name//"/input/seed", name="tissue")

      ! WaterAbsorption
      call obj%seedDomain%open(path=obj%path//obj%name//"/input/seed", name="seedDomain")

      ! Material
      call obj%Permiability%open(path=obj%path//obj%name//"/input/seed/seedDomain", name="Permiability")
      call obj%YoungsModulus%open(path=obj%path//obj%name//"/input/seed/seedDomain", name="YoungsModulus")
      call obj%PoissonRatio%open(path=obj%path//obj%name//"/input/seed/seedDomain", name="PoissonRatio")
      call obj%density%open(path=obj%path//obj%name//"/input/seed/seedDomain", name="density")
      call obj%Cohesion%open(path=obj%path//obj%name//"/input/seed/seedDomain", name="Cohesion")
      call obj%phi%open(path=obj%path//obj%name//"/input/seed/seedDomain", name="phi")
      call obj%psi%open(path=obj%path//obj%name//"/input/seed/seedDomain", name="phi")
      call obj%a_Psi%open(path=obj%path//obj%name//"/input/seed/seedDomain", name="a_Psi")
      call obj%a_P%open(path=obj%path//obj%name//"/input/seed/seedDomain", name="a_P")
      call obj%theta_eq%open(path=obj%path//obj%name//"/input/seed/seedDomain", name="theta_eq")
      call obj%Psi_eq%open(path=obj%path//obj%name//"/input/seed/seedDomain", name="Psi_eq")
      call obj%a_E%open(path=obj%path//obj%name//"/input/seed/seedDomain", name="a_E")
      call obj%a_v%open(path=obj%path//obj%name//"/input/seed/seedDomain", name="a_v")
      call obj%E_eq%open(path=obj%path//obj%name//"/input/seed/seedDomain", name="E_eq")
      call obj%v_eq%open(path=obj%path//obj%name//"/input/seed/seedDomain", name="v_eq")

      ! Boundary
      call obj%disp_x%open(path=obj%path//obj%name//"/input/seed/seedDomain", name="disp_x")
      call obj%disp_y%open(path=obj%path//obj%name//"/input/seed/seedDomain", name="disp_y")
      call obj%disp_z%open(path=obj%path//obj%name//"/input/seed/seedDomain", name="disp_z")
      call obj%const%open(path=obj%path//obj%name//"/input/seed/seedDomain", name="const")

      call obj%traction_x%open(path=obj%path//obj%name//"/input/seed/seedDomain", name="traction_x")
      call obj%traction_y%open(path=obj%path//obj%name//"/input/seed/seedDomain", name="traction_y")
      call obj%traction_z%open(path=obj%path//obj%name//"/input/seed/seedDomain", name="traction_z")
      call obj%flux%open(path=obj%path//obj%name//"/input/seed/seedDomain", name="flux")

      ! link
      call obj%seedDomain%link(Water=obj%water, Tissue=obj%Tissue, &
                               a_Psi=obj%a_Psi, &
                               a_P=obj%a_P, &
                               theta_eq=obj%theta_eq, &
                               Psi_eq=obj%Psi_eq, &
                               a_E=obj%a_E, &
                               a_v=obj%a_v, &
                               E_eq=obj%E_eq, &
                               v_eq=obj%v_eq)
      call obj%water%import(Boundaries=.true., Boundary=obj%const)
      call obj%Tissue%import(Boundaries=.true., Boundary=obj%disp_x)
      call obj%Tissue%import(Boundaries=.true., Boundary=obj%disp_y)
      call obj%Tissue%import(Boundaries=.true., Boundary=obj%disp_z)

      call obj%tissue%bake(template="FiniteDeform_")
      call obj%water%bake(template="DiffusionEq_")

      call obj%seedDomain%import(Water=obj%water, Tissue=obj%tissue)
      call obj%seedDomain%bake()

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

!########################################################
   subroutine flushSeed(obj, path, name)
      class(Seed_), intent(inout)::obj
      character(*), optional, intent(in) :: path
      character(*), optional, intent(in) :: name
      type(IO_) :: f

      obj%step = obj%step + 1

      if (present(path)) then
         obj%path = path
      else
         obj%path = "home/"
      end if

      if (present(name)) then
         obj%name = name
      else
         obj%name = "untitled"
      end if

      ! create directory for save
      call execute_command_line("mkdir -p "//obj%path//obj%name)

      call f%open(obj%path//obj%name//"/type")
      call f%write("seed")
      call f%close()

      call f%open(obj%path//obj%name//"/readme")
      call f%write("./input/ : Data stack for plant-FEM analysis")
      call f%write("           ,which contains mesh and its properties.")
      call f%write("./output/ : Data stack for results of plant-FEM analysis")
      call f%write("           ,which contains field informations.")
      call f%write("./server/ : Server script and options for plant-FEM analysis")
      call f%write("           ,which contains .f90 scripts and options.")
      call f%close()

      call execute_command_line("mkdir -p "//obj%path//obj%name//"/input")
      call execute_command_line("mkdir -p "//obj%path//obj%name//"/output")
      call execute_command_line("mkdir -p "//obj%path//obj%name//"/server")

      ! only for SeedClass
      call execute_command_line("mkdir -p "//obj%path//obj%name//"/output/seed"//str(obj%step)// &
                                "")
      call execute_command_line("mkdir -p "//obj%path//obj%name//"/output/seed"//str(obj%step)// &
                                "")
      call execute_command_line("mkdir -p "//obj%path//obj%name//"/output/seed"//str(obj%step)// &
                                "")
      call obj%water%save(path=obj%path//obj%name//"/output/seed"//str(obj%step)// &
                          "", name="water")
      call obj%tissue%save(path=obj%path//obj%name//"/output/seed"//str(obj%step)// &
                           "", name="tissue")
      call obj%seedDomain%save(path=obj%path//obj%name//"/output/seed"//str(obj%step)// &
                               "", name="seedDomain")

      call obj%Permiability%save(path=obj%path//obj%name//"/output/seed"//str(obj%step)// &
                                 "/seedDomain", name="Permiability")
      call obj%YoungsModulus%save(path=obj%path//obj%name//"/output/seed"//str(obj%step)// &
                                  "/seedDomain", name="YoungsModulus")
      call obj%PoissonRatio%save(path=obj%path//obj%name//"/output/seed"//str(obj%step)// &
                                 "/seedDomain", name="PoissonRatio")
      call obj%density%save(path=obj%path//obj%name//"/output/seed"//str(obj%step)// &
                            "/seedDomain", name="density")
      call obj%Cohesion%save(path=obj%path//obj%name//"/output/seed"//str(obj%step)// &
                             "/seedDomain", name="Cohesion")
      call obj%phi%save(path=obj%path//obj%name//"/output/seed"//str(obj%step)// &
                        "/seedDomain", name="phi")
      call obj%psi%save(path=obj%path//obj%name//"/output/seed"//str(obj%step)// &
                        "/seedDomain", name="phi")
      call obj%a_Psi%save(path=obj%path//obj%name//"/output/seed"//str(obj%step)// &
                          "/seedDomain", name="a_Psi")
      call obj%a_P%save(path=obj%path//obj%name//"/output/seed"//str(obj%step)// &
                        "/seedDomain", name="a_P")
      call obj%theta_eq%save(path=obj%path//obj%name//"/output/seed"//str(obj%step)// &
                             "/seedDomain", name="theta_eq")
      call obj%Psi_eq%save(path=obj%path//obj%name//"/output/seed"//str(obj%step)// &
                           "/seedDomain", name="Psi_eq")
      call obj%a_E%save(path=obj%path//obj%name//"/output/seed"//str(obj%step)// &
                        "/seedDomain", name="a_E")
      call obj%a_v%save(path=obj%path//obj%name//"/output/seed"//str(obj%step)// &
                        "/seedDomain", name="a_v")
      call obj%E_eq%save(path=obj%path//obj%name//"/output/seed"//str(obj%step)// &
                         "/seedDomain", name="E_eq")
      call obj%v_eq%save(path=obj%path//obj%name//"/output/seed"//str(obj%step)// &
                         "/seedDomain", name="v_eq")

      call obj%disp_x%save(path=obj%path//obj%name//"/output/seed"//str(obj%step)// &
                           "/seedDomain", name="disp_x")
      call obj%disp_y%save(path=obj%path//obj%name//"/output/seed"//str(obj%step)// &
                           "/seedDomain", name="disp_y")
      call obj%disp_z%save(path=obj%path//obj%name//"/output/seed"//str(obj%step)// &
                           "/seedDomain", name="disp_z")
      call obj%traction_x%save(path=obj%path//obj%name//"/output/seed"//str(obj%step)// &
                               "/seedDomain", name="traction_x")
      call obj%traction_y%save(path=obj%path//obj%name//"/output/seed"//str(obj%step)// &
                               "/seedDomain", name="traction_y")
      call obj%traction_z%save(path=obj%path//obj%name//"/output/seed"//str(obj%step)// &
                               "/seedDomain", name="traction_z")
      call obj%flux%save(path=obj%path//obj%name//"/output/seed"//str(obj%step)// &
                         "/seedDomain", name="flux")
      call obj%const%save(path=obj%path//obj%name//"/output/seed"//str(obj%step)// &
                          "/seedDomain", name="const")

      call f%open(obj%path//obj%name//"/output/seed"//str(obj%step)// &
                  "/seedDomain/", "seed", ".prop")
      write (f%fh, *) obj%num_of_seed ! num of seed
      write (f%fh, *) obj%mass ! seed mass g/cm^3
      write (f%fh, *) obj%water_content ! seed water_content %
      write (f%fh, *) obj%radius ! seed radius (cm)
      write (f%fh, *) obj%width1
      write (f%fh, *) obj%width2
      write (f%fh, *) obj%width3
      write (f%fh, *) obj%width1_origin
      write (f%fh, *) obj%width2_origin
      write (f%fh, *) obj%width3_origin
      write (f%fh, *) obj%location(1) ! seed location (x,y,z)
      write (f%fh, *) obj%location(2) ! seed location (x,y,z)
      write (f%fh, *) obj%location(3) ! seed location (x,y,z)
      write (f%fh, *) obj%path
      write (f%fh, *) obj%name
      call f%close()

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

!########################################################
   subroutine saveSeed(obj, path, name)
      class(Seed_), intent(inout)::obj
      character(*), optional, intent(in) :: path
      character(*), optional, intent(in) :: name
      type(IO_) :: f

      if (present(path)) then
         obj%path = path
      else
         obj%path = "home/"
      end if

      if (present(name)) then
         obj%name = name
      else
         obj%name = "untitled"
      end if

      ! create directory for save
      call execute_command_line("mkdir -p "//obj%path//obj%name)

      call f%open(obj%path//obj%name//"/type")
      call f%write("seed")
      call f%close()

      call f%open(obj%path//obj%name//"/readme")
      call f%write("./input/ : Data stack for plant-FEM analysis")
      call f%write("           ,which contains mesh and its properties.")
      call f%write("./output/ : Data stack for results of plant-FEM analysis")
      call f%write("           ,which contains field informations.")
      call f%write("./server/ : Server script and options for plant-FEM analysis")
      call f%write("           ,which contains .f90 scripts and options.")
      call f%close()

      call execute_command_line("mkdir -p "//obj%path//obj%name//"/input")
      call execute_command_line("mkdir -p "//obj%path//obj%name//"/output")
      call execute_command_line("mkdir -p "//obj%path//obj%name//"/server")

      ! only for SeedClass
      call execute_command_line("mkdir -p "//obj%path//obj%name//"/input/seed")
      call execute_command_line("mkdir -p "//obj%path//obj%name//"/input/seed/water")
      call execute_command_line("mkdir -p "//obj%path//obj%name//"/input/seed/tissue")
      call execute_command_line("mkdir -p "//obj%path//obj%name//"/input/seed/seedDomain")

      call obj%water%save(path=obj%path//obj%name//"/input/seed", name="water")
      call obj%tissue%save(path=obj%path//obj%name//"/input/seed", name="tissue")
      call obj%seedDomain%save(path=obj%path//obj%name//"/input/seed", name="seedDomain")

      call obj%Permiability%save(path=obj%path//obj%name//"/input/seed/seedDomain", name="Permiability")
      call obj%YoungsModulus%save(path=obj%path//obj%name//"/input/seed/seedDomain", name="YoungsModulus")
      call obj%PoissonRatio%save(path=obj%path//obj%name//"/input/seed/seedDomain", name="PoissonRatio")
      call obj%density%save(path=obj%path//obj%name//"/input/seed/seedDomain", name="density")
      call obj%Cohesion%save(path=obj%path//obj%name//"/input/seed/seedDomain", name="Cohesion")
      call obj%phi%save(path=obj%path//obj%name//"/input/seed/seedDomain", name="phi")
      call obj%psi%save(path=obj%path//obj%name//"/input/seed/seedDomain", name="phi")
      call obj%a_Psi%save(path=obj%path//obj%name//"/input/seed/seedDomain", name="a_Psi")
      call obj%a_P%save(path=obj%path//obj%name//"/input/seed/seedDomain", name="a_P")
      call obj%theta_eq%save(path=obj%path//obj%name//"/input/seed/seedDomain", name="theta_eq")
      call obj%Psi_eq%save(path=obj%path//obj%name//"/input/seed/seedDomain", name="Psi_eq")
      call obj%a_E%save(path=obj%path//obj%name//"/input/seed/seedDomain", name="a_E")
      call obj%a_v%save(path=obj%path//obj%name//"/input/seed/seedDomain", name="a_v")
      call obj%E_eq%save(path=obj%path//obj%name//"/input/seed/seedDomain", name="E_eq")
      call obj%v_eq%save(path=obj%path//obj%name//"/input/seed/seedDomain", name="v_eq")

      call obj%disp_x%save(path=obj%path//obj%name//"/input/seed/seedDomain", name="disp_x")
      call obj%disp_y%save(path=obj%path//obj%name//"/input/seed/seedDomain", name="disp_y")
      call obj%disp_z%save(path=obj%path//obj%name//"/input/seed/seedDomain", name="disp_z")
      call obj%traction_x%save(path=obj%path//obj%name//"/input/seed/seedDomain", name="traction_x")
      call obj%traction_y%save(path=obj%path//obj%name//"/input/seed/seedDomain", name="traction_y")
      call obj%traction_z%save(path=obj%path//obj%name//"/input/seed/seedDomain", name="traction_z")
      call obj%flux%save(path=obj%path//obj%name//"/input/seed/seedDomain", name="flux")
      call obj%const%save(path=obj%path//obj%name//"/input/seed/seedDomain", name="const")

      call f%open(obj%path//obj%name//"/input/seed/seedDomain/", "seed", ".prop")
      write (f%fh, *) obj%num_of_seed ! num of seed
      write (f%fh, *) obj%mass ! seed mass g/cm^3
      write (f%fh, *) obj%water_content ! seed water_content %
      write (f%fh, *) obj%radius ! seed radius (cm)
      write (f%fh, *) obj%width1
      write (f%fh, *) obj%width2
      write (f%fh, *) obj%width3
      write (f%fh, *) obj%width1_origin
      write (f%fh, *) obj%width2_origin
      write (f%fh, *) obj%width3_origin
      write (f%fh, *) obj%location(1) ! seed location (x,y,z)
      write (f%fh, *) obj%location(2) ! seed location (x,y,z)
      write (f%fh, *) obj%location(3) ! seed location (x,y,z)
      write (f%fh, *) obj%path
      write (f%fh, *) obj%name
      call f%close()

      obj%step = obj%step + 1
      call obj%seedDomain%result(path=obj%path//obj%name//"/output/seed", name="seedDomain", step=obj%step)
      call execute_command_line("mv *.pos "//obj%path//obj%name//"/output/")
      call execute_command_line("mv *.txt "//obj%path//obj%name//"/output/")
      call execute_command_line("mv *.scf "//obj%path//obj%name//"/output/")
      call execute_command_line("mv *.msh "//obj%path//obj%name//"/output/")
      call execute_command_line("mv ./a.out "//obj%path//obj%name//"/server/")

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

!########################################################
   subroutine createSeed(obj, Name, MeshType, x_num, y_num, z_num, x_len, y_len, z_len, YoungsModulus, PoissonRatio, Permiability, &
                         a_Psi, a_P, theta_eq, Psi_eq, a_E, a_v, E_eq, v_eq)

      class(Seed_), intent(inout) :: obj

      real(real64), optional, intent(in) :: Permiability, YoungsModulus, PoissonRatio

      real(real64), optional, intent(in) :: a_Psi, a_P, theta_eq, Psi_eq, &
                                            a_E, a_v, E_eq, v_eq
      character(*), optional, intent(in) :: MeshType, Name
      character(:), allocatable :: Mesh_Type, fname
      integer(int32), optional, intent(in) :: x_num, y_num, z_num
      integer(int32) :: xnum, ynum, znum
      real(real64), optional, intent(in) :: x_len, y_len, z_len
      real(real64) :: xlen, ylen, zlen
      real(real64) :: val
      xnum = input(default=10, option=x_num)
      ynum = input(default=10, option=y_num)
      znum = input(default=10, option=z_num)

      xlen = input(default=90.0d0, option=x_len)
      ylen = input(default=80.0d0, option=y_len)
      zlen = input(default=70.0d0, option=z_len)

      !  create mesh
      if (present(MeshType)) then
         Mesh_Type = MeshType
      else
         Mesh_Type = "Sphere3D"
      end if

      if (present(Name)) then
         fname = name//"water"
      else
         fname = "water"
      end if
      call obj%water%create(Name=fname, MeshType=Mesh_Type, x_num=xnum, y_num=ynum, x_len=xlen, y_len=ylen, &
                            thickness=zlen, division=znum)
      call obj%tissue%copy(obj%water, onlyMesh=.true.)

      if (present(Name)) then
         fname = name//"tissue"
      else
         fname = "tissue"
      end if
      call obj%tissue%rename(fname)

      ! create material
      ! for deformation analysis
      val = input(default=100.0d0, option=YoungsModulus)
      call obj%YoungsModulus%create(Name="YoungsModulus", ParaValue=val, Layer=1)
      !call YoungsModulus%create(Name="YoungsModulus",x_max=1500.0d0,x_min=0.0d0,y_max=70.0d0,y_min=0.0d0,&
      !z_max=70.0d0,z_min=0.0d0,ParaValue=6000.0d0,Layer=1)
      val = input(default=0.30d0, option=PoissonRatio)
      call obj%PoissonRatio%create(Name="PoissonRatio", ParaValue=val, Layer=2)
      !call PoissonRatio%create(Name="PoissonRatio",x_max=100.0d0,x_min=0.0d0,y_max=50.0d0,y_min=0.0d0,&
      !z_max=50.0d0,z_min=0.0d0,ParaValue=0.10d0,Layer=2)
      call obj%Density%create(Name="Density", ParaValue=0.00d0, Layer=3)
      call obj%Cohesion%create(Name="Cohesion", ParaValue=100000000.0d0, Layer=4)
      call obj%Phi%create(Name="Psi", ParaValue=0.0d0, Layer=5)
      call obj%psi%create(Name="psi", ParaValue=0.0d0, Layer=6)

      ! for diffusion analysis
      val = input(default=1.0d0, option=Permiability)
      call obj%Permiability%create(Name="Permiability", ParaValue=val, Layer=1)

      val = input(default=20.0d0, option=a_Psi)
      call obj%a_Psi%create(Name="a_Psi", Paravalue=val, Layer=1)
      val = input(default=40.0d0, option=a_P)
      call obj%a_P%create(Name="a_P", Paravalue=val, Layer=1)
      val = input(default=1.0d0, option=theta_eq)
      call obj%theta_eq%create(Name="theta_eq", Paravalue=val, Layer=1)
      val = input(default=20.0d0, option=Psi_eq)
      call obj%Psi_eq%create(Name="Psi_eq", Paravalue=val, Layer=1)

      val = input(default=0.0d0, option=a_E)
      call obj%a_E%create(Name="a_E", Paravalue=val, Layer=1)
      val = input(default=0.0d0, option=a_v)
      call obj%a_v%create(Name="a_v", Paravalue=val, Layer=1)
      val = input(default=100.0d0, option=E_eq)
      call obj%E_eq%create(Name="E_eq", Paravalue=val, Layer=1)
      val = input(default=0.30d0, option=v_eq)
      call obj%v_eq%create(Name="v_eq", Paravalue=val, Layer=1)

      ! import Materials onto FEMDomains
      call obj%tissue%import(Materials=.true., Material=obj%YoungsModulus)
      call obj%tissue%import(Materials=.true., Material=obj%PoissonRatio)
      call obj%tissue%import(Materials=.true., Material=obj%Density)
      call obj%tissue%import(Materials=.true., Material=obj%Cohesion)
      call obj%tissue%import(Materials=.true., Material=obj%Phi)
      call obj%tissue%import(Materials=.true., Material=obj%Psi)

      call obj%water%import(Materials=.true., Material=obj%Permiability)

      call obj%seedDomain%import(a_Psi=obj%a_Psi)
      call obj%seedDomain%import(a_P=obj%a_P)
      call obj%seedDomain%import(theta_eq=obj%theta_eq)
      call obj%seedDomain%import(Psi_eq=obj%Psi_eq)
      call obj%seedDomain%import(a_E=obj%a_E)
      call obj%seedDomain%import(a_v=obj%a_v)
      call obj%seedDomain%import(E_eq=obj%E_eq)
      call obj%seedDomain%import(v_eq=obj%v_eq)

      ! visualize material
      call obj%YoungsModulus%gmsh(Name="YoungsModulus", Tag="M : YoungsModulus (kPa)")
      call obj%PoissonRatio%gmsh(Name="PoissonRatio", Tag="M : PoissonRatio")
      call obj%Permiability%gmsh(Name="Permiability", Tag="M : Permiability (cm/s)")

      call obj%a_Psi%gmsh(Name="a_Psi   ", Tag="a_Psi   ")
      call obj%a_P%gmsh(Name="a_P     ", Tag="a_P     ")
      call obj%theta_eq%gmsh(Name="theta_eq", Tag="theta_eq")
      call obj%Psi_eq%gmsh(Name="Psi_eq  ", Tag="Psi_eq  ")
      call obj%a_E%gmsh(Name="a_E     ", Tag="a_E     ")
      call obj%a_v%gmsh(Name="a_v     ", Tag="a_v     ")
      call obj%E_eq%gmsh(Name="E_eq    ", Tag="E_eq    ")
      call obj%v_eq%gmsh(Name="v_eq    ", Tag="v_eq    ")

      if (present(Name)) then
         fname = name//"water"
      else
         fname = "water"
      end if
      call obj%water%gmsh(Name=fname, Tag="Water Domain")

      if (present(Name)) then
         fname = name//"tissue"
      else
         fname = "tissue"
      end if
      call obj%tissue%gmsh(Name=fname, Tag="Tissue Domain")

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

   subroutine gmshSeed(obj, Name)
      class(Seed_), intent(inout) :: obj
      character(*), optional, intent(in) :: Name
      character(200) :: fname

      if (present(Name)) then
         fname = name//"water"
      else
         fname = "water"
      end if
      call obj%water%gmsh(Name=fname, Tag="Water Domain")

      if (present(Name)) then
         fname = name//"tissue"
      else
         fname = "tissue"
      end if
      call obj%tissue%gmsh(Name=fname, Tag="Tissue Domain")
   end subroutine

!########################################################
   subroutine moveSeed(obj, x, y, z)
      class(Seed_), intent(inout) :: obj
      real(real64), optional, intent(in) :: x, y, z

      call obj%water%move(x=x, y=y, z=z)
      call obj%tissue%move(x=x, y=y, z=z)
   end subroutine
!########################################################

!########################################################
   subroutine rotateSeed(obj, x, y, z)
      class(Seed_), intent(inout) :: obj
      real(real64), optional, intent(in) :: x, y, z

      call obj%water%rotate(x=x, y=y, z=z)
      call obj%tissue%rotate(x=x, y=y, z=z)
   end subroutine
!########################################################

!########################################################
   subroutine materialSeed(obj, YoungsModulus, PoissonRatio, Permiability, &
                           a_Psi, a_P, theta_eq, Psi_eq, a_E, a_v, E_eq, v_eq, x_max, x_min, y_max, y_min, &
                           z_max, z_min)

      class(Seed_), intent(inout) :: obj
      real(real64), optional, intent(in) :: Permiability, YoungsModulus, PoissonRatio

      real(real64), optional, intent(in) :: a_Psi, a_P, theta_eq, Psi_eq, &
                                            a_E, a_v, E_eq, v_eq

      real(real64), optional, intent(in) :: x_max, x_min, y_max, y_min, &
                                            z_max, z_min

      ! create material
      ! for deformation analysis
      if (present(YoungsModulus)) then
         call obj%YoungsModulus%create(Name="YoungsModulus", x_max=x_max, x_min=x_min, y_max=y_max, y_min=y_min, &
                                       z_max=z_max, z_min=z_min, ParaValue=YoungsModulus, Layer=1)
      end if
      if (present(PoissonRatio)) then
         call obj%PoissonRatio%create(Name="PoissonRatio ", x_max=x_max, x_min=x_min, y_max=y_max, y_min=y_min, &
                                      z_max=z_max, z_min=z_min, ParaValue=PoissonRatio, Layer=2)
      end if
      if (present(Permiability)) then
         call obj%Permiability%create(Name="Permiability", x_max=x_max, x_min=x_min, y_max=y_max, y_min=y_min, &
                                      z_max=z_max, z_min=z_min, ParaValue=Permiability, Layer=1)
      end if

      if (present(a_Psi)) then
         call obj%a_Psi%create(Name="a_Psi", x_max=x_max, x_min=x_min, y_max=y_max, y_min=y_min, &
                               z_max=z_max, z_min=z_min, ParaValue=a_Psi, Layer=1)
      end if

      if (present(a_P)) then
         call obj%a_P%create(Name="a_P", x_max=x_max, x_min=x_min, y_max=y_max, y_min=y_min, &
                             z_max=z_max, z_min=z_min, ParaValue=a_P, Layer=1)
      end if

      if (present(theta_eq)) then
         call obj%theta_eq%create(Name="theta_eq", x_max=x_max, x_min=x_min, y_max=y_max, y_min=y_min, &
                                  z_max=z_max, z_min=z_min, ParaValue=theta_eq, Layer=1)
      end if

      if (present(Psi_eq)) then
         call obj%Psi_eq%create(Name="Psi_eq", x_max=x_max, x_min=x_min, y_max=y_max, y_min=y_min, &
                                z_max=z_max, z_min=z_min, ParaValue=Psi_eq, Layer=1)
      end if

      if (present(a_E)) then
         call obj%a_E%create(Name="a_E", x_max=x_max, x_min=x_min, y_max=y_max, y_min=y_min, &
                             z_max=z_max, z_min=z_min, ParaValue=a_E, Layer=1)
      end if

      if (present(a_v)) then
         call obj%a_v%create(Name="a_v", x_max=x_max, x_min=x_min, y_max=y_max, y_min=y_min, &
                             z_max=z_max, z_min=z_min, ParaValue=a_v, Layer=1)
      end if

      if (present(E_eq)) then
         call obj%E_eq%create(Name="E_eq", x_max=x_max, x_min=x_min, y_max=y_max, y_min=y_min, &
                              z_max=z_max, z_min=z_min, ParaValue=E_eq, Layer=1)
      end if

      if (present(v_eq)) then
         call obj%v_eq%create(Name="v_eq", x_max=x_max, x_min=x_min, y_max=y_max, y_min=y_min, &
                              z_max=z_max, z_min=z_min, ParaValue=v_eq, Layer=1)
      end if

      !call obj%Density%create(Name="Density",ParaValue=0.00d0,Layer=3)
      !call obj%Cohesion%create(Name="Cohesion",ParaValue=100000000.0d0,Layer=4)
      !call obj%Phi%create(Name="Psi",ParaValue=0.0d0,Layer=5)
      !call obj%psi%create(Name="psi",ParaValue=0.0d0,Layer=6)

   end subroutine

!########################################################
   subroutine envSeed(obj, disp_x, disp_y, disp_z, WaterContent, x_max, x_min, y_max, y_min, &
                      z_max, z_min)
      class(Seed_), intent(inout) :: obj
      real(real64), optional, intent(in)::disp_x, disp_y, disp_z, WaterContent
      real(real64), optional, intent(in) :: x_max, x_min, y_max, y_min, &
                                            z_max, z_min

      ! please use Layer (1, 2, 3...).
      if (present(disp_x)) then
         call obj%disp_x%create(Category="Dirichlet", x_max=x_max, x_min=x_min, y_max=y_max, y_min=y_min, &
                                z_max=z_max, z_min=z_min, BoundValue=disp_x, Layer=1, Name="disp_x")
         call obj%tissue%import(Boundaries=.true., Boundary=obj%disp_x)
      end if

      if (present(disp_y)) then
         call obj%disp_y%create(Category="Dirichlet", x_max=x_max, x_min=x_min, y_max=y_max, y_min=y_min, &
                                z_max=z_max, z_min=z_min, BoundValue=disp_y, Layer=2, Name="disp_y")
         call obj%tissue%import(Boundaries=.true., Boundary=obj%disp_y)
      end if

      if (present(disp_z)) then
         call obj%disp_z%create(Category="Dirichlet", x_max=x_max, x_min=x_min, y_max=y_max, y_min=y_min, &
                                z_max=z_max, z_min=z_min, BoundValue=disp_z, Layer=3, Name="disp_z")
         call obj%tissue%import(Boundaries=.true., Boundary=obj%disp_z)

      end if

      if (present(WaterContent)) then
         call obj%const%create(Category="Dirichlet", x_max=x_max, x_min=x_min, y_max=y_max, y_min=y_min, &
                               z_max=z_max, z_min=z_min, BoundValue=WaterContent, Layer=1, Name="const")
         call obj%water%import(Boundaries=.true., Boundary=obj%const)

      end if

   end subroutine

!########################################################

!########################################################
   subroutine growSeed(obj, timestep, dt, Display, nr_tol, interval, Name, restart)

      class(Seed_), intent(inout) :: obj
      integer(int32), optional, intent(in) :: timestep, interval
      real(real64), optional, intent(in) :: dt, nr_tol
      logical, optional, intent(in) ::Display, restart
      character(*), optional, intent(in) :: Name
      ! bake data by using templates

      if (present(restart)) then
         if (restart .eqv. .true.) then
            call obj%seedDomain%run(timestep=timestep, dt=dt, SolverType="BiCGSTAB", &
                                  Display=.true., nr_tol=nr_tol, infinitesimal=.true., interval=interval, Name=Name, restart=.true.)
            return
         end if
      end if

      call obj%tissue%bake(template="FiniteDeform_")
      call obj%water%bake(template="DiffusionEq_")

      ! visualize Domains
      call obj%tissue%gmsh(Name="tissue", Tag="Tissue Domain")
      call obj%water%gmsh(Name="water", Tag="Water Domain")

      call obj%tissue%export(Name="tissue")
      call obj%water%export(Name="water")

      ! mount data into a seed (An instance of Water-Absorption Solver)

      call obj%seedDomain%import(Water=obj%water, Tissue=obj%tissue)
      call obj%seedDomain%bake()
      ! run simula

      call obj%seedDomain%run(timestep=timestep, dt=dt, SolverType="BiCGSTAB", &
                              Display=.true., nr_tol=nr_tol, infinitesimal=.true., interval=interval, Name=Name)
      call obj%gmsh()
   end subroutine
!########################################################

!########################################################
   subroutine initSeed(obj, mass, water_content, radius, location, x, y, z, width1, width2, width3)
      class(Seed_), intent(inout) :: obj
      real(real64), optional, intent(in) :: mass, water_content, radius, location(3), x, y, z
      real(real64), optional, intent(in) :: width1, width2, width3
      real(real64) :: loc(3)
      loc(:) = 0.0d0
      if (present(location)) then
         loc(:) = location(:)
      end if
      obj%num_of_seed = 1
      obj%mass = input(default=1.0d0, option=mass)
      obj%water_content = input(default=12.0d0, option=water_content)

      obj%width1 = input(default=obj%radius, option=width1)
      obj%width2 = input(default=obj%radius, option=width2)
      obj%width3 = input(default=obj%radius, option=width3)
      obj%radius = input(default=0.330d0*obj%width1 + 0.330d0*obj%width2 + 0.330d0*obj%width3, option=radius)

      obj%location(:) = loc(:)
      obj%location(1) = obj%location(1) + input(default=0.0d0, option=x)
      obj%location(2) = obj%location(2) + input(default=0.0d0, option=y)
      obj%location(3) = obj%location(3) + input(default=0.0d0 + obj%radius, option=z)

      call obj%FEMDomain%init(simple=.true.)
      obj%width1_origin = obj%width1
      obj%width2_origin = obj%width2
      obj%width3_origin = obj%width3

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

!########################################################
   subroutine importSeed(obj, FileName)
      class(Seed_), intent(inout) :: obj
      character(*), intent(in) :: FileName

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

!########################################################
   subroutine createMeshSeed(obj, FileName, withSTL, ObjType, ElemType)

      class(Seed_), intent(inout):: obj
      type(Mesh_) :: mesh
      character(*), intent(in)    :: FileName, ElemType
      character(*), optional, intent(in) :: ObjType
      character*200   :: meshFileName, meshFileName_m
      character*200   :: command
      character*200   :: strings
      integer(int32) :: intval
      real(real64) :: x_rate, y_rate, z_rate
      logical, optional, intent(in) :: withSTL

      meshFileName = FileName//"createMesh_seed.geo"
      meshFileName_m = FileName//"createMesh_seed.mesh"
      call obj%export(FileName=meshFileName)
      command = "gmsh "//meshFileName//" -3 -format mesh"
      writE (*, '(A)') command
      call execute_command_line(command)

      ! modification
      call mesh%import(FileName=meshFileName_m, extention=".mesh", ElemType=ElemType)
      call obj%FEMDomain%importMesh(mesh)
      x_rate = obj%width1/(maxval(obj%FEMDomain%Mesh%NodCoord(:, 1)) - minval(obj%FEMDomain%Mesh%NodCoord(:, 1)))
      y_rate = obj%width2/(maxval(obj%FEMDomain%Mesh%NodCoord(:, 2)) - minval(obj%FEMDomain%Mesh%NodCoord(:, 2)))
      z_rate = obj%width3/(maxval(obj%FEMDomain%Mesh%NodCoord(:, 3)) - minval(obj%FEMDomain%Mesh%NodCoord(:, 3)))
      call obj%FEMDomain%resize(x_rate=x_rate, y_rate=y_rate, z_rate=z_rate)

      if (present(withSTL)) then
         if (withSTL .eqv. .true.) then
            call obj%export(FileName=FileName, extention=".stl")
         end if
      end if
   end subroutine
!########################################################

!########################################################
   subroutine exportSeed(obj, FileName, SeedID, extention)
      class(Seed_), intent(inout) :: obj
      character(*), intent(in) :: FileName
      character(*), optional, intent(in) ::extention
      character(10) :: ex_format
      integer(int32), optional, intent(in) :: SeedID

      if (.not. present(extention)) then
         ex_format = ".geo"
      else
         ex_format = extention
      end if

      if (ex_format == ".stl" .or. ex_format == "stl") then

         call obj%FEMDomain%export(OptionalFileFormat="stl", &
                                   FileName=FileName//"seed"//adjustl(fstring(input(default=1, option=SeedID)))//".stl", &
                                   MeshDimension=3)

      elseif (ex_format == ".pos" .or. ex_format == "pos") then
         call obj%FEMDomain%GmshPlotMesh(Name=FileName//"seed"//adjustl(fstring(input(default=1, option=SeedID)))//".pos")
      elseif (ex_format == ".geo" .or. ex_format == "geo") then
         open (10, file=FileName)
         write (10, '(A)') "//+"
         write (10, '(A)') 'SetFactory("OpenCASCADE");'
         write (10, *) "Sphere(", input(default=1, option=SeedID), ") = {", &
            obj%location(1), ",", obj%location(2), ",", obj%location(3), ",", &
            obj%radius, ", -Pi/2, Pi/2, 2*Pi};"
         close (10)
      else
         print *, "ERROR :: export formats .pos, .stl, .geo are available."
      end if

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

   subroutine convertMeshTypeSeed(obj, Option)
      class(Seed_), intent(inout)::obj
      character(*), intent(in) :: Option

      call obj%FEMDomain%convertMeshType(Option=Option)
   end subroutine

!########################################################
   subroutine showSeed(obj)
      class(Seed_), intent(in) :: obj

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

!function create_seed_shapeSeed(this,x_num,y_num,z_num,width,length,thickness,seed_type) result(FEMDomain)
!    class(Seed_),intent(in)   :: this
!    integer(int32),intent(in) :: x_num,y_num,z_num
!    real(real64),intent(in)   :: width,length,thickness,
!    character(*),intent(in)   :: seed_type
!    type(FEMDomain_) :: FEMDomain
!
!    select case(seed_type)
!        case("RICE","rice")
!            call FEMDomain%init("Cube3D")
!
!    end select
!end function

end module SeedClass
