module WaterAbsorptionClass
   use fem
   use DiffusionEquationClass
   use FiniteDeformationClass
   implicit none

   type :: WaterAbsorption_
      type(FEMDomain_), pointer:: Water, Tissue
      type(MaterialProp_), pointer:: a_Psi
      type(MaterialProp_), pointer:: a_P
      type(MaterialProp_), pointer:: theta_eq
      type(MaterialProp_), pointer:: Psi_eq
      type(MaterialProp_), pointer:: a_E
      type(MaterialProp_), pointer:: a_v
      type(MaterialProp_), pointer:: E_eq
      type(MaterialProp_), pointer:: v_eq

      type(DiffusionEq_)::DiffusionEq
      type(FiniteDeform_)::FiniteDeform

      real(real64), allocatable :: WaterAbsorbingPower(:)
      real(real64), allocatable :: WaterPotential(:)
      real(real64), allocatable :: TurgorPressure(:)
      real(real64), allocatable :: WaterContent(:)
      real(real64), allocatable :: WaterMass(:)
      real(real64), allocatable :: volume(:)
      real(real64), allocatable :: Conductivity(:)
      real(real64), allocatable :: YoungsModulus(:)
      real(real64), allocatable :: PoissonsRatio(:)
      real(real64), allocatable :: Permiability(:)
      real(real64), allocatable :: PorePressure(:)
      ! Parameters(A,B)
      ! A : Element ID

      real(real64), allocatable :: a_Psi_val(:)
      real(real64), allocatable :: a_P_val(:)
      real(real64), allocatable :: theta_eq_val(:)
      real(real64), allocatable :: Psi_eq_val(:)
      real(real64), allocatable :: a_E_val(:)
      real(real64), allocatable :: a_v_val(:)
      real(real64), allocatable :: E_eq_val(:)
      real(real64), allocatable :: v_eq_val(:)

      real(real64), allocatable :: theta_ps_val(:) ! computed from other variables

      character(200) :: Name = " "
      integer(int32) :: timestep = 0
      integer(int32) :: flushstep = 0
      integer(int32) :: ID = 0
      real(real64) :: dt
   contains
      procedure, public :: import => importWaterAbsorption
      procedure, public :: assign => importWaterAbsorption
      procedure, public :: init => initWaterAbsorption
      procedure, public :: run => runWaterAbsorption
      procedure, public :: update => updateWaterAbsorption
      procedure, public :: export => exportWaterAbsorption
      procedure, public :: display => displayWaterAbsorption
      procedure, public :: bake => bakeWaterAbsorption
      procedure, public :: gnuplot => gnuplotWaterAbsorption
      procedure, public :: updatePermiability => updatePermiabilityWA
      procedure, public :: updateStiffness => updateStiffnessWA
      procedure, public :: UpdateVolume => UpdateVolumeWA

      ! <<<<<<<< IO >>>>>>>>>>
      procedure, public :: save => saveWaterAbsorption
      procedure, public :: open => openWaterAbsorption
      procedure, public :: link => linkWaterAbsorption
      procedure, public :: result => resultWaterAbsorption
      procedure, public :: remove => removeWaterAbsorption
      ! <<<<<<<< IO >>>>>>>>>>

   end type
contains

!#####################################
   subroutine removeWaterAbsorption(obj)
      class(WaterAbsorption_), intent(inout) :: obj

      if (associated(obj%Water)) nullify (obj%Water)
      if (associated(obj%Tissue)) nullify (obj%Tissue)

      if (associated(obj%a_Psi)) nullify (obj%a_Psi)
      if (associated(obj%a_P)) nullify (obj%a_P)
      if (associated(obj%theta_eq)) nullify (obj%theta_eq)
      if (associated(obj%Psi_eq)) nullify (obj%Psi_eq)
      if (associated(obj%a_E)) nullify (obj%a_E)
      if (associated(obj%a_v)) nullify (obj%a_v)
      if (associated(obj%E_eq)) nullify (obj%E_eq)
      if (associated(obj%v_eq)) nullify (obj%v_eq)

      call obj%DiffusionEq%remove()
      call obj%FiniteDeform%remove()

      if (allocated(obj%WaterAbsorbingPower)) deallocate (obj%WaterAbsorbingPower)
      if (allocated(obj%WaterPotential)) deallocate (obj%WaterPotential)
      if (allocated(obj%TurgorPressure)) deallocate (obj%TurgorPressure)
      if (allocated(obj%WaterContent)) deallocate (obj%WaterContent)
      if (allocated(obj%WaterMass)) deallocate (obj%WaterMass)
      if (allocated(obj%Conductivity)) deallocate (obj%Conductivity)
      if (allocated(obj%YoungsModulus)) deallocate (obj%YoungsModulus)
      if (allocated(obj%PoissonsRatio)) deallocate (obj%PoissonsRatio)
      if (allocated(obj%Permiability)) deallocate (obj%Permiability)
      if (allocated(obj%PorePressure)) deallocate (obj%PorePressure)
      ! Parameters(A,B)
      ! A : Element ID

      if (allocated(obj%a_Psi_val)) deallocate (obj%a_Psi_val)
      if (allocated(obj%a_P_val)) deallocate (obj%a_P_val)
      if (allocated(obj%theta_eq_val)) deallocate (obj%theta_eq_val)
      if (allocated(obj%Psi_eq_val)) deallocate (obj%Psi_eq_val)
      if (allocated(obj%a_E_val)) deallocate (obj%a_E_val)
      if (allocated(obj%a_v_val)) deallocate (obj%a_v_val)
      if (allocated(obj%E_eq_val)) deallocate (obj%E_eq_val)
      if (allocated(obj%v_eq_val)) deallocate (obj%v_eq_val)
      if (allocated(obj%theta_ps_val)) deallocate (obj%theta_ps_val)

      obj%Name = " "
      obj%timestep = 0
      obj%flushstep = 0
      obj%dt = 1.0d0

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

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

      call obj%DiffusionEq%display(OptionalStep=step, Name=path//"/output/"//name)
      call obj%FiniteDeform%display(OptionalStep=step, Name=path//"/output/"//name)

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

!#####################################
   subroutine linkWaterAbsorption(obj, Water, Tissue, a_Psi, a_P, theta_eq, Psi_eq, a_E, a_v, E_eq, v_eq)
      class(WaterAbsorption_), intent(inout) :: obj
      type(FEMDomain_), optional, target :: Water, Tissue
      type(MaterialProp_), optional, target :: a_Psi, a_P, theta_eq, Psi_eq, a_E, a_v, E_eq, v_eq

      if (present(Water)) then
         if (associated(obj%Water)) nullify (obj%Water)
         obj%Water => Water
      end if
      if (present(Tissue)) then
         if (associated(obj%Tissue)) nullify (obj%Tissue)
         obj%Tissue => Tissue
      end if
      if (present(a_Psi)) then
         if (associated(obj%a_Psi)) nullify (obj%a_Psi)
         obj%a_Psi => a_Psi
      end if
      if (present(a_P)) then
         if (associated(obj%a_P)) nullify (obj%a_P)
         obj%a_P => a_P
      end if
      if (present(theta_eq)) then
         if (associated(obj%theta_eq)) nullify (obj%theta_eq)
         obj%theta_eq => theta_eq
      end if
      if (present(Psi_eq)) then
         if (associated(obj%Psi_eq)) nullify (obj%Psi_eq)
         obj%Psi_eq => Psi_eq
      end if
      if (present(a_E)) then
         if (associated(obj%a_E)) nullify (obj%a_E)
         obj%a_E => a_E
      end if
      if (present(a_v)) then
         if (associated(obj%a_v)) nullify (obj%a_v)
         obj%a_v => a_v
      end if
      if (present(E_eq)) then
         if (associated(obj%E_eq)) nullify (obj%E_eq)
         obj%E_eq => E_eq
      end if
      if (present(v_eq)) then
         if (associated(obj%v_eq)) nullify (obj%v_eq)
         obj%v_eq => v_eq
      end if

      if (present(Water)) then
         if (associated(obj%DiffusionEq%FEMDomain)) nullify (obj%DiffusionEq%FEMDomain)
         obj%DiffusionEq%FEMDomain => Water
      end if

      if (present(Tissue)) then
         if (associated(obj%FiniteDeform%FEMDomain)) nullify (obj%FiniteDeform%FEMDomain)
         obj%FiniteDeform%FEMDomain => Tissue
      end if
   end subroutine
!#####################################

!#####################################
   subroutine openWaterAbsorption(obj, path, name)
      class(WaterAbsorption_), intent(inout) :: obj
      character(*), optional, intent(in)::path
      character(*), optional, intent(in) :: Name
      character(200) ::fname
      integer(int32) :: fh
      logical :: restart = .true.
      type(IO_) :: f

      if (.not. present(path)) then
         print *, " exportWaterAbsorption ERROR >> .not. present(path)"
         stop
      end if

      if (present(name)) then
         call execute_command_line("mkdir -p "//path//"/"//name)
         fname = path//"/"//name
         call f%open("./", path//"/"//name, "/WaterAbsorption.prop")
         call obj%FiniteDeform%open(path=path//"/"//name, name="FiniteDeform")
         call obj%DiffusionEq%open(path=path//"/"//name, name="DiffusionEq")
      else
         call execute_command_line("mkdir -p "//path//"/WaterAbsorption")
         call f%open("./", path//"/WaterAbsorption", "/WaterAbsorption.prop")
         call obj%FiniteDeform%open(path=path//"/WaterAbsorption", name="FiniteDeform")
         call obj%DiffusionEq%open(path=path//"/WaterAbsorption", name="DiffusionEq")
         fname = path//"/WaterAbsorption"
      end if

      call openArray(f%fh, obj%WaterAbsorbingPower)
      call openArray(f%fh, obj%WaterPotential)
      call openArray(f%fh, obj%TurgorPressure)
      call openArray(f%fh, obj%WaterContent)
      call openArray(f%fh, obj%WaterMass)
      call openArray(f%fh, obj%Conductivity)
      call openArray(f%fh, obj%YoungsModulus)
      call openArray(f%fh, obj%PoissonsRatio)
      call openArray(f%fh, obj%Permiability)
      call openArray(f%fh, obj%PorePressure)
      call openArray(f%fh, obj%a_Psi_val)
      call openArray(f%fh, obj%a_P_val)
      call openArray(f%fh, obj%theta_eq_val)
      call openArray(f%fh, obj%Psi_eq_val)
      call openArray(f%fh, obj%a_E_val)
      call openArray(f%fh, obj%a_v_val)
      call openArray(f%fh, obj%E_eq_val)
      call openArray(f%fh, obj%v_eq_val)
      call openArray(f%fh, obj%theta_ps_val) ! computed from other )
      read (f%fh, '(A)') obj%Name
      read (f%fh, *) obj%timestep
      read (f%fh, *) obj%dt
      call f%close()

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

!#####################################
   subroutine flushWaterAbsorption(obj, path, name)
      class(WaterAbsorption_), intent(inout) :: obj
      character(*), optional, intent(in)::path
      character(*), optional, intent(in) :: Name
      character(200) ::fname
      character(200) ::fs
      integer(int32) :: fh
      logical :: restart = .true.
      type(IO_) :: f

      fs = str(obj%flushstep)
      if (.not. present(path)) then
         print *, " exportWaterAbsorption ERROR >> .not. present(path)"
         stop
      end if

      if (present(name)) then
         fname = path//"/"//name
         call execute_command_line("mkdir -p "//path//"/"//name)

         if (associated(obj%water)) then
            call obj%Water%save(path="./"//fname//"/water", name="water")
         end if

         if (associated(obj%tissue)) then
            call obj%Tissue%save(path="./"//fname//"/", name="tissue")
         end if

         call obj%DiffusionEq%export(path="./"//fname//"/DiffusionEq", restart=restart)
         call obj%FiniteDeform%export(path="./"//fname//"/FiniteDeform", restart=restart)

         if (associated(obj%a_Psi)) then
            call execute_command_line("mkdir -p ./"//fname//"/a_Psi")
            call obj%a_Psi%export(path=fname//"/a_Psi", restart=restart)
         end if

         if (associated(obj%a_P)) then
            call execute_command_line("mkdir -p ./"//fname//"/a_P")
            call obj%a_P%export(path=fname//"/a_P", restart=restart)
         end if

         if (associated(obj%theta_eq)) then
            call execute_command_line("mkdir -p ./"//fname//"/theta_eq")
            call obj%theta_eq%export(path=fname//"/theta_eq", restart=restart)
         end if

         if (associated(obj%Psi_eq)) then
            call execute_command_line("mkdir -p ./"//fname//"/Psi_eq")
            call obj%Psi_eq%export(path=fname//"/Psi_eq", restart=restart)
         end if

         if (associated(obj%a_E)) then
            call execute_command_line("mkdir -p ./"//fname//"/a_E")
            call obj%a_E%export(path=fname//"/a_E", restart=restart)
         end if

         if (associated(obj%a_v)) then
            call execute_command_line("mkdir -p ./"//fname//"/a_v")
            call obj%a_v%export(path=fname//"/a_v", restart=restart)
         end if

         if (associated(obj%E_eq)) then
            call execute_command_line("mkdir -p ./"//fname//"/E_eq")
            call obj%E_eq%export(path=fname//"/E_eq", restart=restart)
         end if

         if (associated(obj%v_eq)) then
            call execute_command_line("mkdir -p ./"//fname//"/v_eq")
            call obj%v_eq%export(path=fname//"/v_eq", restart=restart)
         end if

         call f%open("./", path//"/"//name, "/WaterAbsorption.prop")

         write (f%fh, *) obj%WaterAbsorbingPower(:)
         write (f%fh, *) obj%WaterPotential(:)
         write (f%fh, *) obj%TurgorPressure(:)
         write (f%fh, *) obj%WaterContent(:)
         write (f%fh, *) obj%WaterMass(:)
         write (f%fh, *) obj%Conductivity(:)
         write (f%fh, *) obj%YoungsModulus(:)
         write (f%fh, *) obj%PoissonsRatio(:)
         write (f%fh, *) obj%Permiability(:)
         write (f%fh, *) obj%PorePressure(:)
         write (f%fh, *) obj%a_Psi_val(:)
         write (f%fh, *) obj%a_P_val(:)
         write (f%fh, *) obj%theta_eq_val(:)
         write (f%fh, *) obj%Psi_eq_val(:)
         write (f%fh, *) obj%a_E_val(:)
         write (f%fh, *) obj%a_v_val(:)
         write (f%fh, *) obj%E_eq_val(:)
         write (f%fh, *) obj%v_eq_val(:)
         write (f%fh, *) obj%theta_ps_val(:) ! computed from other variables
         write (f%fh, *) obj%Name
         write (f%fh, *) obj%timestep
         write (f%fh, *) obj%dt
         call f%close()
      else
         call execute_command_line("mkdir -p "//path//"/WaterAbsorption")
         call f%open("./", path//"/WaterAbsorption", "/WaterAbsorption.prop")
         fname = path//"/WaterAbsorption"
         write (f%fh, *) obj%WaterAbsorbingPower(:)
         write (f%fh, *) obj%WaterPotential(:)
         write (f%fh, *) obj%TurgorPressure(:)
         write (f%fh, *) obj%WaterContent(:)
         write (f%fh, *) obj%WaterMass(:)
         write (f%fh, *) obj%Conductivity(:)
         write (f%fh, *) obj%YoungsModulus(:)
         write (f%fh, *) obj%PoissonsRatio(:)
         write (f%fh, *) obj%Permiability(:)
         write (f%fh, *) obj%PorePressure(:)
         write (f%fh, *) obj%a_Psi_val(:)
         write (f%fh, *) obj%a_P_val(:)
         write (f%fh, *) obj%theta_eq_val(:)
         write (f%fh, *) obj%Psi_eq_val(:)
         write (f%fh, *) obj%a_E_val(:)
         write (f%fh, *) obj%a_v_val(:)
         write (f%fh, *) obj%E_eq_val(:)
         write (f%fh, *) obj%v_eq_val(:)
         write (f%fh, *) obj%theta_ps_val(:) ! computed from other variables
         write (f%fh, *) obj%Name
         write (f%fh, *) obj%timestep
         write (f%fh, *) obj%dt
         call f%close()

         if (associated(obj%water)) then
            call obj%Water%save(path="./"//fname//"/water", name="water")
         end if

         if (associated(obj%tissue)) then
            call obj%Tissue%save(path="./"//fname//"/", name="tissue")
         end if

         call obj%DiffusionEq%export(path="./"//fname//"/DiffusionEq", restart=restart)
         call obj%FiniteDeform%export(path="./"//fname//"/FiniteDeform", restart=restart)

         if (associated(obj%a_Psi)) then
            call execute_command_line("mkdir -p ./"//fname//"/a_Psi")
            call obj%a_Psi%export(path=fname//"/a_Psi", restart=restart)
         end if

         if (associated(obj%a_P)) then
            call execute_command_line("mkdir -p ./"//fname//"/a_P")
            call obj%a_P%export(path=fname//"/a_P", restart=restart)
         end if

         if (associated(obj%theta_eq)) then
            call execute_command_line("mkdir -p ./"//fname//"/theta_eq")
            call obj%theta_eq%export(path=fname//"/theta_eq", restart=restart)
         end if

         if (associated(obj%Psi_eq)) then
            call execute_command_line("mkdir -p ./"//fname//"/Psi_eq")
            call obj%Psi_eq%export(path=fname//"/Psi_eq", restart=restart)
         end if

         if (associated(obj%a_E)) then
            call execute_command_line("mkdir -p ./"//fname//"/a_E")
            call obj%a_E%export(path=fname//"/a_E", restart=restart)
         end if

         if (associated(obj%a_v)) then
            call execute_command_line("mkdir -p ./"//fname//"/a_v")
            call obj%a_v%export(path=fname//"/a_v", restart=restart)
         end if

         if (associated(obj%E_eq)) then
            call execute_command_line("mkdir -p ./"//fname//"/E_eq")
            call obj%E_eq%export(path=fname//"/E_eq", restart=restart)
         end if

         if (associated(obj%v_eq)) then
            call execute_command_line("mkdir -p ./"//fname//"/v_eq")
            call obj%v_eq%export(path=fname//"/v_eq", restart=restart)
         end if

      end if

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

!#####################################
   subroutine saveWaterAbsorption(obj, path, name)
      class(WaterAbsorption_), intent(inout) :: obj
      character(*), optional, intent(in)::path
      character(*), optional, intent(in) :: Name
      character(200) ::fname
      integer(int32) :: fh
      logical :: restart = .true.
      type(IO_) :: f

      if (.not. present(path)) then
         print *, " exportWaterAbsorption ERROR >> .not. present(path)"
         stop
      end if

      if (present(name)) then
         call execute_command_line("mkdir -p "//path//"/"//name)
         fname = path//"/"//name
         call f%open("./", path//"/"//name, "/WaterAbsorption.prop")
         call obj%FiniteDeform%save(path=path//"/"//name, name="FiniteDeform")
         call obj%DiffusionEq%save(path=path//"/"//name, name="DiffusionEq")

      else
         call execute_command_line("mkdir -p "//path//"/WaterAbsorption")
         call f%open("./", path//"/WaterAbsorption", "/WaterAbsorption.prop")
         call obj%FiniteDeform%save(path=path//"/WaterAbsorption", name="FiniteDeform")
         call obj%DiffusionEq%save(path=path//"/WaterAbsorption", name="DiffusionEq")
         fname = path//"/WaterAbsorption"
      end if

      call WriteArray(f%fh, obj%WaterAbsorbingPower)
      call WriteArray(f%fh, obj%WaterPotential)
      call WriteArray(f%fh, obj%TurgorPressure)
      call WriteArray(f%fh, obj%WaterContent)
      call WriteArray(f%fh, obj%WaterMass)
      call WriteArray(f%fh, obj%Conductivity)
      call WriteArray(f%fh, obj%YoungsModulus)
      call WriteArray(f%fh, obj%PoissonsRatio)
      call WriteArray(f%fh, obj%Permiability)
      call WriteArray(f%fh, obj%PorePressure)
      call WriteArray(f%fh, obj%a_Psi_val)
      call WriteArray(f%fh, obj%a_P_val)
      call WriteArray(f%fh, obj%theta_eq_val)
      call WriteArray(f%fh, obj%Psi_eq_val)
      call WriteArray(f%fh, obj%a_E_val)
      call WriteArray(f%fh, obj%a_v_val)
      call WriteArray(f%fh, obj%E_eq_val)
      call WriteArray(f%fh, obj%v_eq_val)
      call WriteArray(f%fh, obj%theta_ps_val) ! computed from other )
      write (f%fh, *) obj%Name
      write (f%fh, *) obj%timestep
      write (f%fh, *) obj%dt
      call f%close()

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

! ################################################
   subroutine importWaterAbsorption(obj, Water, Tissue, a_Psi, a_P, theta_eq, theta_ps, &
                                    Psi_eq, a_E, a_v, E_eq, v_eq)
      class(WaterAbsorption_), intent(inout) :: obj
      type(FEMDomain_), target, optional, intent(inout) :: Water, Tissue

      type(MaterialProp_), target, optional, intent(in):: a_Psi
      type(MaterialProp_), target, optional, intent(in):: a_P
      type(MaterialProp_), target, optional, intent(in):: theta_eq
      type(MaterialProp_), target, optional, intent(in):: theta_ps
      type(MaterialProp_), target, optional, intent(in):: Psi_eq
      type(MaterialProp_), target, optional, intent(in):: a_E
      type(MaterialProp_), target, optional, intent(in):: a_v
      type(MaterialProp_), target, optional, intent(in):: E_eq
      type(MaterialProp_), target, optional, intent(in):: v_eq

      if (associated(obj%DiffusionEq%FEMDomain)) then
         nullify (obj%DiffusionEq%FEMDomain)
      end if

      if (associated(obj%FiniteDeform%FEMDomain)) then
         nullify (obj%FiniteDeform%FEMDomain)
      end if
      if (present(water)) then
         obj%DiffusionEq%FEMDomain => Water
         obj%water => water
         Water%SolverType = "DiffusionEq_"
      end if
      if (present(tissue)) then
         obj%FiniteDeform%FEMDomain => Tissue
         Tissue%SolverType = "FiniteDeform_"
         obj%Tissue => Tissue
      end if

      print *, "[ImportFile]>"

      if (present(a_Psi)) then
         if (associated(obj%a_Psi)) then
            nullify (obj%a_Psi)
         end if
         obj%a_Psi => a_Psi
      end if

      if (present(a_P)) then
         if (associated(obj%a_P)) then
            nullify (obj%a_P)
         end if
         obj%a_P => a_P
      end if

      if (present(theta_eq)) then
         if (associated(obj%theta_eq)) then
            nullify (obj%theta_eq)
         end if
         obj%theta_eq => theta_eq
      end if
      if (present(Psi_eq)) then
         if (associated(obj%Psi_eq)) then
            nullify (obj%Psi_eq)
         end if
         obj%Psi_eq => Psi_eq
      end if
      if (present(a_E)) then
         if (associated(obj%a_E)) then
            nullify (obj%a_E)
         end if
         obj%a_E => a_E
      end if
      if (present(a_v)) then
         if (associated(obj%a_v)) then
            nullify (obj%a_v)
         end if
         obj%a_v => a_v
      end if

      if (present(E_eq)) then
         if (associated(obj%E_eq)) then
            nullify (obj%E_eq)
         end if
         obj%E_eq => E_eq
      end if

      if (present(v_eq)) then
         if (associated(obj%v_eq)) then
            nullify (obj%v_eq)
         end if
         obj%v_eq => v_eq
      end if

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

!#####################################
   subroutine runWaterAbsorption(obj, timestep, dt, SolverType, onlyInit, Only1st, Display, nr_tol, ReducedIntegration, &
                                 infinitesimal, interval, Name, restart, ID)
      class(WaterAbsorption_), intent(inout) :: obj
      character(*), intent(in) :: SolverType
      character(*), optional, intent(in) :: Name
      integer(int32) :: i, n, interv, coun
      logical, optional, intent(in) :: onlyInit, Only1st, Display, &
                                       ReducedIntegration, infinitesimal, restart
      real(real64), optional, intent(in) :: nr_tol
      integer(int32), intent(in) :: timestep
      integer(int32), optional, intent(in) :: interval, ID
      real(real64), optional, intent(in) :: dt

      if (present(restart)) then
         if (restart .eqv. .true.) then
            do i = 1, timestep
               call obj%update(SolverType, Display=.true., step=i, nr_tol=nr_tol, restart=.true.)
               print *, "Timestep :: ", i, "/", timestep
            end do
            return
         end if
      end if
      if (present(Name)) then
         obj%water%FileName = Name
         obj%water%Name = Name
         obj%tissue%FileName = Name
         obj%tissue%Name = Name
      end if
      interv = input(default=0, option=interval)
      if (present(ReducedIntegration)) then
         obj%FiniteDeform%ReducedIntegration = ReducedIntegration
      end if

      obj%dt = input(default=1.0d0, option=dt)
      !obj%timestep=timestep

      obj%DiffusionEq%dt = obj%dt
      obj%FiniteDeform%dt = obj%dt

      print *, "[ImportFile]>>>[Initialize]>"

      call obj%init(SolverType, Display=Display, nr_tol=nr_tol, infinitesimal=infinitesimal)
      print *, "[ImportFile]>>>[Initialize]>>>[1st Step]>"
      if (present(onlyInit)) then
         return
      end if
      ! Repeat over time
      coun = 0
      do i = 2, timestep

         if (coun == interv) then
            coun = 0
            call obj%update(SolverType, Display=Display, step=i, nr_tol=nr_tol)
         else
            coun = coun + 1
            call obj%update(SolverType, Display=.false., step=i, nr_tol=nr_tol)
         end if
         print *, "Timestep :: ", i, "/", timestep

         !n=1
         !call showArray(obj%FiniteDeform%DeformVecGloTot,&
         !Name="test"//fstring(n)//".txt")  )

         if (present(only1st)) then
            return
         end if
      end do

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

!#####################################
   subroutine initWaterAbsorption(obj, SolverType, Display, nr_tol, infinitesimal)
      class(WaterAbsorption_), intent(inout) :: obj
      type(Term_)             :: term
      character(*), intent(in) :: SolverType
      logical, optional, intent(in) :: Display, infinitesimal
      real(real64), optional, intent(in) :: nr_tol
      integer(int32) :: n

      call term%init()
      obj%volume = obj%FiniteDeform%getVolume()

      !open(113,file=obj%tissue%name//"x_length.txt",status="replace")
      !open(114,file=obj%tissue%name//"y_length.txt",status="replace")
      !open(115,file=obj%tissue%name//"z_length.txt",status="replace")
      call obj%updatePermiability()
      ! ###### Diffusion Part ###################
      call obj%DiffusionEq%Setup()
      call obj%DiffusionEq%Solve(SolverType=SolverType)
      if (present(Display)) then
         if (Display .eqv. .true.) then
            call DisplayDiffusionEq(obj%DiffusionEq, &
                                    DisplayMode=term%Gmsh, OptionalStep=1)
         end if
      end if

      ! ###### Diffusion Part ###################
      print *, "[ImportFile]>>>[Initialize]>>>[1st Step]>> water >>"

      ! debug
      ! only diffusion

      ! ###### Finite deformation part #############################
      !call obj%updateStiffness()
      !call obj%updateStiffness()
      !call obj%FiniteDeform%import(YoungsModulus=obj%YoungsModulus)
      !call obj%FiniteDeform%import(PoissonsRatio=obj%PoissonsRatio)
      !call obj%FiniteDeform%import(PorePressure=obj%PorePressure)
      ! ###### Finite deformation part #############################
      if (present(infinitesimal)) then
         obj%FiniteDeform%infinitesimal = infinitesimal
      end if
      call obj%FiniteDeform%DivideBC()
      call obj%FiniteDeform%Solve(SolverType=SolverType, nr_tol=nr_tol)
      call DisplayReactionForce(obj%FiniteDeform, obj%ID)
      if (present(Display)) then
         if (Display .eqv. .true.) then
            call DisplayDeformStress(obj%FiniteDeform, &
                                     DisplayMode="gmsh", OptionalStep=1)
            n = obj%FiniteDeform%step
            !call showArray(obj%FiniteDeform%DeformVecGloTot,&
            !Name="test"//fstring(n)//".txt")  )

         end if
      end if
      call obj%UpdateVolume()

      ! update mesh of Diffusion analysis added at 20200827
      obj%DiffusionEq%FEMDomain%Mesh%NodCoord = obj%FiniteDeform%FEMDomain%Mesh%NodCoord
      ! ###### Finite deformation part #############################
      print *, "[ImportFile]>>>[Initialize]>>>[1st Step]>>>>["

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

!#####################################
   subroutine updateWaterAbsorption(obj, SolverType, Display, step, nr_tol, restart)
      class(WaterAbsorption_), intent(inout) :: obj
      character(*), intent(in) :: SolverType
      logical, optional, intent(in) :: Display, restart
      integer(int32), optional, intent(in) :: step
      type(Term_)             :: term
      real(real64), optional, intent(in) :: nr_tol
      integer(int32) :: i, n, m

      call term%init()
      obj%volume = obj%FiniteDeform%getVolume()

      obj%DiffusionEq%FEMDomain%name = "DiffusionEq"
      obj%FiniteDeform%FEMDomain%name = "FiniteDeform"
      obj%DiffusionEq%FEMDomain%filename = "DiffusionEq"
      obj%FiniteDeform%FEMDomain%filename = "FiniteDeform"
      ! ###### update DIffusion parameters ############
      call obj%updatePermiability()
      call obj%DiffusionEq%import(Permiability=obj%Permiability)
      !do i=1,size(obj%water%mesh%nodCoord,1)
      !    write(1234,*) obj%water%mesh%nodCoord(i,1),obj%DiffusionEq%dt*dble(step), obj%DiffusionEq%UnknownVec(i)
      !enddo
      !write(1234,*) "  "
      ! ###### Update Diffusion Field over timesteps ###################
      call obj%DiffusionEq%Update()
      call obj%DiffusionEq%Solve(SolverType=SolverType, restart=restart)
      obj%DiffusionEq%step = step
      if (present(Display)) then
         if (Display .eqv. .true.) then

            call DisplayDiffusionEq(obj%DiffusionEq, &
                                    DisplayMode=term%Gmsh, OptionalStep=step, name="water")
         end if
      end if
      ! ###### Update Diffusion Field over timesteps ###################

      ! only diffusion debug
      call obj%updateStiffness()
      ! ###### update DIffusion parameters ############
      call obj%FiniteDeform%import(YoungsModulus=obj%YoungsModulus)
      call obj%FiniteDeform%import(PoissonsRatio=obj%PoissonsRatio)
      call obj%FiniteDeform%import(PorePressure=obj%PorePressure)
      ! ###### Update Finite Deformation over timesteps ###################
      obj%FiniteDeform%step = step

      call obj%FiniteDeform%UpdateInitConfig()
      call obj%FiniteDeform%UpdateBC()
      call obj%FiniteDeform%Solve(SolverType=SolverType, nr_tol=nr_tol, restart=restart)

      ! update mesh of Diffusion analysis added at 20200827
      obj%DiffusionEq%FEMDomain%Mesh%NodCoord = obj%FiniteDeform%FEMDomain%Mesh%NodCoord

      call DisplayReactionForce(obj%FiniteDeform, obj%ID)
      print *, "Timestep ", step
      if (present(Display)) then
         if (Display .eqv. .true.) then
            call DisplayDeformStress(obj%FiniteDeform, &
                                     DisplayMode="gmsh", OptionalStep=step, name="tissue")

            call obj%export(displacement=.true.)
         end if
      end if

      ! ###### Update Finite Deformation over timesteps ###################

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

   subroutine updatePermiabilityWA(obj)
      class(WaterAbsorption_), intent(inout) :: obj
      integer(int32) :: i, j, n, m
      real(real64) :: k_0, a_hat, val(2), x
      n = size(obj%a_Psi_val)

      if (.not. allocated(obj%Permiability)) then
         allocate (obj%Permiability(n))
      end if

      if (.not. allocated(obj%WaterContent)) then
         allocate (obj%WaterContent(n))
         allocate (obj%WaterMass(n))
         obj%WaterContent(:) = 0.0d0
         obj%WaterMass(:) = 0.0d0
      end if

      if (allocated(obj%DiffusionEq%UnknownValue)) then
         if (size(obj%DiffusionEq%UnknownValue, 1) == n) then
            m = size(obj%DiffusionEq%UnknownValue, 2)
            do i = 1, n
               obj%WaterContent(i) = 0.0d0
               do j = 1, m
                  ! get avarage value of wate content theta
                  obj%WaterContent(i) = obj%WaterContent(i) + obj%DiffusionEq%UnknownValue(i, j)/dble(m)
               end do
            end do
         end if
      end if

      ! determine permiability for each element
      do i = 1, n
         k_0 = obj%Water%MaterialProp%matPara(obj%Water%Mesh%ElemMat(i), 1)
         val(1) = obj%a_Psi_val(i)
         val(2) = obj%a_Psi_val(i) - abs(obj%a_P_val(i))
         ! Hebiside step function
         !if(obj%WaterContent(i) <=obj%theta_ps_val(i) )then
         !    a_hat = val(1)
         !else
         !    a_hat=val(2)
         !endif
         ! use normalization
         x = 2.0d0*(obj%WaterContent(i) - 0.50d0)*100.0d0
         a_hat = obj%a_Psi_val(i) - 1.0d0/(1.0d0 + exp(-x))*abs(obj%a_P_val(i))
         !write(220,*) obj%WaterContent(i), a_hat*k_0
         obj%Permiability(i) = k_0*a_hat
      end do

      !write(100,*) obj%Permiability(:)
   end subroutine updatePermiabilityWA
!#####################################

!#####################################
   subroutine updateStiffnessWA(obj)
      class(WaterAbsorption_), intent(inout) :: obj
      integer(int32) :: i, j, n
      real(real64) :: a_E, E, E_eq, theta, theta_eq, a_P, theta_ps
      real(real64) :: a_v, v, v_eq, p, val(2)
      n = size(obj%water%mesh%elemNod, 1)
      if (.not. allocated(obj%YoungsModulus)) then
         allocate (obj%YoungsModulus(n))
      end if

      if (.not. allocated(obj%PoissonsRatio)) then
         allocate (obj%PoissonsRatio(n))
      end if

      if (.not. allocated(obj%PorePressure)) then
         allocate (obj%PorePressure(n))
         obj%PorePressure = 0.0d0
      end if

      ! update these parameters considering volumetric water content
      do i = 1, n
         ! import data for each element
         a_E = obj%a_E_val(i)
         a_v = obj%a_v_val(i)
         E_eq = obj%E_eq_val(i)
         v_eq = obj%v_eq_val(i)
         theta = obj%WaterContent(i)
         theta_eq = obj%theta_eq_val(i)
         theta_ps = obj%theta_ps_val(i)
         a_P = obj%a_P_val(i)

         ! compute parameters
         if (a_E > 0) then
            a_E = -abs(a_E)
         end if

         if (theta > 1.0d0) then
            theta = 1.0d0
         end if
         E = a_E*(theta - theta_eq) + E_eq
         v = a_v*(theta - theta_eq) + v_eq
         if (theta > theta_ps) then
            p = a_P*(theta - theta_ps)
         else
            p = 0.0d0
         end if

         ! export parameters
         obj%YoungsModulus(i) = E
         obj%PoissonsRatio(i) = v
         obj%PorePressure(i) = p !- obj%PorePressure(i)
      end do
      print *, "maxval(obj%a_P_val(:))", maxval(obj%a_P_val(:)), minval(obj%a_P_val(:))
      print *, "maxval(obj%WaterContent(:))", maxval(obj%WaterContent(:)), minval(obj%WaterContent(:))
      print *, "maxval(obj%theta_ps_val(:))", maxval(obj%theta_ps_val(:)), minval(obj%theta_ps_val(:))
      print *, "maxval(obj%PorePressure(:))", maxval(obj%PorePressure(:)), minval(obj%PorePressure(:))
      call obj%FiniteDeform%import(obj%YoungsModulus, obj%PoissonsRatio, obj%PorePressure)

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

!#####################################
   subroutine exportWaterAbsorption(obj, OptionalFileFormat, OptionalProjectName, FileHandle, &
                                    SolverType, MeshDimension, FileName, Name, regacy, with, displacement, restart, path)
      class(WaterAbsorption_), intent(inout) :: obj
      class(FEMDomain_), optional, intent(inout)::with
      character(*), optional, intent(in)::OptionalFileFormat, path
      character(*), optional, intent(in)::OptionalProjectName, SolverType, FileName
      character*4::FileFormat
      character(*), optional, intent(in) :: Name
      logical, optional, intent(in) :: regacy, displacement, restart
      character*200::ProjectName
      character*200 ::iFileName, fname
      integer(int32), allocatable::IntMat(:, :)
      real(real64), allocatable::RealMat(:, :)
      integer(int32), optional, intent(in)::FileHandle, MeshDimension
      integer(int32) :: fh, i, j, k, NumOfDomain, n, m, DimNum, GpNum, nn, fh_
      character*70 Msg
      type(IO_) :: f

      if (present(restart)) then
         if (.not. present(path)) then
            print *, " exportWaterAbsorption ERROR >> .not. present(path)"
            stop
         end if

         if (present(name)) then
            call execute_command_line("mkdir -p "//path//"/"//name)
            call f%open("./", path//"/"//name, "/WaterAbsorption.prop")
            fname = path//"/"//name
            write (f%fh, *) obj%WaterAbsorbingPower(:)
            write (f%fh, *) obj%WaterPotential(:)
            write (f%fh, *) obj%TurgorPressure(:)
            write (f%fh, *) obj%WaterContent(:)
            write (f%fh, *) obj%Conductivity(:)
            write (f%fh, *) obj%YoungsModulus(:)
            write (f%fh, *) obj%PoissonsRatio(:)
            write (f%fh, *) obj%Permiability(:)
            write (f%fh, *) obj%PorePressure(:)
            write (f%fh, *) obj%a_Psi_val(:)
            write (f%fh, *) obj%a_P_val(:)
            write (f%fh, *) obj%theta_eq_val(:)
            write (f%fh, *) obj%Psi_eq_val(:)
            write (f%fh, *) obj%a_E_val(:)
            write (f%fh, *) obj%a_v_val(:)
            write (f%fh, *) obj%E_eq_val(:)
            write (f%fh, *) obj%v_eq_val(:)
            write (f%fh, *) obj%theta_ps_val(:) ! computed from other variables
            write (f%fh, *) obj%Name
            write (f%fh, *) obj%timestep
            write (f%fh, *) obj%dt
            call f%close()

            call execute_command_line("mkdir -p ./"//fname//"/FEMDomain")
            call f%open("./", fname, "/FEMDomain/Water.prop")
            call obj%Water%export(path="./"//fname//"/FEMDomain", restart=restart)
            call f%close()

            call f%open("./", fname, "/FEMDomain/Tissue.prop")
            call obj%Tissue%export(path="./"//fname//"/FEMDomain", restart=restart)
            call f%close()

            call execute_command_line("mkdir -p ./"//fname//"/DiffusionEq")
            call f%open("./", fname, "/DiffusionEq/DiffusionEq.prop")
            call obj%DiffusionEq%export(path="./"//fname//"/DiffusionEq", restart=restart)
            call f%close()

            call execute_command_line("mkdir -p ./"//fname//"/FiniteDeform")
            call f%open("./", fname, "/FiniteDeform/FiniteDeform.prop")
            call obj%FiniteDeform%export(path="./"//fname//"/FiniteDeform", restart=restart)
            call f%close()

            call execute_command_line("mkdir -p ./"//fname//"/a_Psi")
            call obj%a_Psi%export(path=fname//"/a_Psi", restart=restart)

            call execute_command_line("mkdir -p ./"//fname//"/a_P")
            call obj%a_P%export(path=fname//"/a_P", restart=restart)

            call execute_command_line("mkdir -p ./"//fname//"/theta_eq")
            call obj%theta_eq%export(path=fname//"/theta_eq", restart=restart)

            call execute_command_line("mkdir -p ./"//fname//"/Psi_eq")
            call obj%Psi_eq%export(path=fname//"/Psi_eq", restart=restart)

            call execute_command_line("mkdir -p ./"//fname//"/a_E")
            call obj%a_E%export(path=fname//"/a_E", restart=restart)

            call execute_command_line("mkdir -p ./"//fname//"/a_v")
            call obj%a_v%export(path=fname//"/a_v", restart=restart)

            call execute_command_line("mkdir -p ./"//fname//"/E_eq")
            call obj%E_eq%export(path=fname//"/E_eq", restart=restart)

            call execute_command_line("mkdir -p ./"//fname//"/v_eq")
            call obj%v_eq%export(path=fname//"/v_eq", restart=restart)

         else
            call execute_command_line("mkdir -p "//path//"/WaterAbsorption")
            call f%open("./", path//"/WaterAbsorption", "/WaterAbsorption.prop")
            fname = path//"/WaterAbsorption"
            write (f%fh, *) obj%WaterAbsorbingPower(:)
            write (f%fh, *) obj%WaterPotential(:)
            write (f%fh, *) obj%TurgorPressure(:)
            write (f%fh, *) obj%WaterContent(:)
            write (f%fh, *) obj%Conductivity(:)
            write (f%fh, *) obj%YoungsModulus(:)
            write (f%fh, *) obj%PoissonsRatio(:)
            write (f%fh, *) obj%Permiability(:)
            write (f%fh, *) obj%PorePressure(:)
            write (f%fh, *) obj%a_Psi_val(:)
            write (f%fh, *) obj%a_P_val(:)
            write (f%fh, *) obj%theta_eq_val(:)
            write (f%fh, *) obj%Psi_eq_val(:)
            write (f%fh, *) obj%a_E_val(:)
            write (f%fh, *) obj%a_v_val(:)
            write (f%fh, *) obj%E_eq_val(:)
            write (f%fh, *) obj%v_eq_val(:)
            write (f%fh, *) obj%theta_ps_val(:) ! computed from other variables
            write (f%fh, *) obj%Name
            write (f%fh, *) obj%timestep
            write (f%fh, *) obj%dt
            call f%close()

            call execute_command_line("mkdir -p ./"//fname//"/FEMDomain")
            call f%open("./", fname, "/FEMDomain/Water.prop")
            call obj%Water%export(path="./"//fname//"/FEMDomain", restart=restart)
            call f%close()

            call f%open("./", fname, "/FEMDomain/Tissue.prop")
            call obj%Tissue%export(path="./"//fname//"/FEMDomain", restart=restart)
            call f%close()

            call execute_command_line("mkdir -p ./"//fname//"/DiffusionEq")
            call f%open("./", fname, "/DiffusionEq/DiffusionEq.prop")
            call obj%DiffusionEq%export(path="./"//fname//"/DiffusionEq", restart=restart)
            call f%close()

            call execute_command_line("mkdir -p ./"//fname//"/FiniteDeform")
            call f%open("./", fname, "/FiniteDeform/FiniteDeform.prop")
            call obj%FiniteDeform%export(path="./"//fname//"/FiniteDeform", restart=restart)
            call f%close()

            call execute_command_line("mkdir -p ./"//fname//"/a_Psi")
            call obj%a_Psi%export(path=fname//"/a_Psi", restart=restart)

            call execute_command_line("mkdir -p ./"//fname//"/a_P")
            call obj%a_P%export(path=fname//"/a_P", restart=restart)

            call execute_command_line("mkdir -p ./"//fname//"/theta_eq")
            call obj%theta_eq%export(path=fname//"/theta_eq", restart=restart)

            call execute_command_line("mkdir -p ./"//fname//"/Psi_eq")
            call obj%Psi_eq%export(path=fname//"/Psi_eq", restart=restart)

            call execute_command_line("mkdir -p ./"//fname//"/a_E")
            call obj%a_E%export(path=fname//"/a_E", restart=restart)

            call execute_command_line("mkdir -p ./"//fname//"/a_v")
            call obj%a_v%export(path=fname//"/a_v", restart=restart)

            call execute_command_line("mkdir -p ./"//fname//"/E_eq")
            call obj%E_eq%export(path=fname//"/E_eq", restart=restart)

            call execute_command_line("mkdir -p ./"//fname//"/v_eq")
            call obj%v_eq%export(path=fname//"/v_eq", restart=restart)

         end if

         return
      end if

      if (present(displacement)) then
         if (displacement .eqv. .true.) then
            if (present(Name)) then
               open (12, file="x"//name)
               open (13, file="y"//name)
               open (14, file="z"//name)
            else
               open (12, file="x_Disp.txt")
               open (13, file="y_Disp.txt")
               open (14, file="z_Disp.txt")
            end if
            n = size(obj%water%mesh%nodCoord, 1)
            do i = 1, n
               write (12, *) obj%water%mesh%nodCoord(i, 1), obj%tissue%mesh%nodCoord(i, 1) - obj%water%mesh%nodCoord(i, 1)
               write (13, *) obj%water%mesh%nodCoord(i, 2), obj%tissue%mesh%nodCoord(i, 2) - obj%water%mesh%nodCoord(i, 2)
               write (14, *) obj%water%mesh%nodCoord(i, 3), obj%tissue%mesh%nodCoord(i, 3) - obj%water%mesh%nodCoord(i, 3)
            end do
            close (12)
            close (13)
            close (14)
         end if
      end if
      ! bug exsits
      !call obj%Tissue%export(SolverType="FiniteDeform",Name=Name)
      !call obj%Water%export(SolverType="DiffusionEq",Name=Name)

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

!#####################################
   subroutine displayWaterAbsorption(obj)
      class(WaterAbsorption_), intent(inout) :: obj

      ! implement how to display the results.

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

!#####################################
   subroutine bakeWaterAbsorption(obj)
      class(WaterAbsorption_), intent(inout) :: obj
      integer(int32) :: i, j, k, l, n, m, NumOfMaterial, layer, in_num, NumOfLayer
      real(real64), allocatable :: matPara(:, :), info(:, :)
      integer(int32), allocatable :: key(:)
      type(Rectangle_) :: rect, mrect
      logical :: in_case
      real(real64) :: matparaval, coord(3), x_max(3), x_min(3)

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

      n = size(obj%water%mesh%ElemNod, 1)

      ! initialize all parameters
      if (allocated(obj%WaterAbsorbingPower)) then
         deallocate (obj%WaterAbsorbingPower)
      end if
      allocate (obj%WaterAbsorbingPower(n))
      if (allocated(obj%WaterPotential)) then
         deallocate (obj%WaterPotential)
      end if
      allocate (obj%WaterPotential(n))
      if (allocated(obj%TurgorPressure)) then
         deallocate (obj%TurgorPressure)
      end if
      allocate (obj%TurgorPressure(n))
      if (allocated(obj%WaterContent)) then
         deallocate (obj%WaterContent)
      end if
      allocate (obj%WaterContent(n))
      if (allocated(obj%Conductivity)) then
         deallocate (obj%Conductivity)
      end if
      allocate (obj%Conductivity(n))
      if (allocated(obj%YoungsModulus)) then
         deallocate (obj%YoungsModulus)
      end if
      allocate (obj%YoungsModulus(n))
      if (allocated(obj%PoissonsRatio)) then
         deallocate (obj%PoissonsRatio)
      end if
      allocate (obj%PoissonsRatio(n))

      ! initialize parameters
      obj%WaterAbsorbingPower(:) = 0.0d0 ! kPa
      obj%WaterPotential(:) = 0.0d0 ! water
      obj%TurgorPressure(:) = 0.0d0 ! changes with theta
      obj%WaterContent(:) = 1.0d0 ! theta
      obj%Conductivity(:) = 1.0d0 ! changes with theta
      obj%YoungsModulus(:) = 1.0d0 ! E
      obj%PoissonsRatio(:) = 0.30d0 ! v
      if (.not. associated(obj%a_Psi)) then
         print *, "a_Psi is not associated."
         stop
      end if
      if (.not. associated(obj%a_P)) then
         print *, "a_P is not associated."
         stop
      end if
      if (.not. associated(obj%theta_eq)) then
         print *, "theta_eq is not associated."
         stop
      end if
      if (.not. associated(obj%Psi_eq)) then
         print *, "Psi_eq is not associated."
         stop
      end if
      if (.not. associated(obj%a_E)) then
         print *, "a_E is not associated."
         stop
      end if
      if (.not. associated(obj%a_v)) then
         print *, "a_v is not associated."
         stop
      end if
      if (.not. associated(obj%E_eq)) then
         print *, "E_eq is not associated."
         stop
      end if
      if (.not. associated(obj%v_eq)) then
         print *, "v_eq is not associated."
         stop
      end if

      ! a_Psi
      call obj%a_Psi%getValues(mesh=obj%Tissue%Mesh, Values=obj%a_Psi_val)
      call obj%a_P%getValues(mesh=obj%Tissue%Mesh, Values=obj%a_P_val)
      call obj%theta_eq%getValues(mesh=obj%Tissue%Mesh, Values=obj%theta_eq_val)
      call obj%Psi_eq%getValues(mesh=obj%Tissue%Mesh, Values=obj%Psi_eq_val)
      call obj%a_E%getValues(mesh=obj%Tissue%Mesh, Values=obj%a_E_val)
      call obj%a_v%getValues(mesh=obj%Tissue%Mesh, Values=obj%a_v_val)
      call obj%E_eq%getValues(mesh=obj%Tissue%Mesh, Values=obj%E_eq_val)
      call obj%v_eq%getValues(mesh=obj%Tissue%Mesh, Values=obj%v_eq_val)

      ! a_Psi should be negative values
      n = size(obj%theta_eq_val)
      do i = 1, n
         obj%a_Psi_val(i) = -abs(obj%a_Psi_val(i))
      end do

      ! a_Psi should be negative values
      n = size(obj%theta_eq_val)
      do i = 1, n
         obj%a_E_val(i) = -abs(obj%a_E_val(i))
      end do

      ! a_Psi should be negative values
      n = size(obj%theta_eq_val)
      do i = 1, n
         obj%a_v_val(i) = abs(obj%a_v_val(i))
      end do

      n = size(obj%theta_eq_val)
      if (.not. allocated(obj%theta_ps_val)) then
         allocate (obj%theta_ps_val(n))
      end if

      do i = 1, n
         obj%theta_ps_val(i) = (-obj%Psi_eq_val(i) + obj%a_P_val(i)*obj%theta_eq_val(i))/obj%a_P_val(i)
      end do

      !call showarray(mat=obj%WaterAbsorbingPower, Name="test.txt")

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

!#####################################
   subroutine gnuplotWaterAbsorption(obj, mode, ElemID)
      class(WaterAbsorption_), intent(in) :: obj
      character(*), intent(in) :: mode
      integer(int32), optional, intent(in) :: ElemID
      integer(int32) :: i, j, n, m
      real(real64) :: theta, psi(2)

      n = input(default=1, option=ElemID)
      if (mode == "OWCC" .or. mode == "all") then
         open (20, file="OWCC.txt")
         do i = 0, 100
            theta = dble(i)*obj%theta_eq_val(n)/100.0d0
            write (20, *) theta, (theta - obj%theta_eq_val(n))*obj%a_Psi_val(n) + obj%Psi_eq_val(n)
         end do
         write (20, *) " "
         do i = 0, 100
            theta = dble(i)*obj%theta_eq_val(n)/100.0d0
            psi(1) = 0.0d0
            psi(2) = obj%a_P_val(n)*(theta - obj%theta_ps_val(n))
            write (20, *) theta, maxval(psi)
         end do
         write (20, *) " "
         close (20)
      end if

      if (mode == "YoungsModulus" .or. mode == "all") then
         open (20, file="YoungsModulus.txt")
         do i = 0, 100
            theta = dble(i)*obj%theta_eq_val(n)/100.0d0
            write (20, *) theta, (theta - obj%theta_eq_val(n))*obj%a_E_val(n) + obj%E_eq_val(n)
         end do
         close (20)
      end if

      if (mode == "PoissonsRatio" .or. mode == "all") then
         open (20, file="PoissonsRatio.txt")
         do i = 0, 100
            theta = dble(i)*obj%theta_eq_val(n)/100.0d0
            write (20, *) theta, (theta - obj%theta_eq_val(n))*obj%a_v_val(n) + obj%v_eq_val(n)
         end do
         close (20)
      end if
   end subroutine gnuplotWaterAbsorption
!#####################################

!#####################################
   subroutine UpdateVolumeWA(obj)
      class(WaterAbsorption_), intent(inout) :: obj
      real(real64), allocatable :: volume(:)
      real(real64) :: theta_old, theta_new, mw, vol_old, vol_new
      integer(int32) :: i, j, n, numelem

      numelem = size(obj%FiniteDeform%FEMDomain%Mesh%ElemNod, 1)
      volume = obj%FiniteDeform%getVolume()
      if (.not. allocated(obj%WaterMass)) then
         allocate (obj%WaterMass(numelem))
      end if

      do i = 1, numelem
         theta_old = obj%WaterContent(i)
         vol_old = obj%volume(i)
         mw = theta_old*vol_old
         theta_new = mw/volume(i)
         obj%WaterMass(i) = mw
         obj%WaterContent(i) = theta_new
      end do
      obj%volume = obj%FiniteDeform%getVolume()

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

end module
