
module FiniteDeformationClass
   use, intrinsic :: iso_fortran_env
   use MathClass
   use LinearSolverClass
   use FEMDomainClass
   use PostProcessingClass
   use ConstitutiveModelClass
   implicit none

   type:: FiniteDeform_
      type(FEMDomain_), pointer ::FEMDomain
      real(real64), allocatable ::DeformStress(:, :, :)
      real(real64), allocatable ::DeformStrain(:, :, :)
      real(real64), allocatable ::DeformStressInit(:, :, :)
      real(real64), allocatable ::DeformStressinc(:, :, :)
      real(real64), allocatable ::DeformStressMat(:, :, :)
      real(real64), allocatable ::DeformStressRHS(:, :)
      real(real64), allocatable ::DeformVecEBETot(:, :)
      real(real64), allocatable ::DeformVecEBEInc(:, :)
      real(real64), allocatable ::DeformVecGloTot(:)
      real(real64), allocatable ::DeformVecGloInc(:)
      real(real64), allocatable ::TractionVecGlo(:)
      real(real64), allocatable ::ResidualVecGlo(:)
      real(real64), allocatable ::InternalVecGlo(:)
      real(real64), allocatable ::VolInitCurrEBE(:, :)
      real(real64), allocatable ::YoungsModulus(:) ! directly give parameter #1
      real(real64), allocatable ::PoissonsRatio(:) ! directly give parameter #2
      real(real64), allocatable :: PorePressure(:) ! directly give parameter #2
      real(real64)             ::dt, error, reactionforce
      real(real64)             ::nr_tol = 1.0e-8
      logical :: ReducedIntegration = .false.
      logical :: infinitesimal = .false.

      integer(int32) :: itr
      integer(int32) :: Step = 0
   contains
      procedure :: Solve => SolveFiniteDeformNewton
      procedure :: UpdateSolution => SolveFiniteDeform
      procedure :: DivideBC => DevideBCIntoTimestep
      procedure :: UpdateBC => UpdateBCInTimestep
      procedure :: UpdateInitConfig => UpdateInitConfig
      procedure :: Setup => SetupFiniteDeform
      procedure :: Update => UpdateFiniteDeform
      procedure :: Display => DisplayDeformStress
      procedure :: getDBCVector => getDBCVectorDeform
      procedure :: getDispVector => getDispVectorDeform
      procedure :: getVolume => getVolumeDeform
      procedure :: check => checkFiniteDeform
      procedure :: import => importFiniteDeform
      procedure :: export => exportFiniteDeform
      procedure :: remove => removeFiniteDeform

      procedure :: save => saveFiniteDeform
      procedure :: open => openFiniteDeform
      procedure :: link => linkFiniteDeform

   end type

contains

   subroutine removeFiniteDeform(obj)
      class(FiniteDeform_), intent(inout) :: obj

      if (associated(obj%FEMDomain)) then
         nullify (obj%FEMDomain)
      end if
      if (allocated(obj%DeformStress)) deallocate (obj%DeformStress)
      if (allocated(obj%DeformStrain)) deallocate (obj%DeformStrain)
      if (allocated(obj%DeformStressInit)) deallocate (obj%DeformStressInit)
      if (allocated(obj%DeformStressinc)) deallocate (obj%DeformStressinc)
      if (allocated(obj%DeformStressMat)) deallocate (obj%DeformStressMat)
      if (allocated(obj%DeformStressRHS)) deallocate (obj%DeformStressRHS)
      if (allocated(obj%DeformVecEBETot)) deallocate (obj%DeformVecEBETot)
      if (allocated(obj%DeformVecEBEInc)) deallocate (obj%DeformVecEBEInc)
      if (allocated(obj%DeformVecGloTot)) deallocate (obj%DeformVecGloTot)
      if (allocated(obj%DeformVecGloInc)) deallocate (obj%DeformVecGloInc)
      if (allocated(obj%TractionVecGlo)) deallocate (obj%TractionVecGlo)
      if (allocated(obj%ResidualVecGlo)) deallocate (obj%ResidualVecGlo)
      if (allocated(obj%InternalVecGlo)) deallocate (obj%InternalVecGlo)
      if (allocated(obj%VolInitCurrEBE)) deallocate (obj%VolInitCurrEBE)
      if (allocated(obj%YoungsModulus)) deallocate (obj%YoungsModulus)
      if (allocated(obj%PoissonsRatio)) deallocate (obj%PoissonsRatio)
      if (allocated(obj%PorePressure)) deallocate (obj%PorePressure)
      obj%dt = 1.0d0
      obj%error = 0.0d0
      obj%reactionforce = 0.0d0
      obj%nr_tol = 1.0e-8
      obj%ReducedIntegration = .false.
      obj%infinitesimal = .false.

      obj%itr = 0
      obj%Step = 0
   end subroutine

! #######################################################################
   subroutine linkFiniteDeform(obj, FEMDomain)
      class(FiniteDeform_), intent(inout) :: obj
      type(FEMDomain_), target, intent(in) :: FEMDomain

      if (associated(obj%FEMDomain)) then
         nullify (obj%FEMDomain)
         obj%FEMDomain => FEMDOmain
      end if

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

! #######################################################################
   subroutine openFiniteDeform(obj, path, name)
      class(FiniteDeform_), intent(inout) :: obj
      character(*), intent(in) :: path
      character(*), optional, intent(in) :: name
      character(200) :: pathi
      type(IO_) :: f
      integer(int32) :: n

      if (present(name)) then
         pathi = path
         !if( index(path, "/", back=.true.) == len(path) )then
         !        n=index(path, "/", back=.true.)
         !        pathi(n:n)= " "
         !endif

         call execute_command_line("mkdir -p "//pathi)
         call execute_command_line("mkdir -p "//pathi//"/"//name)
         call f%open(pathi//"/"//name, "/"//"FiniteDeform", ".prop")
      else
         pathi = path
         !if( index(path, "/", back=.true.) == len(path) )then
         !        n=index(path, "/", back=.true.)
         !        pathi(n:n)= " "
         !endif
         call execute_command_line("mkdir -p "//pathi)
         call execute_command_line("mkdir -p "//pathi//"/FiniteDeform")
         call f%open(pathi//"/FiniteDeform", "/FiniteDeform", ".prop")
      end if

      ! write smt at here!

      call openArray(f%fh, obj%DeformStress)
      call openArray(f%fh, obj%DeformStrain)
      call openArray(f%fh, obj%DeformStressInit)
      call openArray(f%fh, obj%DeformStressinc)
      call openArray(f%fh, obj%DeformStressMat)
      call openArray(f%fh, obj%DeformStressRHS)
      call openArray(f%fh, obj%DeformVecEBETot)
      call openArray(f%fh, obj%DeformVecEBEInc)
      call openArray(f%fh, obj%DeformVecGloTot)
      call openArray(f%fh, obj%DeformVecGloInc)
      call openArray(f%fh, obj%TractionVecGlo)
      call openArray(f%fh, obj%ResidualVecGlo)
      call openArray(f%fh, obj%InternalVecGlo)
      call openArray(f%fh, obj%VolInitCurrEBE)
      call openArray(f%fh, obj%YoungsModulus) ! directly give parameter #1
      call openArray(f%fh, obj%PoissonsRatio) ! directly give parameter #2
      call openArray(f%fh, obj%PorePressure) ! directly give parameter #2
      read (f%fh, *) obj%dt, obj%error, obj%reactionforce
      read (f%fh, *) obj%nr_tol
      read (f%fh, *) obj%ReducedIntegration
      read (f%fh, *) obj%infinitesimal

      call f%close()

   end subroutine

! #######################################################################
   subroutine saveFiniteDeform(obj, path, name)
      class(FiniteDeform_), intent(inout) :: obj
      character(*), intent(in) :: path
      character(*), optional, intent(in) :: name
      character(200) :: pathi
      type(IO_) :: f
      integer(int32) :: n

      if (present(name)) then
         pathi = path
         !if( index(path, "/", back=.true.) == len(path) )then
         !        n=index(path, "/", back=.true.)
         !        pathi(n:n)= " "
         !endif

         call execute_command_line("mkdir -p "//pathi)
         call execute_command_line("mkdir -p "//pathi//"/"//name)
         call f%open(pathi//"/"//name, "/"//"FiniteDeform", ".prop")
      else
         pathi = path
         !if( index(path, "/", back=.true.) == len(path) )then
         !        n=index(path, "/", back=.true.)
         !        pathi(n:n)= " "
         !endif
         call execute_command_line("mkdir -p "//pathi)
         call execute_command_line("mkdir -p "//pathi//"/FiniteDeform")
         call f%open(pathi//"/FiniteDeform", "/FiniteDeform", ".prop")
      end if

      ! write smt at here!

      call writeArray(f%fh, obj%DeformStress)
      call writeArray(f%fh, obj%DeformStrain)
      call writeArray(f%fh, obj%DeformStressInit)
      call writeArray(f%fh, obj%DeformStressinc)
      call writeArray(f%fh, obj%DeformStressMat)
      call writeArray(f%fh, obj%DeformStressRHS)
      call writeArray(f%fh, obj%DeformVecEBETot)
      call writeArray(f%fh, obj%DeformVecEBEInc)
      call writeArray(f%fh, obj%DeformVecGloTot)
      call writeArray(f%fh, obj%DeformVecGloInc)
      call writeArray(f%fh, obj%TractionVecGlo)
      call writeArray(f%fh, obj%ResidualVecGlo)
      call writeArray(f%fh, obj%InternalVecGlo)
      call writeArray(f%fh, obj%VolInitCurrEBE)
      call writeArray(f%fh, obj%YoungsModulus) ! directly give parameter #1
      call writeArray(f%fh, obj%PoissonsRatio) ! directly give parameter #2
      call writeArray(f%fh, obj%PorePressure) ! directly give parameter #2
      write (f%fh, *) obj%dt, obj%error, obj%reactionforce
      write (f%fh, *) obj%nr_tol
      write (f%fh, *) obj%ReducedIntegration
      write (f%fh, *) obj%infinitesimal

      call f%close()

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

! ###############################################################
   subroutine importFiniteDeform(obj, YoungsModulus, PoissonsRatio, PorePressure)
      class(FiniteDeform_), intent(inout) :: obj
      real(real64), optional, intent(in) :: YoungsModulus(:), PoissonsRatio(:), &
                                            PorePressure(:)
      integer(int32) :: i, j, n

      if (present(YoungsModulus)) then
         if (allocated(obj%YoungsModulus)) then
            deallocate (obj%YoungsModulus)
         end if
         n = size(YoungsModulus)
         allocate (obj%YoungsModulus(n))
         obj%YoungsModulus(:) = YoungsModulus(:)
      end if

      if (present(PoissonsRatio)) then
         if (allocated(obj%PoissonsRatio)) then
            deallocate (obj%PoissonsRatio)
         end if
         n = size(PoissonsRatio)
         allocate (obj%PoissonsRatio(n))
         obj%PoissonsRatio(:) = PoissonsRatio(:)
      end if

      if (present(PorePressure)) then
         if (allocated(obj%PorePressure)) then
            deallocate (obj%PorePressure)
         end if
         n = size(PorePressure)
         allocate (obj%PorePressure(n))
         obj%PorePressure(:) = PorePressure(:)
      end if

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

! #######################################################
   subroutine checkFiniteDeform(obj)
      class(FiniteDeform_), intent(in) :: obj
      integer(int32) :: i, j, error_counter
      ! check conditions and return alartes.
      error_counter = 0
      print *, "checkFiniteDeform :: Analyzing imported data and checking consistency..."
      if (.not. associated(obj%FEMDomain)) then
         print *, "Alert! checkFiniteDeform >> FEMDomain is not improted."
         stop
      end if
      if (size(obj%FEMDomain%Boundary%DBoundNodID, 2) /= size(obj%FEMDomain%Mesh%NodCoord, 2)) then
         print *, "size(obj%FEMDomain%Boundary%DBoundNodID,2)/=size(obj%FEMDomain%Mesh%NodCoord,2)"
         stop
      end if

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

!######################## Solve deformation by Netwon's method ########################
   subroutine SolveFiniteDeformNewton(obj, OptionItr, Solvertype, nr_tol, infinitesimal, restart)
      class(FiniteDeform_), intent(inout)::obj
      integer(int32), optional, intent(in)::OptionItr
      character(*), optional, intent(in)::Solvertype
      real(real64), optional, intent(in)::nr_tol
      character*70 ::solver, defaultsolver
      logical, optional, intent(in) :: infinitesimal, restart

      real(real64), allocatable::Amat(:, :), bvec(:), xvec(:)
      real(real64)::val, er, residual, tolerance
      integer(int32) ::i, j, n, m, k, l, dim1, dim2, nodeid1, nodeid2, localid, itrmax, SetBC, itr_tol, itr
      integer(int32) :: dim_num, node_num, elem_num, node_num_elmtl
      character*70 :: gmsh, GaussJordan
      logical :: skip = .false.
      gmsh = "Gmsh"
      GaussJordan = "GaussJordan"
      if (present(infinitesimal)) then
         obj%infinitesimal = infinitesimal
      end if

      if (present(OptionItr)) then
         itr_tol = OptionItr
      else
         itr_tol = obj%FEMDomain%ControlPara%ItrTol
      end if
      if (present(nr_tol)) then
         obj%nr_tol = nr_tol
      end if
      if (present(restart)) then
         if (restart .eqv. .true.) then
            skip = .true.
         end if
      end if
      do
         itr = itr + 1
         obj%Itr = itr
         if (itr == 1 .and. skip .eqv. .false.) then

            call SetupFiniteDeform(obj)
            if (obj%Step == 1) then
               !call DisplayDeformStress(obj,DisplayMode=gmsh,OptionalStep=0)
            end if

            call SolveFiniteDeform(obj, SolverType=Solvertype)
            !call DisplayDeformStress(obj,DisplayMode=gmsh,OptionalStep=itr)
            if (obj%infinitesimal .eqv. .true.) then
               exit
            end if

         else

            call UpdateFiniteDeform(obj)

            !call DisplayDeformStress(obj,DisplayMode=gmsh,OptionalStep=0)
            call SolveFiniteDeform(obj, OptionItr=itr, SolverType=Solvertype)
            !call DisplayDeformStress(obj,DisplayMode=gmsh,OptionalStep=itr)

            !print *, "Residual r*r = ",dot_product(obj%ResidualVecGlo,obj%ResidualVecGlo)
            print *, "Step, Itr, ERROR :: ", obj%step, obj%itr, obj%error
            if (obj%error < obj%nr_tol) then
               print *, "itr=", itr, "Netwton's Method is converged !"
               exit
            end if
         end if

         if (itr == itr_tol) then
            print *, "SolveFiniteDeformNewton >> Newton's method did not converge"
            stop "debug"
            exit
         end if

      end do
      call UpdateStressMeasure(obj)
      call UpdateStrainMeasure(obj)

   end subroutine
!######################## Solve deformation by Netwon's method ########################

!########################## Initialize Boundary Conditions of Finite Deformation  ###################################
   subroutine DevideBCIntoTimestep(obj)
      class(FiniteDeform_), intent(inout)::obj

      integer(int32) ::n, m, timestep

      !debug Display Dirichlet Boundary Condition
      !call obj%FEMDomain%GmshPlotMesh(onlyDirichletBC=.true.)
      timestep = obj%FEMDomain%ControlPara%Timestep

      if (obj%FEMDomain%ControlPara%SimMode == 1) then

         n = size(obj%FEMDomain%Boundary%DBoundVal, 1)
         m = size(obj%FEMDomain%Boundary%DBoundVal, 2)

         if (.not. allocated(obj%FEMDomain%Boundary%DBoundValInc)) then
            allocate (obj%FEMDomain%Boundary%DBoundValInc(n, m))
         else
            if (size(obj%FEMDomain%Boundary%DBoundValInc, 1) /= n .or. &
                size(obj%FEMDomain%Boundary%DBoundValInc, 2) /= m) then
               deallocate (obj%FEMDomain%Boundary%DBoundValInc)
            end if
         end if
         obj%FEMDomain%Boundary%DBoundValInc(:, :) = 1.0d0/dble(timestep)*obj%FEMDomain%Boundary%DBoundVal(:, :)

         !obj%FEMDomain%Boundary%DBoundVal(:,:)=obj%FEMDomain%Boundary%DBoundValInc(:,:)
      elseif (obj%FEMDomain%ControlPara%SimMode == 2) then
         n = size(obj%FEMDomain%Boundary%NBoundVal, 1)
         m = size(obj%FEMDomain%Boundary%NBoundVal, 2)
         if (.not. allocated(obj%FEMDomain%Boundary%NBoundValInc)) then
            allocate (obj%FEMDomain%Boundary%NBoundValInc(n, m))
         else
            if (size(obj%FEMDomain%Boundary%NBoundValInc, 1) /= n .or. &
                size(obj%FEMDomain%Boundary%NBoundValInc, 2) /= m) then
               deallocate (obj%FEMDomain%Boundary%DBoundValInc)
            end if
         end if

         obj%FEMDomain%Boundary%NBoundValInc(:, :) = 1.0d0/dble(timestep)*obj%FEMDomain%Boundary%NBoundVal(:, :)
         obj%FEMDomain%Boundary%NBoundVal(:, :) = obj%FEMDomain%Boundary%NBoundValInc(:, :)

      else
         print *, "ERROR :: Displacement Control or Force Control?"
      end if

   end subroutine
!########################## Initialize Boundary Conditions of Finite Deformation  ###################################

!########################## Initialize Boundary Conditions of Finite Deformation  ###################################
   subroutine UpdateBCInTimestep(obj)
      class(FiniteDeform_), intent(inout)::obj

      integer(int32) :: timestep

      if (obj%FEMDomain%ControlPara%SimMode == 1) then

         !obj%FEMDomain%Boundary%DBoundVal(:,:)=obj%FEMDomain%Boundary%DBoundValInc(:,:)

      elseif (obj%FEMDomain%ControlPara%SimMode == 2) then

         !obj%FEMDomain%Boundary%NBoundVal(:,:)=obj%FEMDomain%Boundary%NBoundValInc(:,:)
      else
         print *, "ERROR :: Displacement Control or Force Control?"
      end if

   end subroutine
!########################## Initialize Boundary Conditions of Finite Deformation  ###################################

!#############################################################
   subroutine ImportFEMDomainFiDe(obj, OptionalFileFormat, OptionalProjectName)
      class(FEMDomain_), intent(inout)::obj
      character*4, optional, intent(in)::OptionalFileFormat
      character(*), optional, intent(in)::OptionalProjectName

      character*4::FileFormat
      character*70::ProjectName
      character*74 ::FileName
      integer(int32), allocatable::IntMat(:, :)
      real(real64), allocatable::RealMat(:, :)
      integer(int32) :: fh, i, j, NumOfDomain, n, m, DimNum
      character*70 Msg

      call DeallocateFEMDomain(obj)

      fh = 104
      if (present(OptionalFileFormat)) then
         FileFormat = OptionalFileFormat
      else
         FileFormat = ".scf"
      end if

      if (present(OptionalProjectName)) then
         ProjectName = OptionalProjectName
      else
         ProjectName = "untitled"
      end if

      FileName = ProjectName//FileFormat
      obj%FileName = FileName
      obj%FilePath = "./"
!!print *, "Project : ",ProjectName
!!print *, "is Exported as : ",FileFormat," format"
!!print *, "File Name is : ",FileName

      open (fh, file=FileName, status="old")

      if (FileFormat == ".scf") then

         read (fh, *) NumOfDomain
         allocate (IntMat(NumOfDomain, 2))
         allocate (obj%Mesh%SubMeshNodFromTo(NumOfDomain, 3))
         allocate (obj%Mesh%SubMeshElemFromTo(NumOfDomain, 3))

         do i = 1, NumOfDomain
            obj%Mesh%SubMeshNodFromTo(i, 1) = i
            read (fh, *) obj%Mesh%SubMeshNodFromTo(i, 2), obj%Mesh%SubMeshNodFromTo(i, 3)
         end do

         do i = 1, NumOfDomain
            obj%Mesh%SubMeshElemFromTo(i, 1) = i
            read (fh, *) obj%Mesh%SubMeshElemFromTo(i, 3)
            if (i == 1) then
               obj%Mesh%SubMeshElemFromTo(i, 2) = 1
            else
               obj%Mesh%SubMeshElemFromTo(i, 2) = obj%Mesh%SubMeshElemFromTo(i - 1, 3) + 1
            end if
         end do

         read (fh, *) n, m
         DimNum = m
         allocate (obj%Mesh%NodCoord(n, m))
         call ImportArray(obj%Mesh%NodCoord, OptionalFileHandle=fh)
         call CopyArray(obj%Mesh%NodCoord, obj%Mesh%NodCoordInit)

         read (fh, *) n, m
         read (fh, *) obj%Mesh%ElemType
         obj%ShapeFunction%ElemType = obj%Mesh%ElemType
         allocate (obj%Mesh%ElemNod(n, m))
         allocate (obj%Mesh%ElemMat(n))
         call ImportArray(obj%Mesh%ElemNod, OptionalFileHandle=fh)
         do i = 1, n
            read (fh, *) obj%Mesh%ElemMat(i)
         end do
         read (fh, *) n, m
         allocate (obj%MaterialProp%MatPara(n, m))
         call ImportArray(obj%MaterialProp%MatPara, OptionalFileHandle=fh)

         !######### Dirichlet boundary conditions #################
         DimNum = size(obj%Mesh%NodCoord, 2)
         allocate (obj%Boundary%DBoundNum(DimNum))
         read (fh, *) obj%Boundary%DBoundNum(:)
         allocate (obj%Boundary%DBoundNodID(maxval(obj%Boundary%DBoundNum), size(obj%Boundary%DBoundNum)))
         allocate (obj%Boundary%DBoundVal(maxval(obj%Boundary%DBoundNum), size(obj%Boundary%DBoundNum)))

         obj%Boundary%DBoundNodID(:, :) = -1
         obj%Boundary%DBoundVal(:, :) = 0.0d0

         do i = 1, size(obj%Boundary%DBoundNum, 1)
            do j = 1, obj%Boundary%DBoundNum(i)
               read (fh, *) obj%Boundary%DBoundNodID(j, i)
            end do
            do j = 1, obj%Boundary%DBoundNum(i)
               read (fh, *) obj%Boundary%DBoundVal(j, i)
            end do
         end do
         !######### Dirichlet boundary conditions #################

         !######### Neumann boundary conditions #################
         read (fh, *) n
         if (n == 0) then
            allocate (obj%Boundary%NBoundNum(0))
            allocate (obj%Boundary%NBoundNodID(0, 0))
            allocate (obj%Boundary%NBoundVal(0, 0))
         else
            allocate (obj%Boundary%NBoundNum(DimNum))
            allocate (obj%Boundary%NBoundNodID(n, size(obj%Boundary%NBoundNum)))
            allocate (obj%Boundary%NBoundVal(n, size(obj%Boundary%NBoundNum)))
            obj%Boundary%NBoundNodID(:, :) = -1
            obj%Boundary%NBoundVal(:, :) = 0.0d0

            obj%Boundary%NBoundNum(:) = n
            do i = 1, n
               read (fh, *) m
               obj%Boundary%NBoundNodID(i, :) = m
            end do

            do i = 1, n
               read (fh, *) obj%Boundary%NBoundVal(i, :)
            end do
         end if
         !######### Neumann boundary conditions #################

         !######### Initial conditions #################
         read (fh, *) n
         if (n == 0) then
            allocate (obj%Boundary%TBoundNum(0))
            allocate (obj%Boundary%TBoundNodID(0, 0))
            allocate (obj%Boundary%TBoundVal(0, 0))
         else
            allocate (obj%Boundary%TBoundNodID(n, 1))
            allocate (obj%Boundary%TBoundVal(n, 1))
            allocate (obj%Boundary%TBoundNum(1))

            obj%Boundary%TBoundNum = n

            if (n /= 0) then
               if (n < 0) then
                  print *, "ERROR :: number of initial conditions are to be zero"
               else
                  do i = 1, n
                     read (fh, *) obj%Boundary%TBoundNodID(i, 1)
                  end do
                  do i = 1, n
                     read (fh, *) obj%Boundary%TBoundVal(i, 1)
                  end do
               end if
            end if
         end if
         !######### Initial conditions #################

         read (fh, *) obj%ControlPara%SimMode, obj%ControlPara%ItrTol, obj%ControlPara%Timestep
         close (fh)
      else
    !!print *, "ERROR :: ExportFEMDomain >> only .scf file can be exported."
      end if

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

!#############################################################
   subroutine SetupFiniteDeform(obj, tol)
      class(FiniteDeform_), intent(inout)::obj
      integer(int32), optional, intent(in) :: tol

      if (obj%dt == 0.0d0 .or. obj%dt /= obj%dt) then
         obj%dt = 1.0d0
      end if
      !if(present(tol) )then
      !        obj%nr_tol = tol
      !else
      !        obj%nr_tol = 1.0e-08
      !endif

      call UpdateCurrConfig(obj)
      call GetDeformStressMatAndVector(obj)
   end subroutine
!#############################################################

!#############################################################
   subroutine UpdateFiniteDeform(obj, restart)
      class(FiniteDeform_), intent(inout)::obj
      logical, optional, intent(in) :: restart

      call UpdateCurrConfig(obj)

      call GetDeformStressMatAndVector(obj)

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

!#############################################################
   subroutine UpdateCurrConfig(obj, restart)
      class(FiniteDeform_), intent(inout)::obj
      logical, optional, intent(in) :: restart
      integer(int32) :: i, j, n, m
      integer(int32) :: num_node, num_elem, num_dim

      num_node = size(obj%FEMDomain%Mesh%NodCoord, 1)
      num_dim = size(obj%FEMDomain%Mesh%NodCoord, 2)
      num_elem = size(obj%FEMDomain%Mesh%ElemNod, 1)
      ! Displacement :: u
      if (.not. allocated(obj%DeformVecGloTot)) then
         allocate (obj%DeformVecGloTot(num_node*num_dim))
         obj%DeformVecGloTot(:) = 0.0d0
      end if
      ! ⊿u = v
      if (.not. allocated(obj%DeformVecGloInc)) then
         allocate (obj%DeformVecGloInc(num_node*num_dim))
         obj%DeformVecGloInc(:) = 0.0d0
      end if
      if (.not. allocated(obj%FEMDomain%Mesh%NodCoordInit)) then
         allocate (obj%FEMDomain%Mesh%NodCoordInit(num_node, num_dim))
         obj%FEMDomain%Mesh%NodCoordInit(:, :) = obj%FEMDomain%Mesh%NodCoord(:, :)
      end if

      !call showArraySIze(obj%FEMDomain%Mesh%NodCoord)
      !call showArraySIze(obj%FEMDomain%Mesh%NodCoordInit)
      !call showArraySize(obj%DeformVecGloInc)
      !call showArraySize(obj%DeformVecGloTot)
      !stop
      do i = 1, num_node
      do j = 1, num_dim
         obj%FEMDomain%Mesh%NodCoord(i, j) = obj%FEMDomain%Mesh%NodCoordInit(i, j) + &
                                             obj%DeformVecGloInc(num_dim*(i - 1) + j) + obj%DeformVecGloTot(num_dim*(i - 1) + j)
      end do
      end do

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

!#############################################################
   subroutine UpdateInitConfig(obj)
      class(FiniteDeform_), intent(inout)::obj

      integer(int32) :: i, j, n, m
      integer(int32) :: num_node, num_elem, num_dim

      num_node = size(obj%FEMDomain%Mesh%NodCoord, 1)
      num_dim = size(obj%FEMDomain%Mesh%NodCoord, 2)
      num_elem = size(obj%FEMDomain%Mesh%ElemNod, 1)
      if (.not. allocated(obj%DeformVecGloTot)) then
         allocate (obj%DeformVecGloTot(num_node*num_dim))
         obj%DeformVecGloTot(:) = 0.0d0
      end if
      if (.not. allocated(obj%DeformVecGloInc)) then
         allocate (obj%DeformVecGloInc(num_node*num_dim))
         obj%DeformVecGloInc(:) = 0.0d0
      end if

      do i = 1, num_node
         do j = 1, num_dim
            ! Updated Lagrangian
            obj%FEMDomain%Mesh%NodCoordInit(i, j) = obj%FEMDomain%Mesh%NodCoord(i, j)
         end do
      end do
   end subroutine
!#############################################################

!================================================================================
! �S�̍����}�g���N�X�̌v�Z
   subroutine GetDeformStressMatAndVector(obj, OptionalStep, restart)
      class(FiniteDeform_), intent(inout)::obj
      type(ConstModel_)        ::mdl
      logical, optional, intent(in) :: restart

      integer(int32), optional, intent(in) :: OptionalStep
      integer(int32) :: itr_rm, itr, itr_contact
      integer(int32) :: nod_num, dim_num, elemnod_num, elem_num
      real(real64), allocatable :: g(:, :), Bmat(:, :), Ce_neoHK(:, :), BTmat(:, :), Psymat(:, :), &
                                xymat(:, :), xymat_c(:, :), Jmat(:, :), s(:), Kmat_e(:, :), c_nod_coord(:, :), cc_nod_coord(:, :), &
                           dNdxi(:, :), F_iJ_n(:, :), F_iJ(:, :), C_IJ(:, :), Cp_IJ_n(:, :), Cp_IJ(:, :), b_ij(:, :), F_inv(:, :), &
                                   F_T_inv(:, :), F_T(:, :), dNdx(:, :), M_IJ(:, :), Cp_IJ_inv(:, :), gvec_e(:)
      integer(int32), allocatable::ij(:, :)
      integer(int32) i, j, k, m, p, q, gp_num, NumOfStrainMeasure
      real(real64) detJ, Lamda, mu, c, phy, psy, E, v, xx, Tol
      real(real64), allocatable::MatPara(:)

      itr_rm = obj%FEMDomain%ControlPara%ItrTol
      itr = obj%itr
      itr_contact = obj%FEMDomain%ControlPara%ItrTol
      Tol = obj%FEMDomain%ControlPara%Tol
      !if( abs(tol)<1.0e-15 )then
      tol = 1.0e-15
      !endif

      nod_num = size(obj%FEMDomain%Mesh%NodCoord, 1)
      elem_num = size(obj%FEMDomain%Mesh%ElemNod, 1)
      elemnod_num = size(obj%FEMDomain%Mesh%ElemNod, 2)
      dim_num = size(obj%FEMDomain%Mesh%NodCoord, 2)

      if (.not. allocated(obj%DeformVecEBEInc)) then
         allocate (obj%DeformVecEBEInc(elem_num, elemnod_num*dim_num))
         obj%DeformVecEBEInc(:, :) = 0.0d0
      end if
      if (.not. allocated(obj%DeformVecEBETot)) then
         allocate (obj%DeformVecEBETot(elem_num, elemnod_num*dim_num))
         obj%DeformVecEBETot = 0.0d0
      end if
      if (.not. allocated(obj%DeformStressRHS)) then
         allocate (obj%DeformStressRHS(elem_num, elemnod_num*dim_num))
         obj%DeformStressRHS(:, :) = 0.0d0
      end if
      if (.not. allocated(obj%DeformStressMat)) then
         allocate (obj%DeformStressMat(elem_num, elemnod_num*dim_num, elemnod_num*dim_num))
         obj%DeformStressMat(:, :, :) = 0.0d0
      end if
      if (.not. allocated(obj%VolInitCurrEBE)) then
         allocate (obj%VolInitCurrEBE(elem_num, 3)) !initial, current
         obj%VolInitCurrEBE(:, :) = 1.0d0
      end if

      obj%DeformStressRHS(:, :) = 0.0d0
      obj%DeformStressMat(:, :, :) = 0.0d0
      obj%FEMDomain%ShapeFunction%ElemType = obj%FEMDomain%Mesh%GetElemType()
      call SetShapeFuncType(obj%FEMDomain%ShapeFunction, ReducedIntegration=obj%ReducedIntegration)

      gp_num = obj%FEMDomain%ShapeFunction%NumOfGp
      if (dim_num <= 2) then
         NumOfStrainMeasure = 15
      else
         NumOfStrainMeasure = 49
      end if
      allocate (MatPara(7))
      if (.not. allocated(obj%DeformStrain)) then
         allocate (obj%DeformStrain(elem_num, gp_num, NumOfStrainMeasure))
         obj%DeformStrain(:, :, :) = 0.0d0
         obj%DeformStrain(:, :, 1) = 1.0d0  ! Cp_iJ_n(1,1)
         obj%DeformStrain(:, :, 2) = 1.0d0  ! Cp_iJ_n(2,2)
         obj%DeformStrain(:, :, 4) = 1.0d0  ! Cp_iJ(1,1)
         obj%DeformStrain(:, :, 5) = 1.0d0  ! Cp_iJ(2,2)
         obj%DeformStrain(:, :, 7) = 1.0d0  ! F_iJ_n(1,1)
         obj%DeformStrain(:, :, 8) = 1.0d0  ! F_iJ_n(2,2)
         obj%DeformStrain(:, :, 11) = 1.0d0  ! F_iJ(1,1)
         obj%DeformStrain(:, :, 12) = 1.0d0  ! F_iJ(2,2)
         if (NumOfStrainMeasure == 49) then
            obj%DeformStrain(:, :, 16) = 1.0d0 ! Cp_iJ_n(3,3)
            obj%DeformStrain(:, :, 19) = 1.0d0 ! Cp_iJ(3,3)
            obj%DeformStrain(:, :, 22) = 1.0d0 ! F_iJ_n(3,3)
            obj%DeformStrain(:, :, 27) = 1.0d0 ! F_iJ(3,3)
         end if

      end if
      !initialize
      if (.not. allocated(obj%DeformStress)) then
         allocate (obj%DeformStress(elem_num, gp_num, 6))
         obj%DeformStress(:, :, :) = 0.0d0
      end if
      if (.not. allocated(obj%DeformStressinc)) then
         allocate (obj%DeformStressinc(elem_num, gp_num, 6))
         obj%DeformStressinc(:, :, :) = 0.0d0
      end if
      m = elemnod_num*dim_num
      allocate (Kmat_e(elemnod_num*dim_num, elemnod_num*dim_num), gvec_e(elemnod_num*dim_num))
      obj%VolInitCurrEBE(:, 2) = 0.0d0
      do i = 1, elem_num !�v�f���ƃ��[�
         Kmat_e(:, :) = 0.0d0
         gvec_e(:) = 0.0d0
         E = obj%FEMDomain%MaterialProp%MatPara(obj%FEMDomain%Mesh%ElemMat(i), 1)
         v = obj%FEMDomain%MaterialProp%MatPara(obj%FEMDomain%Mesh%ElemMat(i), 2)

         ! allow direct-import of Youngs modulus and Poissons ratio
         if (allocated(obj%YoungsModulus)) then
            if (size(obj%YoungsModulus) /= elem_num) then
               print *, "ERROR :: FiniteDeform :: size(obj%YoungsModulus/=elem_num)"
            else
               E = obj%YoungsModulus(i)
            end if
         end if

         ! allow direct-import of Youngs modulus and Poissons ratio
         if (allocated(obj%PoissonsRatio)) then
            if (size(obj%PoissonsRatio) /= elem_num) then
               print *, "ERROR :: FiniteDeform :: size(obj%PoissonsRatio/=elem_num)"
            else
               v = obj%PoissonsRatio(i)
            end if
         end if
         mdl%Lamda = v*E/(1.0d0 + v)/(1.0d0 - 2.0d0*v)
         mdl%mu = E/2.0d0/(1.0d0 + v)
         c = obj%FEMDomain%MaterialProp%MatPara(obj%FEMDomain%Mesh%ElemMat(i), 4)
         phy = obj%FEMDomain%MaterialProp%MatPara(obj%FEMDomain%Mesh%ElemMat(i), 5)
         psy = obj%FEMDomain%MaterialProp%MatPara(obj%FEMDomain%Mesh%ElemMat(i), 6)

         MatPara(1) = E
         MatPara(2) = v
         MatPara(3) = Lamda
         MatPara(4) = mu
         MatPara(5) = c
         MatPara(6) = phy
         MatPara(7) = psy

         do j = 1, obj%FEMDomain%ShapeFunction%NumOfGp !�K�E�X�ϕ����ƃ��[�v
            ! -----J�}�g���N�X�̌v�Z-----------------------------------------
            call GetAllShapeFunc(obj%FEMDomain%ShapeFunction, elem_id=i, nod_coord=obj%FEMDomain%Mesh%NodCoord, &
                                 elem_nod=obj%FEMDomain%Mesh%ElemNod, OptionalGpID=j, ReducedIntegration=obj%ReducedIntegration)

            if (obj%infinitesimal .eqv. .false.) then
               ! ---Compute stress/strain measures of the Finite Elasto-Plast-
               call F_tensor_ICU(obj, i, j, mdl%F_iJ_n, mdl%F_iJ)

               call trans_rank_2(mdl%F_iJ, mdl%F_T)

               call C_tensor(mdl%F_iJ, mdl%C_IJ, mdl%b_ij, itr, dim_num)

               !call Cp_tensor(i,j,obj%DeformStrain,mdl%Cp_IJ_n,mdl%Cp_IJ,mdl%Cp_IJ_inv,dim_num)

               call inverse_rank_2(mdl%F_iJ, mdl%F_inv_iJ)
               call inverse_rank_2(mdl%F_T, mdl%F_T_inv_iJ)

               !call M_neo_Hookean(mdl%C_IJ,mdl%Cp_IJ,mdl%Cp_IJ_inv,mdl%M_IJ,mdl%Lamda,mdl%mu,i,j)

               !----Return-Mapping algorithm----------------------------------

               !call Return_Mapping_MCDP(dim_num,i,j,mdl%C_IJ,mdl%Cp_IJ,mdl%Cp_IJ_n,mdl%Cp_IJ_inv,mdl%M_IJ,MatPara,&
               !itr_rm,tol,obj%DeformStress,mdl%F_T,mdl%F_T_inv_ij,itr,itr_contact,obj%DeformStrain,OptionalStep)

               ! -----current config. Ce matrix-----------------------------------------
               !call Ce_neoHK_current(dim_num, i,j,mdl%Lamda,mdl%mu,mdl%C_IJ,mdl%Cp_IJ,mdl%b_ij,mdl%M_IJ,Ce_neoHK,mdl%F_T,mdl%F_T_inv,ij)

               ! -----B,W�}�g���N�X�̌v�Z----------------------------------------
               call B_mat(dim_num, obj%FEMDomain%ShapeFunction%dNdgzi, &
                          obj%FEMDomain%ShapeFunction%Jmat, obj%FEMDomain%ShapeFunction%detJ, mdl%Bmat, m)

               ! -----BT�}�g���N�X�̌v�Z-----------------------------------------
               !call trans_rank_2(Bmat,BTmat)
               !-----2-D neo-Hookean current stifness matrix-------------------

               ! Get Volumetric change
               obj%VolInitCurrEBE(i, 2) = obj%VolInitCurrEBE(i, 2) + &
                                          det_mat(mdl%F_iJ, size(mdl%F_iJ, 1))/dble(obj%FEMDomain%ShapeFunction%NumOfGp)
               mdl%ModelType = "NeoHookean"
               call GetKmat(obj, mdl, obj%FEMDomain%ShapeFunction, Kmat_e, gvec_e, dim_num, elemnod_num, i, j)
            else
               ! infinitesimal strain theory

               ! ---Compute stress/strain measures of the Finite Elasto-Plast-
               call F_tensor_ICU(obj, i, j, mdl%F_iJ_n, mdl%F_iJ)
               ! -----B,W�}�g���N�X�̌v�Z----------------------------------------
               call B_mat(dim_num, obj%FEMDomain%ShapeFunction%dNdgzi, &
                          obj%FEMDomain%ShapeFunction%Jmat, obj%FEMDomain%ShapeFunction%detJ, mdl%Bmat, m)
               call GetKmat(obj, mdl, obj%FEMDomain%ShapeFunction, Kmat_e, gvec_e, dim_num, elemnod_num, i, j)

            end if
            !do k=1,size(Kmat_e,1)
            !        write(200,*) Kmat_e(k,:),"|",gvec_e(k)
            !enddo
            !stop "debug"
            !call K_mat_e(j,obj%FEMDomain%ShapeFunction%GaussIntegWei, &
            !        BTmat, Ce_neoHK, Bmat, obj%FEMDomain%ShapeFunction%detJ, Kmat_e,F_iJ)
            !call GetGvec(obj,mdl,obj%FEMDomain%ShapeFunction,gvec_e,dim_num,elemnod_num,i)
            !call g_vector_e(i,j,obj%FEMDomain%ShapeFunction%GaussIntegWei,&
            !  BTmat,obj%DeformStress, obj%FEMDomain%ShapeFunction%detJ, gvec_e)

            ! -----�v�f�����}�g���N�X�̑�������--------------------------------
         end do

         ! ------�S�̍����}�g���N�X�ւ̏d�ˍ��킹------------------------------

         call K_mat_ICU(obj%DeformStressMat, obj%FEMDomain%Mesh%ElemNod, i, Kmat_e)

         call g_vector_ICU(i, obj%FEMDomain%Mesh%ElemNod, gvec_e, obj%DeformStressRHS)

      end do

   end subroutine
!================================================================================

!================================================================================
   subroutine GetKmat(obj, mdl, sf, Kmat_e, gvec_e, dim_num, elemnod_num, elem, gp)
      type(FiniteDeform_), intent(inout) :: obj
      type(ConstModel_), intent(inout)::mdl
      type(ShapeFunction_), intent(in)::sf
      real(real64), intent(inout) :: Kmat_e(:, :), gvec_e(:)
      integer(int32), intent(in)::dim_num, elemnod_num, elem, gp
      real(real64):: a0(3, 3), dXdgzi(3, 3)
      real(real64):: Xmat(elemnod_num, 3)
      real(real64), allocatable:: a0_inv(:, :), dgzidX(:, :), Jgzimat_inv(:, :), Dmat(:, :), Sigma(:), Sigma_i(:, :), &
                                  dumat_i(:, :), Sigma_tot(:, :)
      real(real64):: a1(elemnod_num, 3)
      real(real64):: K1(elemnod_num*dim_num, elemnod_num*dim_num)
      real(real64):: a2(elemnod_num, 3)
      real(real64):: A3(elemnod_num, elemnod_num)
      real(real64)::dumat(elemnod_num, 3)
      real(real64):: a4(elemnod_num, 3, elemnod_num, 3)
      real(real64)::Jgzimat(3, 3), dxxdgzi(3, 3), dNdX(elemnod_num, 3)

      integer(int32) ::alpha, beta, ganma, ii, I, J, jj, K, kk, L, ll, n, m, o, p, q, rr, R, s, loc_1, loc_2
      real(real64) :: detJgzi
      real(real64) :: PorePressure
      character*70 :: DerType
      !DerType="F_iJ"
      DerType = "c_current"

      allocate (dumat_i(elemnod_num*3, 1))
      !mdl%infinitesimal = obj%infinitesimal
      n = size(sf%ElemCoord, 1)
      m = size(sf%ElemCoord, 2)
      do i = 1, n !num of node per elems
         do j = 1, m ! num of dim
            dumat(i, j) = obj%DeformVecEBEInc(elem, m*(i - 1) + j)
            dumat_i((i - 1)*3 + j, 1) = dumat(i, j)
         end do
      end do
      !Xmat(:,:)   =sf%ElemCoord(:,:) !- dumat(:,:)
      dXdgzi(:, :) = matmul(sf%dNdgzi, sf%ElemCoord)
      !dXdgzi(:,:)=matmul(sf%dNdgzi,Xmat)
      call inverse_rank_2(dXdgzi, dgzidX)
      dNdX(:, :) = matmul(transpose(sf%dNdgzi), dgzidX)
      !a0(:,:)=matmul(sf%dNdgzi,sf%ElemCoord) !dxdgzi
      Jgzimat(:, :) = matmul(mdl%F_iJ, dXdgzi)
      call inverse_rank_2(Jgzimat, Jgzimat_inv)
      detJgzi = det_mat(Jgzimat, size(Jgzimat, 1))

      if (obj%infinitesimal .eqv. .false.) then

         !call inverse_rank_2(a0,a0_inv)
         call HyperElasticStress(mdl)
         call HyperElasticDer(mdl, DerType) !get cc mat
         call GetDmat(Dmat, mdl%StressDer, dim_num)
         call GetSigmaVec(Sigma, mdl%sigma, dim_num)

         K1(:, :) = 0.0d0
         A3(:, :) = matmul(matmul(dNdX, mdl%sigma), transpose(dNdX))
         do i = 1, elemnod_num
            do j = 1, elemnod_num
               K1((I - 1)*dim_num + 1, (J - 1)*dim_num + 1) = K1((I - 1)*dim_num + 1, (J - 1)*dim_num + 1) + A3(I, J)
               K1((I - 1)*dim_num + 2, (J - 1)*dim_num + 2) = K1((I - 1)*dim_num + 2, (J - 1)*dim_num + 2) + A3(I, J)
               K1((I - 1)*dim_num + 3, (J - 1)*dim_num + 3) = K1((I - 1)*dim_num + 3, (J - 1)*dim_num + 3) + A3(I, J)
            end do
         end do
         ! debug
         !print *, "size(obj%PorePressure)" ,size(obj%PorePressure)
         !write(3001,*) maxval(obj%PorePressure),minval(obj%PorePressure)
                !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         ! For water-soil coupling analysis

         if (.not. allocated(obj%PorePressure)) then
            PorePressure = 0.0d0
         elseif (size(obj%PorePressure) == 0) then
            PorePressure = 0.0d0
         else
            PorePressure = obj%PorePressure(elem)
         end if
         Sigma(1) = Sigma(1) - PorePressure/dble(obj%FEMDomain%ShapeFunction%NumOfGp)
         Sigma(2) = Sigma(2) - PorePressure/dble(obj%FEMDomain%ShapeFunction%NumOfGp)
         Sigma(3) = Sigma(3) - PorePressure/dble(obj%FEMDomain%ShapeFunction%NumOfGp)
                !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    Kmat_e(:, :) = Kmat_e(:, :) + matmul(matmul(transpose(mdl%Bmat), Dmat), mdl%Bmat)/det_mat(mdl%F_ij, size(mdl%F_ij, 1))*sf%detJ &
                        + K1(:, :)*sf%detJ
         !stop "debug l3"
         gvec_e(:) = gvec_e(:) + matmul(transpose(mdl%Bmat), Sigma)*sf%detJ
      else
         !call inverse_rank_2(a0,a0_inv)
         !call HyperElasticStress(mdl)
         !call HyperElasticDer(mdl,DerType) !get cc mat
         !call GetDmat(Dmat,mdl%StressDer,dim_num)
         !call GetSigmaVec(Sigma,mdl%sigma,dim_num)

         ! clear Dmat and import St. Venant stiffness matrix
         allocate (Dmat(6, 6))
         Dmat(:, :) = 0.0d0
         Dmat(1, 1) = 2.0d0*mdl%mu + mdl%lamda
         Dmat(1, 2) = mdl%lamda
         Dmat(1, 3) = mdl%lamda
         Dmat(2, 1) = mdl%lamda
         Dmat(2, 2) = 2.0d0*mdl%mu + mdl%lamda
         Dmat(2, 3) = mdl%lamda
         Dmat(3, 1) = mdl%lamda
         Dmat(3, 2) = mdl%lamda
         Dmat(3, 3) = 2.0d0*mdl%mu + mdl%lamda
         Dmat(4, 4) = mdl%mu
         Dmat(5, 5) = mdl%mu
         Dmat(6, 6) = mdl%mu

         ! stress dsigma = [D][B]{du}
         allocate (Sigma_i(6, 1))
         allocate (Sigma_tot(6, 1))
         Sigma_i(:, :) = matmul(Dmat, matmul(mdl%Bmat, dumat_i)) !\partial \sigma / \artial t * dt
         Sigma_tot(1:6, 1) = obj%DeformStress(elem, gp, 1:6)
         ! debug
         !print *, "size(obj%PorePressure)" ,size(obj%PorePressure)
         !write(3001,*) maxval(obj%PorePressure),minval(obj%PorePressure)
                !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         ! For water-soil coupling analysis

         if (.not. allocated(obj%PorePressure)) then
            PorePressure = 0.0d0
         elseif (size(obj%PorePressure) == 0) then
            PorePressure = 0.0d0
         else
            PorePressure = obj%PorePressure(elem)
         end if
         Sigma_tot(1, 1) = Sigma_tot(1, 1) - PorePressure!/dble(obj%FEMDomain%ShapeFunction%NumOfGp)
         Sigma_tot(2, 1) = Sigma_tot(2, 1) - PorePressure!/dble(obj%FEMDomain%ShapeFunction%NumOfGp)
         Sigma_tot(3, 1) = Sigma_tot(3, 1) - PorePressure!/dble(obj%FEMDomain%ShapeFunction%NumOfGp)
                !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

         Kmat_e(:, :) = Kmat_e(:, :) + matmul(matmul(transpose(mdl%Bmat), Dmat), mdl%Bmat)*sf%detJ
         !stop "debug l3"
         gvec_e(:) = gvec_e(:) + matmul(transpose(mdl%Bmat), Sigma_tot(:, 1))*sf%detJ &
                     + matmul(transpose(mdl%Bmat), Sigma_i(:, 1))*sf%detJ
         ! update d \sigma
         obj%DeformStressinc(elem, gp, 1:6) = Sigma_i(1:6, 1)
      end if

      !
      !do I=1,elemnod_num
      !        do ii=1,3
      !                do K=1,elemnod_num
      !                        do kk=1,3
      !                                do R=1,3
      !                                        do j=1,3
      !                                                a4(II,ii,K,kk)=a4(II,ii,K,kk)&
      !                                                +a1(I,j)*mdl%StressDer(ii,j,kk,R)*dNdX(K,R)*detJgzi
      !                                        enddo
      !                                enddo
      !                        enddo
      !                enddo
      !        enddo
      !enddo
!
      !do I=1,elemnod_num
      !        do ii=1,3
      !                do K=1,elemnod_num
      !                        do kk=1,3
      !                                loc_1=(I - 1)*3+ii
      !                                loc_2=(K - 1)*3+kk
      !                                a0(:,:)=matmul(mdl%F_inv_iJ,mdl%sigma)
      !                                a1(:,:)=matmul(dNdX,a0)
      !                                a2(:,:)=matmul(dNdX,mdl%F_inv_iJ)
      !                                Kmat_e(loc_1,loc_2)=Kmat_e(loc_1,loc_2)&
      !                                -a1(K,ii)*a2(I,kk)*detJgzi
!
      !                                a0(:,:)=matmul(mdl%F_inv_iJ,mdl%sigma)
      !                                a1(:,:)=matmul(dNdX,a0)
      !                                a2(:,:)=matmul(dNdX,dXdgzi)
      !                                a3(:,:)=matmul(a2,Jgzimat_inv)
      !                                Kmat_e(loc_1,loc_2)=Kmat_e(loc_1,loc_2)&
      !                                +a1(I,ii)*a3(K,kk)*detJgzi
      !
      !                                Kmat_e(loc_1,loc_2)=Kmat_e(loc_1,loc_2)&
      !                                +a4(II,ii,K,kk)*detJgzi
      !                        enddo
      !                enddo
      !        enddo
      !enddo
      !
      !do i=1,size(Kmat_e,1)
      !        write(10,*) Kmat_e(i,:)
      !enddo
      !write (10,*) " "
!        do alpha=1,elemnod_num
!                do ganma=1,3
!                        do I=1,elemnod_num
!                                do ii=1,3
!                                        do m=1,3
!                                                do j=1,3
!                                                        a4(alpha,ganma,I,ii)=a4(alpha,ganma,I,ii)&
!                                                        +a1(alpha,j)*a3(I,m)*mdl%StressDer(ii,j,ganma,m)*sf%detJ
!                                                enddo
!                                        enddo
!                                enddo
!                        enddo
!                enddo
!        enddo
!
!        do alpha=1,elemnod_num
!                do ganma=1,3
!                        do I=1,elemnod_num
!                                do ii=1,3
!                                        loc_1=(alpha-1)*3+ii
!                                        loc_2=(I - 1)*3+ganma
!                                        Kmat_e(loc_1,loc_2)=Kmat_e(loc_1,loc_2)&
!                                        -a1(alpha,ganma)*a2(I,ii)*sf%detJ&
!                                        +a2(alpha,ii)*a1(I,ganma)*sf%detJ&
!                                        -a4(alpha,ganma,I,ii)
!                                enddo
!                        enddo
!                enddo
!        enddo
!

   end subroutine
!=================================================================================

!================================================================================
   subroutine GetGvec(obj, mdl, sf, gvec_e, dim_num, elemnod_num, elem)
      type(FiniteDeform_), intent(in) :: obj
      type(ConstModel_), intent(inout)::mdl
      type(ShapeFunction_), intent(in)::sf
      real(real64), intent(inout) :: gvec_e(:)
      integer(int32), intent(in)::dim_num, elemnod_num, elem

      real(real64):: a0(3, 3), Jgzimat(3, 3), dXdgzi(3, 3)
      real(real64), allocatable:: a0_inv(:, :), dgzidX(:, :)
      real(real64):: a1(elemnod_num, 3)
      real(real64):: a2(elemnod_num, 3)
      real(real64):: dumat(elemnod_num, 3)
      real(real64):: Xmat(elemnod_num, 3)
      real(real64):: a3(elemnod_num, 3)
      real(real64):: dNdX(elemnod_num, 3)
      real(real64):: a4(elemnod_num, 3, elemnod_num, 3), detJgzi

      integer(int32) ::alpha, beta, ganma, ii, I, j, k, l, m, o, p, q, r, s, loc_1, n
      character*70 :: DerType
      DerType = "F_iJ"

      n = size(sf%ElemCoord, 1)
      m = size(sf%ElemCoord, 2)
      do i = 1, n !num of node per elems
         do j = 1, m ! num of dim
            dumat(i, j) = obj%DeformVecEBEInc(elem, m*(i - 1) + j)
         end do
      end do
      Xmat(:, :) = sf%ElemCoord(:, :) - dumat(:, :)
      dXdgzi(:, :) = matmul(sf%dNdgzi, Xmat)

      a0(:, :) = matmul(sf%dNdgzi, sf%ElemCoord)
      call inverse_rank_2(a0, a0_inv)
      call HyperElasticStress(mdl)
      call HyperElasticDer(mdl, DerType)
      a1(:, :) = matmul(transpose(sf%dNdgzi), a0)
      a2(:, :) = matmul(a1, mdl%sigma)
      a3(:, :) = matmul(a1, mdl%F_ij)
      a4(:, :, :, :) = 0.0d0
      Jgzimat(:, :) = matmul(mdl%F_iJ, dXdgzi)
      detJgzi = det_mat(Jgzimat, size(Jgzimat, 1))
      dXdgzi(:, :) = matmul(sf%dNdgzi, Xmat)
      call inverse_rank_2(dXdgzi, dgzidX)
      dNdX(:, :) = matmul(transpose(sf%dNdgzi), dgzidX)

      a1(:, :) = matmul(dNdX, mdl%F_inv_iJ)

      do I = 1, elemnod_num
         do ii = 1, 3
            do j = 1, 3
               loc_1 = (I - 1)*3 + ii
               gvec_e(loc_1) = gvec_e(loc_1) &
                               + a1(I, j)*mdl%sigma(ii, j)*detJgzi
            end do
         end do
      end do

   end subroutine

!=================================================================================

!=================================================================================
! K�}�g���N�X�ւ̏d�ˍ��킹
   subroutine K_mat_ICU(Kmat, elem_nod, i, Kemat)
      integer(int32), intent(in) :: i, elem_nod(:, :)
      real(real64), intent(in) :: Kemat(:, :)
      real(real64), intent(inout) :: Kmat(:, :, :)

      Kmat(i, :, :) = Kmat(i, :, :) + Kemat(:, :)

   end subroutine
!=================================================================================

!=================================================================================
   subroutine g_vector_ICU(elem, elem_nod, gvec_e, gvec)

      integer(int32), intent(in) :: elem_nod(:, :), elem
      real(real64), intent(in):: gvec_e(:)
      real(real64), intent(inout)::gvec(:, :)

      gvec(elem, :) = gvec_e(:)

   end subroutine
!=================================================================================

!=================================================================================
   subroutine F_tensor_ICU(obj, elem, gauss, F_iJ_n, F_iJ)
      class(FiniteDeform_), intent(inout)::obj
      integer(int32), intent(in)::elem, gauss

      real(real64), allocatable:: F_iJ_n(:, :), F_iJ(:, :), dNdx(:, :), x_dNdxi_inv(:, :), x_dNdxi(:, :), &
                                  dumat_t(:, :), f_n_n_1(:, :), dumat(:, :), x_u(:, :), tr_dNdgzi(:, :), dF_iJ(:, :)
      real(real64)::LungeKutta(3, 3), LungeKutta_1(3, 3), LungeKutta_2(3, 3), LungeKutta_3(3, 3), &
                     LungeKutta_4(3, 3), delta(3, 3), dt
      integer(int32) n, m, i, j

      n = size(obj%FEMDomain%ShapeFunction%ElemCoord, 1)
      m = size(obj%FEMDomain%ShapeFunction%ElemCoord, 2)
      if (.not. allocated(F_iJ_n)) allocate (F_iJ_n(3, 3))
      if (.not. allocated(f_n_n_1)) allocate (f_n_n_1(3, 3))
      if (.not. allocated(F_iJ)) allocate (F_iJ(3, 3))
      if (.not. allocated(dF_iJ)) then
         allocate (dF_iJ(3, 3))
         dF_iJ(:, :) = 0.0d0
      end if
      if (.not. allocated(dNdx)) allocate (dNdx(n, m))
      if (.not. allocated(x_dNdxi)) allocate (x_dNdxi(m, m))
      if (.not. allocated(x_dNdxi_inv)) allocate (x_dNdxi_inv(m, m))
      if (.not. allocated(dumat_t)) allocate (dumat_t(m, n))
      if (.not. allocated(dumat)) allocate (dumat(n, m))
      if (.not. allocated(x_u)) allocate (x_u(n, m))
      if (.not. allocated(tr_dNdgzi)) allocate (tr_dNdgzi(n, m))

      n = size(obj%FEMDomain%ShapeFunction%ElemCoord, 1)
      m = size(obj%FEMDomain%ShapeFunction%ElemCoord, 2)
      do i = 1, n !num of node per elems
         do j = 1, m ! num of dim
            dumat(i, j) = obj%DeformVecEBEInc(elem, m*(i - 1) + j)

         end do
      end do
      dumat_t(:, :) = transpose(dumat)
      delta(:, :) = 0.0d0
      delta(1, 1) = 1.0d0
      delta(2, 2) = 1.0d0
      delta(3, 3) = 1.0d0

      x_dNdxi(:, :) = 0.0d0
      tr_dNdgzi(:, :) = transpose(obj%FEMDomain%ShapeFunction%dNdgzi)
      x_u(:, :) = obj%FEMDomain%ShapeFunction%ElemCoord(:, :) - dumat(:, :)
      x_dNdxi(:, :) = transpose(matmul(obj%FEMDomain%ShapeFunction%dNdgzi, x_u)) !dxndgzi at tn

      call inverse_rank_2(x_dNdxi, x_dNdxi_inv)!dgzidxn at tn

      dNdx(:, :) = matmul(tr_dNdgzi, x_dNdxi_inv)

      if (m == 2) then
         f_n_n_1(:, :) = 0.0d0
         f_n_n_1(3, 3) = 1.0d0

         F_iJ(:, :) = 0.0d0
         F_iJ(3, 3) = 1.0d0

         F_iJ(1, 1) = 1.0d0
         F_iJ(1, 2) = 0.0d0
         F_iJ(2, 1) = 0.0d0
         F_iJ(2, 2) = 1.0d0

         F_iJ_n(:, :) = 0.0d0
         F_iJ_n(1, 1) = obj%DeformStrain(elem, gauss, 7)
         F_iJ_n(1, 2) = obj%DeformStrain(elem, gauss, 9)
         F_iJ_n(1, 3) = 0.0d0
         F_iJ_n(2, 1) = obj%DeformStrain(elem, gauss, 10)
         F_iJ_n(2, 2) = obj%DeformStrain(elem, gauss, 8)
         F_iJ_n(2, 3) = 0.0d0
         F_iJ_n(3, 1) = 0.0d0
         F_iJ_n(3, 2) = 0.0d0
         F_iJ_n(3, 3) = 1.0d0

         f_n_n_1(1:2, 1:2) = F_iJ(1:2, 1:2) + matmul(transpose(dumat), dNdx)
         F_iJ(:, :) = matmul(f_n_n_1, F_iJ_n)

         obj%DeformStrain(elem, gauss, 11) = F_iJ(1, 1)
         obj%DeformStrain(elem, gauss, 12) = F_iJ(2, 2)
         obj%DeformStrain(elem, gauss, 13) = F_iJ(1, 2)
         obj%DeformStrain(elem, gauss, 14) = F_iJ(2, 1)
      elseif (m == 3) then

         f_n_n_1(:, :) = 0.0d0

         F_iJ(:, :) = 0.0d0
         F_iJ(1, 1) = 1.0d0
         F_iJ(2, 2) = 1.0d0
         F_iJ(3, 3) = 1.0d0

         F_iJ_n(:, :) = 0.0d0
         F_iJ_n(1, 1) = obj%DeformStrain(elem, gauss, 7)
         F_iJ_n(1, 2) = obj%DeformStrain(elem, gauss, 9)
         F_iJ_n(1, 3) = obj%DeformStrain(elem, gauss, 23)
         F_iJ_n(2, 1) = obj%DeformStrain(elem, gauss, 10)
         F_iJ_n(2, 2) = obj%DeformStrain(elem, gauss, 8)
         F_iJ_n(2, 3) = obj%DeformStrain(elem, gauss, 24)
         F_iJ_n(3, 1) = obj%DeformStrain(elem, gauss, 25)
         F_iJ_n(3, 2) = obj%DeformStrain(elem, gauss, 26)
         F_iJ_n(3, 3) = obj%DeformStrain(elem, gauss, 22)

         dF_iJ(1, 1) = obj%DeformStrain(elem, gauss, 32)
         dF_iJ(1, 2) = obj%DeformStrain(elem, gauss, 33)
         dF_iJ(1, 3) = obj%DeformStrain(elem, gauss, 34)
         dF_iJ(2, 1) = obj%DeformStrain(elem, gauss, 35)
         dF_iJ(2, 2) = obj%DeformStrain(elem, gauss, 36)
         dF_iJ(2, 3) = obj%DeformStrain(elem, gauss, 37)
         dF_iJ(3, 1) = obj%DeformStrain(elem, gauss, 38)
         dF_iJ(3, 2) = obj%DeformStrain(elem, gauss, 39)
         dF_iJ(3, 3) = obj%DeformStrain(elem, gauss, 40)

         ! Clank-Nicolson
         f_n_n_1(:, :) = F_iJ(:, :) + 0.50d0*matmul(transpose(dumat), dNdx) + 0.50d0*dF_iJ(:, :)
         ! 4th order Lunge-Kutta
         !dt=1.0d0
         !f_n_n_1(:,:)=matmul(transpose(dumat),dNdx)
         !LungeKutta_1(:,:)=f_n_n_1(:,:)
         !LungeKutta_2(:,:)=f_n_n_1(:,:)+dt/2.0d0*matmul(f_n_n_1,LungeKutta_1)
         !LungeKutta_3(:,:)=f_n_n_1(:,:)+dt/2.0d0*matmul(f_n_n_1,LungeKutta_2)
         !LungeKutta_4(:,:)=f_n_n_1(:,:)+dt*matmul(f_n_n_1,LungeKutta_3)
         !LungeKutta(:,:)=dt/6.0d0*( LungeKutta_1(:,:)+2.0d0*LungeKutta_2(:,:)+2.0d0*LungeKutta_3(:,:)+LungeKutta_4(:,:) )
         !F_iJ(:,:)=F_iJ_n(:,:) + matmul(LungeKutta,F_iJ_n)
         F_iJ(:, :) = matmul(f_n_n_1, F_iJ_n)

         obj%DeformStrain(elem, gauss, 11) = F_iJ(1, 1)
         obj%DeformStrain(elem, gauss, 12) = F_iJ(2, 2)
         obj%DeformStrain(elem, gauss, 13) = F_iJ(1, 2)
         obj%DeformStrain(elem, gauss, 14) = F_iJ(2, 1)
         obj%DeformStrain(elem, gauss, 27) = F_iJ(3, 3)
         obj%DeformStrain(elem, gauss, 28) = F_iJ(1, 3)
         obj%DeformStrain(elem, gauss, 29) = F_iJ(2, 3)
         obj%DeformStrain(elem, gauss, 30) = F_iJ(3, 1)
         obj%DeformStrain(elem, gauss, 31) = F_iJ(3, 2)

         dF_iJ(:, :) = matmul(transpose(dumat), dNdx)
         obj%DeformStrain(elem, gauss, 41) = dF_iJ(1, 1)
         obj%DeformStrain(elem, gauss, 42) = dF_iJ(2, 2)
         obj%DeformStrain(elem, gauss, 43) = dF_iJ(1, 2)
         obj%DeformStrain(elem, gauss, 44) = dF_iJ(2, 1)
         obj%DeformStrain(elem, gauss, 45) = dF_iJ(3, 3)
         obj%DeformStrain(elem, gauss, 46) = dF_iJ(1, 3)
         obj%DeformStrain(elem, gauss, 47) = dF_iJ(2, 3)
         obj%DeformStrain(elem, gauss, 48) = dF_iJ(3, 1)
         obj%DeformStrain(elem, gauss, 49) = dF_iJ(3, 2)

      else
         print *, "FiniteDeformationClass/F_tensor_ICU :: >> Dimension", m, "is not Supported."
         stop
      end if

   end subroutine
!=================================================================================

!=================================================================================
   subroutine C_tensor(F, C_IJ, b_ij, itr, dim_num)
      real(real64), allocatable::C_IJ(:, :), b_ij(:, :), F_T(:, :)
      real(real64), intent(in)::F(:, :)
      integer(int32), intent(in)::itr, dim_num

      integer(int32) i

      if (.not. allocated(C_IJ)) allocate (C_IJ(3, 3))
      if (.not. allocated(b_ij)) allocate (b_ij(3, 3))

      call trans_rank_2(F, F_T)

      if (dim_num == 2) then
         C_IJ(:, :) = 0.0d0
         b_ij(:, :) = 0.0d0
         C_IJ(3, 3) = 1.0d0
         b_ij(3, 3) = 1.0d0

         C_IJ(1:2, 1:2) = matmul(F_T(1:2, 1:2), F(1:2, 1:2))
         b_ij(1:2, 1:2) = matmul(F(1:2, 1:2), F_T(1:2, 1:2))
      elseif (dim_num == 3) then

         C_IJ(:, :) = matmul(F_T, F)
         b_ij(:, :) = matmul(F, transpose(F))

      else
         print *, "C_tensor :: dimension : ", dim_num, "is not supported"
         stop
      end if

   end subroutine C_tensor
!=================================================================================
   subroutine Cp_tensor(elem, gauss, strain_measure, Cp_IJ_n, Cp_IJ, Cp_IJ_inv, dim_num)
      real(real64), allocatable:: Cp_iJ_n(:, :), Cp_iJ(:, :), Cp_IJ_inv(:, :)
      real(real64), intent(in)::strain_measure(:, :, :)
      integer(int32), intent(in)::elem, gauss, dim_num

      integer(int32) i

      if (.not. allocated(Cp_iJ_n)) allocate (Cp_iJ_n(3, 3))
      if (.not. allocated(Cp_iJ)) allocate (Cp_iJ(3, 3))

      Cp_iJ(:, :) = 0.0d0
      Cp_iJ_n(:, :) = 0.0d0
      Cp_iJ(3, 3) = 1.0d0
      Cp_iJ_n(3, 3) = 1.0d0

      Cp_iJ(1, 1) = strain_measure(elem, gauss, 4)
      Cp_iJ(1, 2) = strain_measure(elem, gauss, 6)
      Cp_iJ(2, 1) = strain_measure(elem, gauss, 6)
      Cp_iJ(2, 2) = strain_measure(elem, gauss, 5)

      Cp_iJ_n(1, 1) = strain_measure(elem, gauss, 1)
      Cp_iJ_n(1, 2) = strain_measure(elem, gauss, 3)
      Cp_iJ_n(2, 1) = strain_measure(elem, gauss, 3)
      Cp_iJ_n(2, 2) = strain_measure(elem, gauss, 2)

      if (dim_num == 3) then

         Cp_iJ_n(3, 3) = strain_measure(elem, gauss, 16)
         Cp_iJ_n(1, 3) = strain_measure(elem, gauss, 17)
         Cp_iJ_n(2, 3) = strain_measure(elem, gauss, 18)
         Cp_iJ_n(3, 1) = strain_measure(elem, gauss, 17)
         Cp_iJ_n(3, 2) = strain_measure(elem, gauss, 18)

         Cp_iJ(3, 3) = strain_measure(elem, gauss, 19)
         Cp_iJ(1, 3) = strain_measure(elem, gauss, 20)
         Cp_iJ(2, 3) = strain_measure(elem, gauss, 21)
         Cp_iJ(3, 1) = strain_measure(elem, gauss, 20)
         Cp_iJ(3, 2) = strain_measure(elem, gauss, 21)

      end if

      call inverse_rank_2(Cp_IJ, Cp_IJ_inv)

   end subroutine Cp_tensor
!=================================================================================
   subroutine M_neo_Hookean(C_IJ, Cp_IJ, Cp_IJ_inv, M_IJ, Lamda, mu, elem, gauss)
      real(real64), intent(in)::C_IJ(:, :), Lamda, mu, Cp_IJ(:, :), Cp_IJ_inv(:, :)
      real(real64), allocatable::M_IJ(:, :), G_IJ(:, :), C_Cp_1(:, :)
      integer(int32), intent(in):: elem, gauss
      integer(int32) i, j, n

      if (.not. allocated(M_IJ)) allocate (M_IJ(3, 3))

      allocate (C_Cp_1(3, 3), G_IJ(3, 3)) !3-D for constitutive equations

      n = size(C_IJ, 1)

      !Metric tensor
      G_IJ(:, :) = 0.0d0
      G_IJ(1, 1) = 1.0d0
      G_IJ(2, 2) = 1.0d0
      G_IJ(3, 3) = 1.0d0

      C_Cp_1(:, :) = matmul(C_IJ, Cp_IJ_inv)

      M_IJ(:, :) = Lamda/2.0d0*(det_mat(C_IJ, n)/det_mat(Cp_IJ, n) - 1.0d0)*G_IJ(:, :) + mu*(C_Cp_1(:, :) - G_IJ(:, :))

   end subroutine M_neo_Hookean
!=================================================================================
   subroutine Return_Mapping_MCDP(dim_num, elem, gauss, C_IJ, Cp_IJ, Cp_IJ_n, Cp_IJ_inv, M_IJ, MatPara, &
                                  itr_rm, tol, sigma, F_T, F_T_inv, itr, itr_contact, strain_measure, step)
      real(real64), intent(in)::C_IJ(:, :), Cp_IJ_n(:, :), F_T(:, :), F_T_inv(:, :), MatPara(:)
      real(real64), intent(inout)::Cp_IJ(:, :), sigma(:, :, :), strain_measure(:, :, :)
      real(real64), allocatable, intent(inout)::M_IJ(:, :), Cp_IJ_inv(:, :)
      integer(int32), intent(in)::elem, gauss, itr, itr_rm, step, itr_contact, dim_num

real(real64), allocatable::G_IJ(:, :), Jmat(:, :), Xvec(:), Yvec(:), dXvec(:), Zmat(:, :), sigm(:, :), M_FT(:, :), C_IJ_inv(:, :), &
    M_1(:,:),M_2(:,:),M_2_T(:,:),M_3(:,:),M_4(:,:),M_5(:,:),C_1(:,:),C_2(:,:),C_3(:,:),B_6(:,:),B_7(:,:),B_8(:,:),B_9(:,:),B_10(:,:)
  real(real64) detF,I1_M,J2_M,J3_M,Theta_M,BI_MC,BI_DP,fc_MC,er,tol,residual_0,yield_function_mc,detC,detCp,alpha,beta,gunma,omega,&
         a11,a12,a13,a14,a15,a16,a17,a18,a21,a22,a23,a24,a25,a26,MM,xx,c1,c2,c3,c4,c5,c6,c7,val,E,v,Lamda,mu,c,phy,psy
      integer(int32) n, A, B, I, J, K, L, R, M, P, Q, nn, itrmax, retmap_itr, mesh

      allocate (G_IJ(3, 3), Jmat(4, 4), Xvec(4), dXvec(4), Yvec(4), Zmat(3, 3), sigm(3, 3), M_FT(3, 3))
      allocate (M_1(3, 3), M_2(3, 3), M_2_T(3, 3), M_3(3, 3), M_4(3, 3), M_5(3, 3), C_1(3, 3), C_2(3, 3), C_3(3, 3), &
                B_6(3, 3), B_7(3, 3), B_8(3, 3), B_9(3, 3), B_10(3, 3))

      E = MatPara(1)
      v = MatPara(2)
      Lamda = MatPara(3)
      mu = MatPara(4)
      c = MatPara(5)
      phy = MatPara(6)
      psy = MatPara(7)
      !compute vriables for R-M algorithm

      G_IJ(:, :) = 0.0d0
      G_IJ(1, 1) = 1.0d0
      G_IJ(2, 2) = 1.0d0
      G_IJ(3, 3) = 1.0d0

      n = size(C_IJ, 1)
      detF = sqrt(det_mat(C_IJ, n))
      I1_M = M_iJ(1, 1) + M_iJ(2, 2) + M_iJ(3, 3)

      J2_M = 0.0d0
      do i = 1, 3
         do j = 1, 3

            J2_M = J2_M + 1.0d0/2.0d0*(M_IJ(i, j) - I1_M/3.0d0*G_IJ(i, j))*(M_IJ(i, j) - I1_M/3.0d0*G_IJ(i, j))
         end do
      end do

      J3_M = 0.0d0
      do I = 1, 3
         do J = 1, 3
            do K = 1, 3
               J3_M = J3_M + 1.0d0/3.0d0*(M_IJ(i, j) - I1_M/3.0d0*G_IJ(i, j))*(M_IJ(j, k) - I1_M/3.0d0*G_IJ(j, k)) &
                      *(M_IJ(k, i) - I1_M/3.0d0*G_IJ(k, i))

            end do
         end do
      end do

      if (J2_M == 0.0d0) then
         Theta_M = 0.0d0
      else
         val = 1.50d0*sqrt(3.0d0)*J3_M/(J2_M*sqrt(J2_M))
         if (val < -1.0d0) then
            val = -1.0d0
         elseif (val > 1.0d0) then
            val = 1.0d0
         end if
         Theta_M = 1.0d0/3.0d0*asin(val)
      end if

      BI_MC = (1.0d0 + sin(phy))/(1.0d0 - sin(phy))
      BI_DP = (1.0d0 + sin(psy))/(1.0d0 - sin(psy))
      fc_MC = (2.0d0*c*cos(phy))/(1.0d0 - sin(phy))

      if ((1.0d0 - sin(phy)) == 0.0d0) stop "  ( 1.0d0-sin(phy) )  =0"

      Yield_function_MC = 1.0d0/detF*(sqrt(J2_M) + ((BI_MC - 1.0d0)*I1_M - 3.0d0*detF*fc_MC)/(3.0d0*(BI_MC + 1.0d0)*cos(Theta_M) + &
                                                                                          sqrt(3.0d0)*(BI_MC - 1.0d0)*sin(Theta_M)))
      !write(54,*)&
      !                sqrt(2.0d0/3.0d0)*sqrt(2.0d0*J2_M)*sin(Theta_M + 2.0d0/3.0d0*3.141590d0) + I1_M/3.0d0,&
      !                sqrt(2.0d0/3.0d0)*sqrt(2.0d0*J2_M)*sin(Theta_M                       ) + I1_M/3.0d0,&
      !                sqrt(2.0d0/3.0d0)*sqrt(2.0d0*J2_M)*sin(Theta_M-2.0d0/3.0d0*3.141590d0) + I1_M/3.0d0
      !Plastic corrector algorithm

      !print *, sqrt(J2_M),Yield_function_MC,(BI_MC-1.0d0),I1_M,3.0d0*fc_MC

      !print *, Yield_function_MC
      !write(109,*) step,itr,Yield_function_MC,detF,I1_M,fc_Mc

      if (Yield_function_MC > 0.0d0 .and. itr_contact >= 3) then

         !initialization
         retmap_itr = 1

         Xvec(1) = Cp_IJ_n(1, 1)
         Xvec(2) = Cp_IJ_n(2, 2)
         Xvec(3) = Cp_IJ_n(1, 2)

         Xvec(4) = 0.0d0

         Yvec(1) = 0.0d0
         Yvec(2) = 0.0d0
         Yvec(3) = 0.0d0
         Yvec(4) = yield_function_mc

         call inverse_rank_2(C_IJ, C_IJ_inv)

         !Loop for Return-Mapping algorithm
         do
            dXvec(:) = 0.0d0
            Zmat(:, :) = 0.0d0
            Jmat(:, :) = 0.0d0

            I1_M = M_iJ(1, 1) + M_iJ(2, 2) + M_iJ(3, 3)

            J2_M = 0.0d0
            do i = 1, 3
               do j = 1, 3
                  J2_M = J2_M + 1.0d0/2.0d0*(M_IJ(i, j) - I1_M/3.0d0*G_IJ(i, j))*(M_IJ(i, j) - I1_M/3.0d0*G_IJ(i, j))
               end do
            end do

            J3_M = 0.0d0
            do I = 1, 3
               do J = 1, 3
                  do K = 1, 3
                     J3_M = J3_M + 1.0d0/3.0d0*(M_IJ(i, j) - I1_M/3.0d0*G_IJ(i, j))*(M_IJ(j, k) - I1_M/3.0d0*G_IJ(j, k)) &
                            *(M_IJ(k, i) - I1_M/3.0d0*G_IJ(k, i))

                  end do
               end do
            end do

            MM = 0.0d0
            do i = 1, 3
               do j = 1, 3
                  MM = MM + M_IJ(I, J)*M_IJ(J, I)
               end do
            end do

            if (J2_M == 0.0d0) then
               Theta_M = 0.0d0
            else
               Theta_M = 1.0d0/3.0d0*asin(1.50d0*sqrt(3.0d0)*J3_M/(J2_M*sqrt(J2_M)))
               if (J2_M == 0.0d0) stop "J2_M=0"
            end if

      Yield_function_MC = 1.0d0/detF*(sqrt(J2_M) + ((BI_MC - 1.0d0)*I1_M - 3.0d0*detF*fc_MC)/(3.0d0*(BI_MC + 1.0d0)*cos(Theta_M) + &
                                                                                          sqrt(3.0d0)*(BI_MC - 1.0d0)*sin(Theta_M)))
            Yvec(4) = yield_function_mc
            if (detF == 0.0d0) stop "  detF   =0"

            !debug
            print *, "YC", yield_function_mc

            if (retmap_itr == 1) then
               residual_0 = sqrt(dot_product(Yvec, Yvec))
            end if

            print *, "itr, residual R-M", retmap_itr, sqrt(dot_product(Yvec, Yvec))

            detC = det_mat(C_IJ, n)
            detCp = det_mat(Cp_IJ, n)

            !update dYij/dCpkl
            a18 = -Lamda/2*detC/detCp

            c1 = -4.0d0*(1.0d0/J2_M*(MM - I1_M*I1_M/3.0d0) + ((BI_DP - 1.0d0)/(BI_DP + 2.0d0))**2.0d0)**(-3.0d0/2.0d0)
            c2 = 2.0d0*(1.0d0/J2_M*(MM - I1_M*I1_M/3.0d0) + ((BI_DP - 1.0d0)/(BI_DP + 2.0d0))**2.0d0)**(-1.0d0/2.0d0)
            c3 = -I1_M*I1_M/3.0d0*a18*xvec(4)*c1/J2_M/sqrt(J2_M)
            c4 = -mu/2.0d0*xvec(4)*c1/J2_M/sqrt(J2_M)
            c5 = xvec(4)*c1/sqrt(3.0d0)*(BI_DP - 1.0d0)/(BI_DP + 2.0d0)/J2_M
            c6 = c3 + a18*c5*(3.0d0 - I1_M)

            M_1(:, :) = M_IJ(:, :) - I1_M/3.0d0*G_IJ(:, :)

            !M_2(:,:)=matmul(matmul(C_IJ,M_IJ),Cp_IJ_inv)

            !M_3(:,:)=matmul(Cp_IJ_inv,M_2)

            M_4(:, :) = matmul(M_1, Cp_IJ_n)

            M_5(:, :) = matmul(matmul(Cp_IJ_inv, C_IJ), matmul(M_1, Cp_IJ_inv))

            C_1(:, :) = matmul(C_IJ, Cp_IJ_inv)

            C_2(:, :) = matmul(Cp_IJ_inv, Cp_IJ_n)

            C_3(:, :) = matmul(Cp_IJ_inv, matmul(C_IJ, Cp_IJ_inv))

            alpha = -0.50d0/J2_M/sqrt(J2_M)
            beta = 1.0d0/sqrt(J2_M)
            gunma = -xvec(4)*mu*c2

            c7 = c4 + alpha*gunma

            !in case of the edge
            if ((3.0d0*(BI_MC + 1.0d0)*cos(Theta_M) + sqrt(3.0d0)*(BI_MC - 1.0d0)*sin(Theta_M)) == 0.0d0) then
               stop "singular stress of Mohr-Coulomb"
            end if

            a11 = 0.50d0/detF/sqrt(J2_M)
            a12 = -(BI_MC - 1.0d0)/detF/(3.0d0*(BI_MC + 1.0d0)*cos(Theta_M) + sqrt(3.0d0)*(BI_MC - 1.0d0)*sin(Theta_M))
            a13 = -1.0d0/detF*((BI_MC - 1.0d0)*I1_M + 3.0d0*detF*fc_MC)/ &
                  (3.0d0*(BI_MC + 1.0d0)*cos(Theta_M) + sqrt(3.0d0)*(BI_MC - 1.0d0)*sin(Theta_M))**2.0d0 &
                  *(3.0d0*(BI_MC + 1.0d0)*sin(Theta_M) - sqrt(3.0d0)*(BI_MC - 1.0d0)*cos(Theta_M))
            a14 = sqrt(3.0d0)/2.0d0/sqrt(J2_M)/J2_M/ &
                  sqrt(1.0d0 - (27.0d0/4.0d0*J3_M*J3_M/J2_M/J2_M/J2_M)**2.0d0)
            a15 = -3.0d0*sqrt(3.0d0)/4.0d0/sqrt(J2_M)/J2_M/J2_M*J3_M/ &
                  sqrt(1.0d0 - (27.0d0/4.0d0*J3_M*J3_M/J2_M/J2_M/J2_M)**2.0d0)
            a16 = -2.0d0*I1_M/3.0d0
            a17 = 2.0d0*I1_M*I1_M/9.0d0
            a18 = -Lamda/2*detC/detCp

            if ((3.0d0*(BI_MC + 1.0d0)*cos(Theta_M) + sqrt(3.0d0)*(BI_MC - 1.0d0)*sin(Theta_M)) == 0.0d0) &
               stop "  (3.0d0*(BI_MC + 1.0d0)*cos(Theta_M) + sqrt(3.0d0)*(BI_MC-1.0d0)*sin(Theta_M) )  =0"

            a21 = -a13*a14
            a22 = a11 - a13*a15 - a13*a14*a16
            a23 = -I1_M/3.0d0*a11 - a12 + a13*a15*I1_M/3.0d0 + a13*a14/3.0d0*MM - a13*a14*a17

            a25 = a18*(a21*MM + a22*I1_M + 3.0d0*a23)

            B_6(:, :) = matmul(C_IJ, matmul(M_IJ, Cp_IJ_inv))
            B_7(:, :) = matmul(Cp_IJ_inv, B_6)
            B_8(:, :) = matmul(Cp_IJ_inv, matmul(C_IJ, Cp_IJ_inv))

            B_6(:, :) = matmul(Cp_IJ_inv, matmul(M_IJ, C_IJ))
            B_9(:, :) = matmul(B_6, Cp_IJ_inv)

            B_6(:, :) = matmul(Cp_IJ_inv, matmul(M_IJ, M_IJ))
            B_10(:, :) = matmul(matmul(B_6, C_IJ), Cp_IJ_inv)

            Jmat(:, :) = 0.0d0
            I = 1
            J = 1
            K = 1
            L = 1
            Jmat(1, 1) = 1.0d0 - c6*Cp_IJ_n(I, J)*Cp_IJ_inv(L, K) - c7*M_4(I, J)*M_5(K, L) + mu*c5*Cp_IJ_n(I, J)*M_5(K, L) &
                         - beta*gunma*C_1(I, K)*C_2(L, J) + 1.0d0/3.0d0*beta*gunma*C_3(L, K)*Cp_IJ_n(I, J)

            I = 1
            J = 1
            K = 2
            L = 2
            Jmat(1, 2) = -c6*Cp_IJ_n(I, J)*Cp_IJ_inv(L, K) - c7*M_4(I, J)*M_5(K, L) + mu*c5*Cp_IJ_n(I, J)*M_5(K, L) &
                         - beta*gunma*C_1(I, K)*C_2(L, J) + 1.0d0/3.0d0*beta*gunma*C_3(L, K)*Cp_IJ_n(I, J)

            I = 1
            J = 1
            K = 1
            L = 2
            Jmat(1, 3) = -c6*Cp_IJ_n(I, J)*Cp_IJ_inv(L, K) - c7*M_4(I, J)*M_5(K, L) + mu*c5*Cp_IJ_n(I, J)*M_5(K, L) &
                         - beta*gunma*C_1(I, K)*C_2(L, J) + 1.0d0/3.0d0*beta*gunma*C_3(L, K)*Cp_IJ_n(I, J)

            I = 2
            J = 2
            K = 1
            L = 1
            Jmat(2, 1) = -c6*Cp_IJ_n(I, J)*Cp_IJ_inv(L, K) - c7*M_4(I, J)*M_5(K, L) + mu*c5*Cp_IJ_n(I, J)*M_5(K, L) &
                         - beta*gunma*C_1(I, K)*C_2(L, J) + 1.0d0/3.0d0*beta*gunma*C_3(L, K)*Cp_IJ_n(I, J)

            I = 2
            J = 2
            K = 2
            L = 2
            Jmat(2, 2) = 1.0d0 - c6*Cp_IJ_n(I, J)*Cp_IJ_inv(L, K) - c7*M_4(I, J)*M_5(K, L) + mu*c5*Cp_IJ_n(I, J)*M_5(K, L) &
                         - beta*gunma*C_1(I, K)*C_2(L, J) + 1.0d0/3.0d0*beta*gunma*C_3(L, K)*Cp_IJ_n(I, J)

            I = 2
            J = 2
            K = 1
            L = 2
            Jmat(2, 3) = -c6*Cp_IJ_n(I, J)*Cp_IJ_inv(L, K) - c7*M_4(I, J)*M_5(K, L) + mu*c5*Cp_IJ_n(I, J)*M_5(K, L) &
                         - beta*gunma*C_1(I, K)*C_2(L, J) + 1.0d0/3.0d0*beta*gunma*C_3(L, K)*Cp_IJ_n(I, J)

            I = 1
            J = 2
            K = 1
            L = 1
            Jmat(3, 1) = -c6*Cp_IJ_n(I, J)*Cp_IJ_inv(L, K) - c7*M_4(I, J)*M_5(K, L) + mu*c5*Cp_IJ_n(I, J)*M_5(K, L) &
                         - beta*gunma*C_1(I, K)*C_2(L, J) + 1.0d0/3.0d0*beta*gunma*C_3(L, K)*Cp_IJ_n(I, J)

            I = 1
            J = 2
            K = 2
            L = 2
            Jmat(3, 2) = -c6*Cp_IJ_n(I, J)*Cp_IJ_inv(L, K) - c7*M_4(I, J)*M_5(K, L) + mu*c5*Cp_IJ_n(I, J)*M_5(K, L) &
                         - beta*gunma*C_1(I, K)*C_2(L, J) + 1.0d0/3.0d0*beta*gunma*C_3(L, K)*Cp_IJ_n(I, J)

            I = 1
            J = 2
            K = 1
            L = 2
            Jmat(3, 3) = 1.0d0 - c6*Cp_IJ_n(I, J)*Cp_IJ_inv(L, K) - c7*M_4(I, J)*M_5(K, L) + mu*c5*Cp_IJ_n(I, J)*M_5(K, L) &
                         - beta*gunma*C_1(I, K)*C_2(L, J) + 1.0d0/3.0d0*beta*gunma*C_3(L, K)*Cp_IJ_n(I, J)

            do i = 1, size(zmat, 1)
               do j = 1, size(zmat, 1)
                  zmat(i, j) = zmat(i, j) + (1.0d0/sqrt(J2_M)*(M_IJ(i, j) - I1_M/3.0d0*G_IJ(i, j)) + &
                                             2.0d0*(BI_DP - 1.0d0)/sqrt(3.0d0)/(BI_DP + 2.0d0)*G_IJ(i, j))
                  if ((BI_DP + 2.0d0) == 0.0d0) stop "  (BI_DP + 2.0d0)  =0"

               end do
            end do
            zmat(:, :) = 2.0d0*((MM - I1_M*I1_M/3.0d0)/J2_M + &
                                ((BI_DP - 1.0d0)/(BI_DP + 2.0d0))**2.0d0)**(-1.0d0/2.0d0)*zmat(:, :)

            do M = 1, 3
               Jmat(1, 4) = Jmat(1, 4) - Zmat(1, M)*Cp_IJ_n(M, 1)
               Jmat(2, 4) = Jmat(2, 4) - Zmat(2, M)*Cp_IJ_n(M, 2)
               Jmat(3, 4) = Jmat(3, 4) - Zmat(1, M)*Cp_IJ_n(M, 2)
            end do
            I = 1
            J = 1
            Jmat(4, 1) = a25*Cp_IJ_inv(J, I) - mu*(a21*B_10(J, I) + a22*B_7(I, J) + a23*B_8(J, I))

            I = 2
            J = 2
            Jmat(4, 2) = a25*Cp_IJ_inv(J, I) - mu*(a21*B_10(J, I) + a22*B_7(I, J) + a23*B_8(J, I))

            I = 1
            J = 2
            Jmat(4, 3) = a25*Cp_IJ_inv(J, I) - mu*(a21*B_10(J, I) + a22*B_7(I, J) + a23*B_8(J, I))

            Jmat(4, 4) = 0.0d0

            !control parameters
            !write(53,*)"I,yield_function_mc,norm",retmap_itr,yield_function_mc,sqrt(dot_product(Yvec,Yvec))

            ! stop "debug r-m"
            itrmax = itr_rm
            er = 1.0e-14
            nn = size(dXvec)
            print *, "detJmat=", det_mat(Jmat, 4), det_mat(zmat, 2)

            !do i=1,size(Jmat,1)
            !        write(52,*)Jmat(i,:)
            !enddo

            Yvec(1) = Cp_IJ(1, 1)
            Yvec(2) = Cp_IJ(2, 2)
            Yvec(3) = Cp_IJ(1, 2)
            Yvec(4) = yield_function_mc
            do K = 1, 3
               Yvec(1) = Yvec(1) - (G_IJ(1, K) + Xvec(4)*Zmat(1, K))*Cp_IJ_n(K, 1)
               Yvec(2) = Yvec(2) - (G_IJ(2, K) + Xvec(4)*Zmat(2, K))*Cp_IJ_n(K, 2)
               Yvec(3) = Yvec(3) - (G_IJ(1, K) + Xvec(4)*Zmat(1, K))*Cp_IJ_n(K, 2)
            end do

            !solve
            call gauss_jordan_pv(Jmat, dXvec, Yvec, nn)

            Xvec(:) = Xvec(:) - dXvec(:)

            !update variables
            Cp_IJ(:, :) = 0.0d0
            Cp_IJ(3, 3) = 1.0d0
            Cp_IJ(1, 1) = Xvec(1)
            Cp_IJ(2, 2) = Xvec(2)
            Cp_IJ(1, 2) = Xvec(3)
            Cp_IJ(2, 1) = Xvec(3)

            call inverse_rank_2(Cp_IJ, Cp_IJ_inv)
            call M_neo_Hookean(C_IJ, Cp_IJ, Cp_IJ_inv, M_IJ, Lamda, mu, elem, gauss)

            !Judgement of convergence
            if (sqrt(dot_product(Yvec, Yvec)) <= TOL) then
               !converged
               if (residual_0 == 0.0d0) stop "  residual_0  =0"

               print *, "R-M converge itr=", itr
               exit
            elseif (retmap_itr >= itrmax) then
               print *, "Return-Mapping (E-P) did not converge"

               !debug: print yield surface
               xx = 0.40d0
               mesh = 20
               do i = 1, mesh

                  do j = 1, 12                !1
                     I1_M = xx*(1.0d0 - 2.0d0/dble(mesh)*dble(i))
                     Theta_M = 3.14159d0/6.0d0 - 3.14159d0/3.0d0/12.0d0*dble(j)
                     J2_M = (-(BI_MC - 1.0d0)*I1_M + 3.0d0*detF*fc_MC)/(3.0d0*(BI_MC + 1.0d0)*cos(Theta_M) + &
                                                                        sqrt(3.0d0)*(BI_MC - 1.0d0)*sin(Theta_M))

                     write (55, *) &
                        2.0d0/sqrt(3.0d0)*J2_M*sin(Theta_M + 2.0d0/3.0d0*3.141590d0) + I1_M/3.0d0, &
                        2.0d0/sqrt(3.0d0)*J2_M*sin(Theta_M) + I1_M/3.0d0, &
                        2.0d0/sqrt(3.0d0)*J2_M*sin(Theta_M - 2.0d0/3.0d0*3.141590d0) + I1_M/3.0d0
                  end do
                  do j = 1, 12                        !6
                     I1_M = xx*(1.0d0 - 2.0d0/dble(mesh)*dble(i))
                     Theta_M = -3.14159d0/6.0d0 + 3.14159d0/3.0d0/12.0d0*dble(j)
                     J2_M = (-(BI_MC - 1.0d0)*I1_M + 3.0d0*detF*fc_MC)/(3.0d0*(BI_MC + 1.0d0)*cos(Theta_M) + &
                                                                        sqrt(3.0d0)*(BI_MC - 1.0d0)*sin(Theta_M))
                     write (55, *) &
                        2.0d0/sqrt(3.0d0)*J2_M*sin(Theta_M + 2.0d0/3.0d0*3.141590d0) + I1_M/3.0d0, &
                        2.0d0/sqrt(3.0d0)*J2_M*sin(Theta_M - 2.0d0/3.0d0*3.141590d0) + I1_M/3.0d0, &
                        2.0d0/sqrt(3.0d0)*J2_M*sin(Theta_M) + I1_M/3.0d0
                  end do

                  do j = 1, 12                !3
                     I1_M = xx*(1.0d0 - 2.0d0/dble(mesh)*dble(i))
                     Theta_M = 3.14159d0/6.0d0 - 3.14159d0/3.0d0/12.0d0*dble(j)
                     J2_M = (-(BI_MC - 1.0d0)*I1_M + 3.0d0*detF*fc_MC)/(3.0d0*(BI_MC + 1.0d0)*cos(Theta_M) + &
                                                                        sqrt(3.0d0)*(BI_MC - 1.0d0)*sin(Theta_M))
                     write (55, *) &
                        2.0d0/sqrt(3.0d0)*J2_M*sin(Theta_M) + I1_M/3.0d0, &
                        2.0d0/sqrt(3.0d0)*J2_M*sin(Theta_M - 2.0d0/3.0d0*3.141590d0) + I1_M/3.0d0, &
                        2.0d0/sqrt(3.0d0)*J2_M*sin(Theta_M + 2.0d0/3.0d0*3.141590d0) + I1_M/3.0d0
                  end do

                  do j = 1, 12                        !5
                     I1_M = xx*(1.0d0 - 2.0d0/dble(mesh)*dble(i))
                     Theta_M = -3.14159d0/6.0d0 + 3.14159d0/3.0d0/12.0d0*dble(j)
                     J2_M = (-(BI_MC - 1.0d0)*I1_M + 3.0d0*detF*fc_MC)/(3.0d0*(BI_MC + 1.0d0)*cos(Theta_M) + &
                                                                        sqrt(3.0d0)*(BI_MC - 1.0d0)*sin(Theta_M))
                     write (55, *) &
                        2.0d0/sqrt(3.0d0)*J2_M*sin(Theta_M - 2.0d0/3.0d0*3.141590d0) + I1_M/3.0d0, &
                        2.0d0/sqrt(3.0d0)*J2_M*sin(Theta_M) + I1_M/3.0d0, &
                        2.0d0/sqrt(3.0d0)*J2_M*sin(Theta_M + 2.0d0/3.0d0*3.141590d0) + I1_M/3.0d0
                  end do

                  do j = 1, 12                        !4
                     I1_M = xx*(1.0d0 - 2.0d0/dble(mesh)*dble(i))
                     Theta_M = 3.14159d0/6.0d0 - 3.14159d0/3.0d0/12.0d0*dble(j)
                     J2_M = (-(BI_MC - 1.0d0)*I1_M + 3.0d0*detF*fc_MC)/(3.0d0*(BI_MC + 1.0d0)*cos(Theta_M) + &
                                                                        sqrt(3.0d0)*(BI_MC - 1.0d0)*sin(Theta_M))
                     write (55, *) &
                        2.0d0/sqrt(3.0d0)*J2_M*sin(Theta_M - 2.0d0/3.0d0*3.141590d0) + I1_M/3.0d0, &
                        2.0d0/sqrt(3.0d0)*J2_M*sin(Theta_M + 2.0d0/3.0d0*3.141590d0) + I1_M/3.0d0, &
                        2.0d0/sqrt(3.0d0)*J2_M*sin(Theta_M) + I1_M/3.0d0
                  end do
                  do j = 1, 12                !2
                     I1_M = xx*(1.0d0 - 2.0d0/dble(mesh)*dble(i))
                     Theta_M = -3.14159d0/6.0d0 + 3.14159d0/3.0d0/12.0d0*dble(j)
                     J2_M = (-(BI_MC - 1.0d0)*I1_M + 3.0d0*detF*fc_MC)/(3.0d0*(BI_MC + 1.0d0)*cos(Theta_M) + &
                                                                        sqrt(3.0d0)*(BI_MC - 1.0d0)*sin(Theta_M))
                     write (55, *) &
                        2.0d0/sqrt(3.0d0)*J2_M*sin(Theta_M) + I1_M/3.0d0, &
                        2.0d0/sqrt(3.0d0)*J2_M*sin(Theta_M + 2.0d0/3.0d0*3.141590d0) + I1_M/3.0d0, &
                        2.0d0/sqrt(3.0d0)*J2_M*sin(Theta_M - 2.0d0/3.0d0*3.141590d0) + I1_M/3.0d0
                  end do
               end do
               stop
            else
               retmap_itr = retmap_itr + 1
               cycle
            end if
         end do
      else
         !print *, "elastic"
      end if

      M_FT(:, :) = matmul(M_IJ, F_T)
      sigm(:, :) = 1.0d0/detF*matmul(F_T_inv, M_FT)

      if (dim_num == 2) then
         sigma(elem, gauss, 1) = sigm(1, 1)
         sigma(elem, gauss, 2) = sigm(2, 2)
         sigma(elem, gauss, 3) = sigm(1, 2)
         sigma(elem, gauss, 4) = sigm(3, 3)
      elseif (dim_num == 3) then

         sigma(elem, gauss, 1) = sigm(1, 1)
         sigma(elem, gauss, 2) = sigm(2, 2)
         sigma(elem, gauss, 3) = sigm(3, 3)
         sigma(elem, gauss, 4) = sigm(1, 2)
         sigma(elem, gauss, 5) = sigm(2, 3)
         sigma(elem, gauss, 6) = sigm(3, 1)
      else
         stop "Return_Mapping_MCDP >> dimension is to be 2 or 3"
      end if
      strain_measure(elem, gauss, 4) = Cp_IJ(1, 1)
      strain_measure(elem, gauss, 6) = Cp_IJ(1, 2)
      strain_measure(elem, gauss, 6) = Cp_IJ(2, 1)
      strain_measure(elem, gauss, 5) = Cp_IJ(2, 2)

      if (dim_num == 3) then
         strain_measure(elem, gauss, 19) = Cp_iJ(3, 3)
         strain_measure(elem, gauss, 20) = Cp_iJ(1, 3)
         strain_measure(elem, gauss, 21) = Cp_iJ(2, 3)
      end if

   end subroutine Return_Mapping_MCDP
!=================================================================================
! D�}�g���N�X�̌v�Z
   subroutine Ce_neoHK_current(dim_num, elem, gauss, Lame1, Lame2, C_IJ, Cp_IJ, b_ij, M_IJ, Ce_neoHK, F_T, F_T_inv, ij)
      integer(int32), intent(in) :: dim_num, elem, gauss
      real(real64), intent(in) :: Lame1, Lame2, C_IJ(:, :), Cp_IJ(:, :), b_ij(:, :), F_T(:, :), F_T_inv(:, :), M_IJ(:, :)
      real(real64), allocatable, intent(out) :: Ce_neoHK(:, :)
      integer(int32), allocatable, intent(out)::ij(:, :)
      real(real64), allocatable::t_ij(:, :), G_IJ(:, :)
      integer(int32) n, m, p, q, i, j, k, l
      real(real64) detF, detCp, detC

      allocate (G_IJ(3, 3))

      G_IJ(:, :) = 0.0d0
      G_IJ(1, 1) = 1.0d0
      G_IJ(2, 2) = 1.0d0
      G_IJ(3, 3) = 1.0d0
      if (allocated(ij)) deallocate (ij)
      if (dim_num == 2) then
         n = 3 ! �Ђ��݂���11,��22,��12��3����

         allocate (ij(3, 1))
         ij(1, 1) = 1; ij(1, 2) = 1; 
         ij(2, 1) = 2; ij(2, 2) = 2; 
         ij(3, 1) = 1; ij(3, 2) = 2; 
         m = size(C_IJ, 1)
         if (.not. allocated(Ce_neoHK)) allocate (Ce_neoHK(n, n))

         allocate (t_ij(n, n))

         detC = det_mat(C_IJ, m)
         detCp = det_mat(Cp_IJ, m)
         detF = sqrt(detC)

         t_ij = matmul(F_T_inv, matmul(M_IJ, F_T))

         do p = 1, 3
            do q = 1, 3
               i = ij(p, 1)
               j = ij(p, 2)
               k = ij(q, 1)
               l = ij(q, 2)

               Ce_neoHK(p, q) = Lame1*(detC/detCp*b_ij(i, j)*b_ij(l, k) &
                                       - Lame1*(detC/detCp - 1)*b_ij(i, k)*b_ij(l, j)) &
                                + 2.0d0*Lame2*b_ij(i, k)*b_ij(l, j) &
                                + G_IJ(i, k)*t_ij(l, j)
            end do
         end do
         Ce_neoHK(:, :) = 1.0d0/detF*Ce_neoHK(:, :)

      elseif (dim_num == 3) then
         n = 6 ! �Ђ��݂���11,��22,��12��3����
         allocate (ij(6, 2))
         ij(1, 1) = 1; ij(1, 2) = 1; 
         ij(2, 1) = 2; ij(2, 2) = 2; 
         ij(3, 1) = 3; ij(3, 2) = 3; 
         ij(4, 1) = 1; ij(4, 2) = 2; 
         ij(5, 1) = 2; ij(5, 2) = 3; 
         ij(6, 1) = 3; ij(6, 2) = 1; 
         m = size(C_IJ, 1)

         if (.not. allocated(Ce_neoHK)) allocate (Ce_neoHK(n, n))
         allocate (t_ij(n, n))

         detC = det_mat(C_IJ, m)
         detCp = det_mat(Cp_IJ, m)
         detF = sqrt(detC)

         t_ij = matmul(F_T_inv, matmul(M_IJ, F_T))

         do p = 1, 6
            do q = 1, 6
               i = ij(p, 1)
               j = ij(p, 2)
               k = ij(q, 1)
               l = ij(q, 2)

               Ce_neoHK(p, q) = Lame1*(detC/detCp*b_ij(i, j)*b_ij(l, k) &
                                       - Lame1*(detC/detCp - 1)*b_ij(i, k)*b_ij(l, j)) &
                                + 2.0d0*Lame2*b_ij(i, k)*b_ij(l, j) !+  G_IJ(i,k)*t_ij(l,j)
            end do
         end do
         Ce_neoHK(:, :) = 1.0d0/detF*Ce_neoHK(:, :)

      else
         print *, "Ce_neoHK_current >> Dimension, ", dim_num, " is not supported"
         stop
      end if

   end subroutine Ce_neoHK_current
!=================================================================================

!=================================================================================
! D�}�g���N�X�̌v�Z
   subroutine GetSigmaVec(Sigma, Sigma_ij, dim_num)
      real(real64), intent(in) :: Sigma_ij(:, :)
      real(real64), allocatable, intent(inout) :: Sigma(:)
      integer(int32), allocatable::ij(:, :)
      integer(int32), intent(in)::dim_num
      integer(int32) n, m, p, q, i, j, k, l
      real(real64) detF, detCp, detC

      if (allocated(ij)) deallocate (ij)
      if (dim_num == 2) then
         n = 3 ! �Ђ��݂���11,��22,��12��3����

         allocate (ij(3, 1))
         ij(1, 1) = 1; ij(1, 2) = 1; 
         ij(2, 1) = 2; ij(2, 2) = 2; 
         ij(3, 1) = 1; ij(3, 2) = 2; 
         if (.not. allocated(Sigma)) allocate (Sigma(n))
         do p = 1, 3
            i = ij(p, 1)
            j = ij(p, 2)
            Sigma(p) = Sigma_ij(i, j)

         end do

      elseif (dim_num == 3) then
         n = 6 ! �Ђ��݂���11,��22,��12��3����
         allocate (ij(6, 2))
         ij(1, 1) = 1; ij(1, 2) = 1; 
         ij(2, 1) = 2; ij(2, 2) = 2; 
         ij(3, 1) = 3; ij(3, 2) = 3; 
         ij(4, 1) = 1; ij(4, 2) = 2; 
         ij(5, 1) = 2; ij(5, 2) = 3; 
         ij(6, 1) = 3; ij(6, 2) = 1; 
         if (.not. allocated(Sigma)) allocate (Sigma(n))

         do p = 1, 6
            Sigma(p) = Sigma_ij(ij(p, 1), ij(p, 2))
         end do

      else
         print *, "Dmat >> Dimension, ", dim_num, " is not supported"
         stop
      end if

   end subroutine
!=================================================================================
!=================================================================================
! D�}�g���N�X�̌v�Z
   subroutine GetDmat(Dmat, c_ijkl, dim_num)
      real(real64), intent(in) :: c_ijkl(:, :, :, :)
      real(real64), allocatable, intent(inout) :: Dmat(:, :)
      integer(int32), allocatable::ij(:, :)
      integer(int32), intent(in)::dim_num
      integer(int32) n, m, p, q, i, j, k, l
      real(real64) detF, detCp, detC

      if (allocated(ij)) deallocate (ij)
      if (dim_num == 2) then
         n = 3 ! �Ђ��݂���11,��22,��12��3����

         allocate (ij(3, 1))
         ij(1, 1) = 1; ij(1, 2) = 1; 
         ij(2, 1) = 2; ij(2, 2) = 2; 
         ij(3, 1) = 1; ij(3, 2) = 2; 
         if (.not. allocated(Dmat)) allocate (Dmat(n, n))

         do p = 1, 3
            do q = 1, 3
               i = ij(p, 1)
               j = ij(p, 2)
               k = ij(q, 1)
               l = ij(q, 2)

               Dmat(p, q) = c_ijkl(i, j, k, l)
            end do
         end do

      elseif (dim_num == 3) then
         n = 6 ! �Ђ��݂���11,��22,��12��3����
         allocate (ij(6, 2))
         ij(1, 1) = 1; ij(1, 2) = 1; 
         ij(2, 1) = 2; ij(2, 2) = 2; 
         ij(3, 1) = 3; ij(3, 2) = 3; 
         ij(4, 1) = 1; ij(4, 2) = 2; 
         ij(5, 1) = 2; ij(5, 2) = 3; 
         ij(6, 1) = 3; ij(6, 2) = 1; 
         if (.not. allocated(Dmat)) allocate (Dmat(n, n))

         do p = 1, 6
            do q = 1, 6
               i = ij(p, 1)
               j = ij(p, 2)
               k = ij(q, 1)
               l = ij(q, 2)

               Dmat(p, q) = c_ijkl(i, j, k, l)
            end do
         end do

      else
         print *, "Dmat >> Dimension, ", dim_num, " is not supported"
         stop
      end if

   end subroutine
!=================================================================================
   !B�}�g���N�X�̌v�Z
   subroutine B_mat(dim_num, Psymat, Jmat, detJ, Bmat, mm)
      real(real64), intent(in) :: Psymat(:, :), Jmat(:, :), detJ ! J�̋t�s��
      real(real64), allocatable, intent(inout) :: Bmat(:, :)
      integer(int32), intent(in)::dim_num
      real(real64), allocatable :: JPsy(:, :), Jin(:, :)
      integer(int32) k, l, m, n, a, b, p, mm, i, j, q

      if (dim_num == 2) then
         k = 3
      elseif (dim_num == 3) then
         k = 6
      else
         stop "B_mat >> dim_num = tobe 2 or 3 "
      end if
      !k = size(ij,1)   ! �Ђ��݂���11,��22,��12��3����
      l = mm  ! 8�ړ_�v�f*2����
      m = size(Psymat, 1)
      n = size(Psymat, 2)

      if (allocated(Bmat)) deallocate (Bmat)
      allocate (Bmat(k, l))
      allocate (JPsy(m, n))
      allocate (Jin(m, m))

! J:Psymat�̌v�Z
      if (mm/2 == 3) then
         stop 'Now Constructing'
      elseif (mm/2 == 4) then
         if (detJ == 0.0d0) stop "Bmat,detJ=0"
         Jin(1, 1) = (1.0d0/detJ)*Jmat(2, 2)
         Jin(2, 2) = (1.0d0/detJ)*Jmat(1, 1)
         Jin(1, 2) = (-1.0d0/detJ)*Jmat(1, 2)
         Jin(2, 1) = (-1.0d0/detJ)*Jmat(2, 1)
         JPsy(:, :) = matmul(Jin, Psymat)

         Bmat(1, 1) = JPsy(1, 1)
         Bmat(1, 2) = 0.0d0
         Bmat(1, 3) = 1.0d0*JPsy(1, 2)
         Bmat(1, 4) = 0.0d0
         Bmat(1, 5) = 1.0d0*JPsy(1, 3)
         Bmat(1, 6) = 0.0d0
         Bmat(1, 7) = 1.0d0*JPsy(1, 4)
         Bmat(1, 8) = 0.0d0
         Bmat(2, 1) = 0.0d0
         Bmat(2, 2) = 1.0d0*JPsy(2, 1)
         Bmat(2, 3) = 0.0d0
         Bmat(2, 4) = 1.0d0*JPsy(2, 2)
         Bmat(2, 5) = 0.0d0
         Bmat(2, 6) = 1.0d0*JPsy(2, 3)
         Bmat(2, 7) = 0.0d0
         Bmat(2, 8) = 1.0d0*JPsy(2, 4)
         Bmat(3, 1) = Bmat(2, 2)
         Bmat(3, 2) = Bmat(1, 1)
         Bmat(3, 3) = Bmat(2, 4)
         Bmat(3, 4) = Bmat(1, 3)
         Bmat(3, 5) = Bmat(2, 6)
         Bmat(3, 6) = Bmat(1, 5)
         Bmat(3, 7) = Bmat(2, 8)
         Bmat(3, 8) = Bmat(1, 7)

      elseif (mm/2 == 8) then
         Jin(1, 1) = (1.0d0/detJ)*Jmat(2, 2)
         Jin(2, 2) = (1.0d0/detJ)*Jmat(1, 1)
         Jin(1, 2) = (-1.0d0/detJ)*Jmat(2, 1)
         Jin(2, 1) = (-1.0d0/detJ)*Jmat(1, 2)
         JPsy(:, :) = matmul(Jin, Psymat)

         Bmat(1, 1) = -JPsy(1, 1)
         Bmat(1, 2) = 0.0d0
         Bmat(1, 3) = -1.0d0*JPsy(1, 2)
         Bmat(1, 4) = 0.0d0
         Bmat(1, 5) = -1.0d0*JPsy(1, 3)
         Bmat(1, 6) = 0.0d0
         Bmat(1, 7) = -1.0d0*JPsy(1, 4)
         Bmat(1, 8) = 0.0d0
         Bmat(1, 9) = -1.0d0*JPsy(1, 5)
         Bmat(1, 10) = 0.0d0
         Bmat(1, 11) = -1.0d0*JPsy(1, 6)
         Bmat(1, 12) = 0.0d0
         Bmat(1, 13) = -1.0d0*JPsy(1, 7)
         Bmat(1, 14) = 0.0d0
         Bmat(1, 15) = -1.0d0*JPsy(1, 8)
         Bmat(1, 16) = 0.0d0
         Bmat(2, 1) = 0.0d0
         Bmat(2, 2) = -1.0d0*JPsy(2, 1)
         Bmat(2, 3) = 0.0d0
         Bmat(2, 4) = -1.0d0*JPsy(2, 2)
         Bmat(2, 5) = 0.0d0
         Bmat(2, 6) = -1.0d0*JPsy(2, 3)
         Bmat(2, 7) = 0.0d0
         Bmat(2, 8) = -1.0d0*JPsy(2, 4)
         Bmat(2, 9) = 0.0d0
         Bmat(2, 10) = -1.0d0*JPsy(2, 5)
         Bmat(2, 11) = 0.0d0
         Bmat(2, 12) = -1.0d0*JPsy(2, 6)
         Bmat(2, 13) = 0.0d0
         Bmat(2, 14) = -1.0d0*JPsy(2, 7)
         Bmat(2, 15) = 0.0d0
         Bmat(2, 16) = -1.0d0*JPsy(2, 8)
         Bmat(3, 1) = Bmat(2, 2)
         Bmat(3, 2) = Bmat(1, 1)
         Bmat(3, 3) = Bmat(2, 4)
         Bmat(3, 4) = Bmat(1, 3)
         Bmat(3, 5) = Bmat(2, 6)
         Bmat(3, 6) = Bmat(1, 5)
         Bmat(3, 7) = Bmat(2, 8)
         Bmat(3, 8) = Bmat(1, 7)
         Bmat(3, 9) = Bmat(2, 10)
         Bmat(3, 10) = Bmat(1, 9)
         Bmat(3, 11) = Bmat(2, 12)
         Bmat(3, 12) = Bmat(1, 11)
         Bmat(3, 13) = Bmat(2, 14)
         Bmat(3, 14) = Bmat(1, 13)
         Bmat(3, 15) = Bmat(2, 16)
         Bmat(3, 16) = Bmat(1, 15)
      elseif (k == 6) then

         if (detJ == 0.0d0) stop "Bmat,detJ=0"

         call inverse_rank_2(Jmat, Jin)

         JPsy(:, :) = transpose(matmul(transpose(Psymat), Jin)) !dNdgzi* dgzidx
         Bmat(:, :) = 0.0d0
         do q = 1, size(JPsy, 2)
            do p = 1, dim_num
               Bmat(p, dim_num*(q - 1) + p) = JPsy(p, q)
            end do
            Bmat(4, dim_num*(q - 1) + 1) = JPsy(2, q);
             Bmat(4, dim_num*(q - 1) + 2) = JPsy(1, q);
             Bmat(4, dim_num*(q - 1) + 3) = 0.0d0;
             
            Bmat(5, dim_num*(q - 1) + 1) = 0.0d0;
             Bmat(5, dim_num*(q - 1) + 2) = JPsy(3, q);
             Bmat(5, dim_num*(q - 1) + 3) = JPsy(2, q);
             
            Bmat(6, dim_num*(q - 1) + 1) = JPsy(3, q);
             Bmat(6, dim_num*(q - 1) + 2) = 0.0d0;
             Bmat(6, dim_num*(q - 1) + 3) = JPsy(1, q);
             
         end do

         Bmat(4:6, :) = 0.50d0*Bmat(4:6, :)

      else
         stop "Bmat >> The element is not supported."
      end if

   end subroutine B_mat
!=================================================================================
! �K�E�X�ϕ�
   subroutine K_mat_e(j, s, BTmat, Ce_neoHK, Bmat, detJ, Kmat_e, F_iJ)
      integer(int32), intent(in) :: j
      real(real64), intent(in) :: BTmat(:, :), Ce_neoHK(:, :), Bmat(:, :), detJ, s(:), F_iJ(:, :)
      real(real64), intent(out) :: Kmat_e(:, :)
      real(real64), allocatable :: DBmat(:, :)
      integer(int32) nm, e, n

      n = size(F_iJ, 1)
      nm = size(Bmat, 2)
      e = size(Bmat, 1)
      allocate (DBmat(e, nm))

      DBmat = matmul(Ce_neoHK, Bmat)
      Kmat_e(:, :) = Kmat_e(:, :) + s(j)*detJ*matmul(BTmat, DBmat) !/det_mat(F_iJ,n)

   end subroutine K_mat_e
!=================================================================================
   subroutine g_vector_e(elem, gauss, s, BTmat, sigma, detJ, gvec_e)

      integer(int32), intent(in) :: elem, gauss
      real(real64), intent(in) :: BTmat(:, :), sigma(:, :, :), detJ, s(:)
      real(real64), intent(inout) :: gvec_e(:)
      real(real64), allocatable :: sigm(:)
      integer(int32) nm, i

      if (size(BTmat, 2) == 3) then
         allocate (sigm(3))
         nm = size(BTmat, 1)

         do i = 1, size(BTmat, 2)
            sigm(i) = sigma(elem, gauss, i)
         end do
      elseif (size(BTmat, 2) == 6) then
         allocate (sigm(6))
         nm = size(BTmat, 1)

         do i = 1, size(BTmat, 2)
            sigm(i) = sigma(elem, gauss, i)
         end do
      else
         stop "Error :: g_vec_e :: Invalid size of Bmat"
      end if

      gvec_e(:) = gvec_e(:) + s(gauss)*detJ*matmul(BTmat, sigm)

   end subroutine g_vector_e
!=================================================================================

!######################## Solve DiffusionEq ########################
   subroutine SolveFiniteDeform(obj, OptionItr, Solvertype, nr_tol)
      class(FiniteDeform_), intent(inout)::obj
      integer(int32), optional, intent(in)::OptionItr
      character(*), optional, intent(in)::Solvertype
      real(real64), optional, intent(in) :: nr_tol
      character(:), allocatable ::solver, defaultsolver
      type(IO_) :: f

      real(real64), allocatable::Amat(:, :), bvec(:), xvec(:)
      real(real64)::val, er
      integer(int32) ::i, j, n, m, k, l, dim1, dim2, nodeid1, nodeid2, localid, itrmax, SetBC, int1, int2
      integer(int32) :: dim_num, node_num, elem_num, node_num_elmtl
      !obj%nr_tol=1.0e-08
      if (present(nr_tol)) then
         obj%nr_tol = nr_tol
      end if

      node_num = size(obj%FEMDomain%Mesh%NodCoord, 1)
      dim_num = size(obj%FEMDomain%Mesh%NodCoord, 2)
      elem_num = size(obj%FEMDomain%Mesh%ElemNod, 1)
      node_num_elmtl = size(obj%FEMDomain%Mesh%ElemNod, 2)
      if (present(OptionItr)) then
         ! Interation in Netwon's Loop >> Set all BCs for zero.
         if (OptionItr <= 0) then
            SetBC = 1
         elseif (OptionItr >= 1) then
            SetBC = 0
         else
            print *, "ERROR :: FiniteDeformationClass.f90 >> SolveDeformStress :: Please check OptionalItr"
         end if
      else
         SetBC = 1
      end if

      !if sorving initial value
      if (SetBC == 1) then
         obj%DeformVecEBEInc(:, :) = 0.0d0
         obj%DeformVecGloInc(:) = 0.0d0
      end if

      if (present(SolverType)) then
         solver = SolverType
      else
         solver = "BiCGSTAB"
      end if

      n = node_num
      m = dim_num
      allocate (Amat(n*m, n*m), bvec(n*m), xvec(n*m))
      Amat(:, :) = 0.0d0
      bvec(:) = 0.0d0

      !print *, "stop debug"
      !call showArray(Amat, Name="Amat.txt")
      !stop

      !===================================
      ! assemble matrix
      do i = 1, elem_num
         do j = 1, node_num_elmtl
            nodeid1 = obj%FEMDomain%Mesh%ElemNod(i, j)
            do k = 1, node_num_elmtl
               nodeid2 = obj%FEMDomain%Mesh%ElemNod(i, k)
               do dim1 = 1, dim_num
                  do dim2 = 1, dim_num
                     int1 = (dim_num*(j - 1)) + dim1
                     int2 = (dim_num*(k - 1)) + dim2
                     Amat(dim_num*(nodeid1 - 1) + dim1, dim_num*(nodeid2 - 1) + dim2) = &
                        Amat(dim_num*(nodeid1 - 1) + dim1, dim_num*(nodeid2 - 1) + dim2) + &
                        (obj%DeformStressMat(i, int1, int2))
                  end do
               end do
            end do
         end do
      end do
      !===================================

      call GetTractionVector(obj)

      call GetInternalVector(obj)
      call GetResidualVector(obj)
      bvec(:) = obj%ResidualVecGlo(:)

      !===================================
      ! introduce D.B.C
      do k = 1, size(obj%FEMDomain%Boundary%DBoundNum)
         do i = 1, size(obj%FEMDomain%Boundary%DBoundNodID, 1)
            nodeid1 = obj%FEMDomain%Boundary%DBoundNodID(i, k)
            if (nodeid1 < 1) then
               cycle
            else

               val = obj%FEMDomain%Boundary%DBoundValInc(i, k)
               do j = 1, size(bvec)
                  bvec(j) = bvec(j) - Amat(j, dim_num*(nodeid1 - 1) + k)*val*dble(SetBC)
               end do
               Amat(dim_num*(nodeid1 - 1) + k, :) = 0.0d0
               Amat(:, dim_num*(nodeid1 - 1) + k) = 0.0d0
               Amat(dim_num*(nodeid1 - 1) + k, dim_num*(nodeid1 - 1) + k) = 1.0d0
               bvec(dim_num*(nodeid1 - 1) + k) = val*dble(SetBC)
            end if
         end do
      end do
      !===================================

      obj%ResidualVecGlo(:) = bvec(:)

      itrmax = 1000
      er = 1.0e-15
      n = size(bvec)
      xvec(:) = 0.0d0

      if (solver == "GaussJordan") then
         print *, "Solver type :: GaussJordan"
         call gauss_jordan_pv(Amat, xvec, bvec, n)
      elseif (solver == "BiCGSTAB") then
         print *, "Solver type :: BiCGSTAB"
         call bicgstab_dirichlet(Amat, bvec, xvec, size(xvec), itrmax, er, &
                                 obj%FEMDomain%Boundary%DBoundNodID, obj%FEMDomain%Boundary%DBoundValInc, SetBC)
      else
         print *, "Solver Name :", solver
         print *, "Critical ERROR :: No solvers are selected in FiniteDeform_"
         stop
      end if

!        print *, size(bvec)
!        call showArray(Amat,name="Amat.txt")
!        call f%open("Bvec.txt")
!        do i=1,size(bvec)
!                writE(f%fh,*) bvec(:)
!        enddo
!        call f%close()
!
!        stop

      !=====================================
      ! Export Values to Element-by-Element form
      do i = 1, elem_num
         do j = 1, node_num_elmtl
            nodeid1 = obj%FEMDomain%Mesh%ElemNod(i, j)
            do dim1 = 1, dim_num
               obj%DeformVecEBEInc(i, dim_num*(j - 1) + dim1) = &
                  obj%DeformVecEBEInc(i, dim_num*(j - 1) + dim1) + &
                  (xvec((dim_num*(nodeid1 - 1)) + dim1))
            end do
         end do
      end do
      obj%DeformVecGloInc(:) = obj%DeformVecGloInc(:) + xvec(:)

      if (dot_product(obj%DeformVecGloInc, obj%DeformVecGloInc) == 0.0d0) then
         obj%error = dot_product(xvec, xvec)
      else
         obj%error = abs(1.0d0 - dot_product(obj%DeformVecGloInc - xvec, obj%DeformVecGloInc - xvec)/ &
                         dot_product(obj%DeformVecGloInc, obj%DeformVecGloInc))
      end if

      !=====================================

   end subroutine

!######################## Solve DiffusionEq ########################

!######################## Display Finite Deformation ########################
   subroutine resultFiniteDeform(obj, path, name, step)
      class(FiniteDeform_), intent(inout)::obj
      character(*), intent(in) :: Name, path
      integer(int32), optional, intent(in)::step
      real(real64), allocatable::DBCvec(:, :), DispVector(:, :)
      integer(int32) :: i, j, n, m, fstep, dim_num

      if (size(obj%FEMDomain%Mesh%NodCoord, 2) == 2) then
         ! 2-D condition
         if (present(step)) then
            fstep = step
         else
            fstep = 1
         end if

         call GmshExportStress(obj%FEMDomain, obj%DeformVecGloTot, obj%DeformStress, &
                               obj%DeformStrain, fstep, Name=path//"/"//name)
         call obj%getDispVector(DispVector)
         call GmshPlotVector(obj%FEMDomain, Vector=DispVector, Name=path//"/" &
                             //name, FieldName="Displacement", &
                             Step=step, withMsh=.true., NodeWize=.true., onlyDirichlet=.true.)

      elseif (size(obj%FEMDomain%Mesh%NodCoord, 2) == 3) then

         ! 3-D condition
         if (present(step)) then
            fstep = step
         else
            fstep = 1
         end if
         call GmshExportStress(obj%FEMDomain, obj%DeformVecGloTot, obj%DeformStress, &
                               obj%DeformStrain, step, Name=path//"/"//name)
         call obj%getDBCVector(DBCvec)
         call obj%FEMDomain%GmshPlotVector(Vector=DBCvec, Name=path//"/" &
                                           //name, Step=fstep, FieldName="DispBoundary", &
                                           NodeWize=.true.)
         !call obj%exportAsPly()

      end if
   end subroutine
!######################## Display Finite Deformation ########################

   subroutine exportAsPlyFiniteDeform(obj)
      class(FiniteDeform_), intent(inout) :: obj
      real(real64), allocatable :: scalar(:)
      integer(int32) :: n, m
      n = size(obj%FEMDomain%Mesh%NodCoord, 1)

      allocate (scalar(n))

   end subroutine

!######################## Display Finite Deformation ########################
   subroutine DisplayDeformStress(obj, DisplayMode, OptionalStep, Name, withDirichlet)
      class(FiniteDeform_), intent(inout)::obj
      character(*), optional, intent(in) :: Name, DisplayMode
      logical, optional, intent(in) :: withDirichlet
      integer(int32), optional, intent(in)::OptionalStep
      real(real64), allocatable::DBCvec(:, :), DispVector(:, :)
      integer(int32) :: i, j, n, m, step, dim_num

      if (.not. associated(obj%FEMDomain)) then
         return
      end if
      if (obj%FEMDomain%Mesh%empty() .eqv. .true.) then
         return
      end if
      if (size(obj%FEMDomain%Mesh%NodCoord, 2) == 2) then
         ! 2-D condition
         if (present(OptionalStep)) then
            step = OptionalStep
         else
            step = 1
         end if

         if (present(DisplayMode)) then
            if (DisplayMode == "Terminal" .or. DisplayMode == "terminal") then
               do i = 1, size(obj%DeformVecEBEInc, 1)
                  print *, obj%DeformVecEBEInc(i, :)
               end do
            elseif (DisplayMode == "gmsh" .or. DisplayMode == "Gmsh") then
               call GmshExportStress(obj%FEMDomain, obj%DeformVecGloTot, obj%DeformStress, obj%DeformStrain, step, Name=Name)
               call obj%getDispVector(DispVector)
               call GmshPlotVector(obj%FEMDomain, Vector=DispVector, Name=Name, FieldName="Displacement", &
                                   Step=step, withMsh=.true., NodeWize=.true., onlyDirichlet=.true.)

            elseif (DisplayMode == "gnuplot" .or. DisplayMode == "Gnuplot") then
               call GnuplotExportStress(obj%FEMDomain, obj%DeformVecGloTot, obj%DeformStress, obj%DeformStrain, step)
            else
               print *, "Invalid DisplayMode >> DisplayDiffusionEq "
               print *, "DisplayMode '", DisplayMode, "' is not defined"
            end if
            return
         end if
      elseif (size(obj%FEMDomain%Mesh%NodCoord, 2) == 3) then

         ! 3-D condition
         if (present(OptionalStep)) then
            step = OptionalStep
         else
            step = 1
         end if

         if (present(DisplayMode)) then
            if (DisplayMode == "Terminal" .or. DisplayMode == "terminal") then
               do i = 1, size(obj%DeformVecEBEInc, 1)
                  print *, obj%DeformVecEBEInc(i, :)
               end do
            elseif (DisplayMode == "gmsh" .or. DisplayMode == "Gmsh") then

               call GmshExportStress(obj%FEMDomain, obj%DeformVecGloTot, obj%DeformStress, obj%DeformStrain, step, Name=Name)
               if (present(withDirichlet)) then
                  if (withDirichlet .eqv. .true.) then
                     ! Export dirichlet (Deformation) boundary conditions
                     call obj%getDBCVector(DBCvec)
                   call obj%FEMDomain%GmshPlotVector(Vector=DBCvec, Name=Name, Step=step, FieldName="DispBoundary", NodeWize=.true.)
                  end if
               end if
            elseif (DisplayMode == "gnuplot" .or. DisplayMode == "Gnuplot") then
               call GnuplotExportStress(obj%FEMDomain, obj%DeformVecGloTot, obj%DeformStress, obj%DeformStrain, step)
            else
               print *, "Invalid DisplayMode >> DisplayDiffusionEq "
               print *, "DisplayMode '", DisplayMode, "' is not defined"
            end if
            return
         end if
      else
         stop "only 2D and 3D can be displayed"
      end if
   end subroutine
!######################## Display Finite Deformation ########################

!############# Get Traction Vector ######################
   subroutine GetTractionVector(obj)
      class(FiniteDeform_), intent(inout)::obj

      integer(int32) :: i, j, dim_num, node

      if (allocated(obj%TractionVecGlo) .and. size(obj%DeformVecGloTot) /= size(obj%TractionVecGlo)) then
         deallocate (obj%TractionVecGlo)
      end if
      if (.not. allocated(obj%TractionVecGlo)) then
         allocate (obj%TractionVecGlo(size(obj%DeformVecGloTot)))
      end if

      obj%TractionVecGlo(:) = 0.0d0
      if (allocated(obj%FEMDomain%Boundary%NBoundNodID)) then
         dim_num = size(obj%FEMDomain%Boundary%NBoundNodID, 2)
         if (dim_num == 0) then
            return
         end if

         do i = 1, size(obj%FEMDomain%Boundary%NBoundNodID, 1)
            do j = 1, size(obj%FEMDomain%Boundary%NBoundNodID, 2)
               if (obj%FEMDomain%Boundary%NBoundNodID(i, j) > 0) then
                  node = obj%FEMDomain%Boundary%NBoundNodID(i, j)
                  obj%TractionVecGlo(dim_num*(node - 1) + j) = obj%TractionVecGlo(dim_num*(node - 1) + j) + &
                                                               obj%FEMDomain%Boundary%NBoundVal(i, j)
               end if
            end do
         end do
      end if
   end subroutine
!############# Get Traction Vector ######################

!############# Get Internal Vector ######################
   subroutine GetInternalVector(obj)
      class(FiniteDeform_), intent(inout)::obj

      integer(int32) :: i, j, dim_num, node_num, nodeid1, node_num_elmtl, elem_num, dim1

      node_num = size(obj%FEMDomain%Mesh%NodCoord, 1)
      dim_num = size(obj%FEMDomain%Mesh%NodCoord, 2)
      elem_num = size(obj%FEMDomain%Mesh%ElemNod, 1)
      node_num_elmtl = size(obj%FEMDomain%Mesh%ElemNod, 2)
      if (allocated(obj%InternalVecGlo) .and. size(obj%DeformVecGloTot) /= size(obj%InternalVecGlo)) then
         deallocate (obj%InternalVecGlo)
      end if
      if (.not. allocated(obj%InternalVecGlo)) then
         allocate (obj%InternalVecGlo(size(obj%DeformVecGloTot)))
      end if
      obj%InternalVecGlo(:) = 0.0d0
      do i = 1, elem_num
      do j = 1, node_num_elmtl
         nodeid1 = obj%FEMDomain%Mesh%ElemNod(i, j)
         do dim1 = 1, dim_num
            obj%InternalVecGlo(dim_num*(nodeid1 - 1) + dim1) = &
               obj%InternalVecGlo(dim_num*(nodeid1 - 1) + dim1) + &
               (obj%DeformStressRHS(i, (dim_num*(j - 1)) + dim1))
         end do
      end do
      end do

   end subroutine
!############# Get Internal Vector ######################

!############# Get Residual Vector ######################
   subroutine GetResidualVector(obj)
      class(FiniteDeform_), intent(inout)::obj

      integer(int32) :: i, j, dim_num, node

      if (allocated(obj%ResidualVecGlo) .and. size(obj%DeformVecGloTot) /= size(obj%ResidualVecGlo)) then
         deallocate (obj%ResidualVecGlo)
      end if
      if (.not. allocated(obj%ResidualVecGlo)) then
         allocate (obj%ResidualVecGlo(size(obj%DeformVecGloTot)))
      end if

      obj%ResidualVecGlo(:) = obj%TractionVecGlo(:) - obj%InternalVecGlo(:)

   end subroutine
!############# Get Residual Vector ######################

   subroutine UpdateStressMeasure(obj)
      class(FiniteDeform_), intent(inout)::obj
      if (obj%infinitesimal .eqv. .true.) then
         obj%DeformStress(:, :, :) = obj%DeformStress(:, :, :) + obj%DeformStressinc(:, :, :)
         obj%DeformStressinc(:, :, :) = 0.0d0
      end if
   end subroutine

!############# Uodate strain variables ######################
   subroutine UpdateStrainMeasure(obj)
      class(FiniteDeform_), intent(inout)::obj

      integer(int32) :: i, j, dim_num

      dim_num = size(obj%FEMDomain%Mesh%NodCoord, 2)

      obj%VolInitCurrEBE(:, 3) = obj%VolInitCurrEBE(:, 2) - obj%VolInitCurrEBE(:, 1)
      obj%VolInitCurrEBE(:, 1) = obj%VolInitCurrEBE(:, 2)
      do i = 1, size(obj%DeformStrain, 1)

         do j = 1, size(obj%DeformStrain, 2)
            obj%DeformStrain(i, j, 1) = obj%DeformStrain(i, j, 4)
            obj%DeformStrain(i, j, 2) = obj%DeformStrain(i, j, 5)
            obj%DeformStrain(i, j, 3) = obj%DeformStrain(i, j, 6)
            obj%DeformStrain(i, j, 7) = obj%DeformStrain(i, j, 11)
            obj%DeformStrain(i, j, 8) = obj%DeformStrain(i, j, 12)
            obj%DeformStrain(i, j, 9) = obj%DeformStrain(i, j, 13)
            obj%DeformStrain(i, j, 10) = obj%DeformStrain(i, j, 14)

            if (dim_num == 3) then
               obj%DeformStrain(:, :, 16) = obj%DeformStrain(:, :, 19)
               obj%DeformStrain(:, :, 17) = obj%DeformStrain(:, :, 20)
               obj%DeformStrain(:, :, 18) = obj%DeformStrain(:, :, 21)
               obj%DeformStrain(:, :, 22) = obj%DeformStrain(:, :, 27)
               obj%DeformStrain(:, :, 23) = obj%DeformStrain(:, :, 28)
               obj%DeformStrain(:, :, 24) = obj%DeformStrain(:, :, 29)
               obj%DeformStrain(:, :, 25) = obj%DeformStrain(:, :, 30)
               obj%DeformStrain(:, :, 26) = obj%DeformStrain(:, :, 31)

               obj%DeformStrain(:, :, 32) = obj%DeformStrain(:, :, 41)
               obj%DeformStrain(:, :, 33) = obj%DeformStrain(:, :, 42)
               obj%DeformStrain(:, :, 34) = obj%DeformStrain(:, :, 43)
               obj%DeformStrain(:, :, 35) = obj%DeformStrain(:, :, 44)
               obj%DeformStrain(:, :, 36) = obj%DeformStrain(:, :, 45)
               obj%DeformStrain(:, :, 37) = obj%DeformStrain(:, :, 46)
               obj%DeformStrain(:, :, 38) = obj%DeformStrain(:, :, 47)
               obj%DeformStrain(:, :, 39) = obj%DeformStrain(:, :, 48)
               obj%DeformStrain(:, :, 40) = obj%DeformStrain(:, :, 49)

            end if

         end do
      end do

   end subroutine
!############# Uodate strain variables ######################

!############# Reaction Force at Loading Dirichlet Boundary ######################
   subroutine DisplayReactionForce(obj, id)
      class(FiniteDeform_), intent(in)::obj
      integer(int32), optional, intent(in) :: id

      integer(int32) :: i, j, k, dim_num, dbc_num
      character(200) :: filename
      real(real64), allocatable :: ReactionForce(:), ReactionForce_wall(:)
      real(real64) :: val

      dim_num = size(obj%FEMDomain%Boundary%DBoundNodID, 2)
      dbc_num = size(obj%FEMDomain%Boundary%DBoundNodID, 1)
      allocate (ReactionForce(dim_num), ReactionForce_wall(dim_num))
      ReactionForce(:) = 0.0d0
      ReactionForce_wall(:) = 0.0d0
      do i = 1, dim_num
         do j = 1, dbc_num
            k = obj%FEMDomain%Boundary%DBoundNodID(j, i)
            val = obj%FEMDomain%Boundary%DBoundValInc(j, i)
            if (k < 1) then
               cycle
            elseif (k >= 1 .and. Val /= 0.0d0) then
               ReactionForce(i) = ReactionForce(i) + obj%InternalVecGlo(dim_num*(k - 1) + i)
            elseif (k >= 1 .and. Val == 0.0d0) then
               ReactionForce_wall(i) = ReactionForce_wall(i) + abs(obj%InternalVecGlo(dim_num*(k - 1) + i))
            else
               cycle
            end if
         end do
      end do

      k = input(default=0, option=id)
      filename = "ReactionForce"//str(k)//"_wall.txt"
      if (obj%Step == 1) then
         open (101, file="all_"//FileName, status="replace")
      else
         open (101, file="all_"//FileName, position="append")
      end if
      write (101, *) obj%Step, ReactionForce(:)
      close (101)

      if (obj%Step == 1) then
         open (101, file=FileName, status="replace")
      else
         open (101, file=FileName, position="append")
      end if
      write (101, *) obj%Step, ReactionForce_wall(:)
      close (101)

   end subroutine
!############# Reaction Force at Loading Dirichlet Boundary ######################

! ##################################################
   subroutine getDBCVectorDeform(obj, DBCvec)
      class(FiniteDeform_), intent(in)::obj
      real(real64), allocatable, intent(inout)::DBCvec(:, :)
      integer(int32) :: i, j, n, m, k, l
      n = size(obj%FEMDomain%Mesh%NodCoord, 1)
      m = size(obj%FEMDomain%Mesh%NodCoord, 2)
      if (.not. allocated(DBCvec)) then
         allocate (DBCvec(n, m))
         DBCvec(:, :) = 0.0d0
      end if

      ! check number of DBC
      do i = 1, size(obj%FEMDomain%Boundary%DBoundNum)
         k = countif(Array=obj%FEMDomain%Boundary%DBoundNodID(:, i), Value=-1, notEqual=.true.)
         l = obj%FEMDomain%Boundary%DBoundNum(i)
         if (k /= l) then
            print *, "Caution :: FiniteDeformationClass::getDBCVector :: check number of DBC :: k /= l"
         end if
      end do

      do i = 1, size(obj%FEMDomain%Boundary%DBoundNodID, 1)
         do j = 1, size(obj%FEMDomain%Boundary%DBoundNodID, 2)
            if (obj%FEMDomain%Boundary%DBoundNodID(i, j) <= 0) then
               cycle
            end if
            DBCvec(obj%FEMDomain%Boundary%DBoundNodID(i, j), j) = obj%FEMDomain%Boundary%DBoundVal(i, j)
         end do
      end do

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

! ##################################################
   subroutine getDispVectorDeform(obj, Vector)
      class(FiniteDeform_), intent(in)::obj
      real(real64), allocatable, intent(inout)::Vector(:, :)
      integer(int32) :: i, j, n, m

      n = size(obj%FEMDomain%Mesh%NodCoord, 1)
      m = size(obj%FEMDomain%Mesh%NodCoord, 2)
      if (.not. allocated(Vector)) then
         allocate (Vector(n, m))
         Vector(:, :) = 0.0d0
      end if

      do i = 1, n
         do j = 1, m
            Vector(i, j) = obj%DeformVecGloTot(m*(i - 1) + j)
         end do
      end do
   end subroutine
! ##################################################

! ##################################################
   subroutine exportFiniteDeform(obj, restart, path)
      class(FiniteDeform_), intent(inout) :: obj
      logical, optional, intent(in) :: restart
      character(*), intent(in) :: path
      type(IO_) :: f, ff
      integer(int32) :: fh_

      if (present(restart)) then
         ! make new dir
         call execute_command_line("mkdir -p "//path//"/FiniteDeform")

         if (associated(obj%FEMDomain)) then
            call obj%FEMDomain%export(path=path//"/FiniteDeform", restart=restart)
         end if

         call f%open("./", path//"/FiniteDeform", "/FiniteDeform.prop")
         write (f%fh, *) obj%DeformStress(:, :, :)
         write (f%fh, *) obj%DeformStrain(:, :, :)
         write (f%fh, *) obj%DeformStressInit(:, :, :)
         write (f%fh, *) obj%DeformStressinc(:, :, :)
         write (f%fh, *) obj%DeformStressMat(:, :, :)
         write (f%fh, *) obj%DeformStressRHS(:, :)
         write (f%fh, *) obj%DeformVecEBETot(:, :)
         write (f%fh, *) obj%DeformVecEBEInc(:, :)
         write (f%fh, *) obj%DeformVecGloTot(:)
         write (f%fh, *) obj%DeformVecGloInc(:)
         write (f%fh, *) obj%TractionVecGlo(:)
         write (f%fh, *) obj%ResidualVecGlo(:)
         write (f%fh, *) obj%InternalVecGlo(:)
         write (f%fh, *) obj%VolInitCurrEBE(:, :)
         write (f%fh, *) obj%YoungsModulus(:)
         write (f%fh, *) obj%PoissonsRatio(:)
         write (f%fh, *) obj%PorePressure(:)
         write (f%fh, *) obj%dt, obj%error, obj%reactionforce
         write (f%fh, *) obj%nr_tol
         write (f%fh, *) obj%ReducedIntegration
         write (f%fh, *) obj%infinitesimal
         write (f%fh, *) obj%itr
         write (f%fh, *) obj%Step
         call f%close()
      end if

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

! ##################################################
   function getVolumeDeform(obj) result(volume)
      class(FiniteDeform_), intent(inout) :: obj
      real(real64), allocatable :: volume(:)
      integer(int32) :: numnode, numelem, i, j, gp_num

      obj%FEMDomain%ShapeFunction%ElemType = obj%FEMDomain%Mesh%GetElemType()
      call SetShapeFuncType(obj%FEMDomain%ShapeFunction)
      gp_num = obj%FEMDomain%ShapeFunction%NumOfGp
      allocate (volume(size(obj%FEMDomain%Mesh%ElemNod, 1)))
      do i = 1, size(obj%FEMDomain%Mesh%ElemNod, 1)
         do j = 1, obj%FEMDomain%ShapeFunction%NumOfGp !�K�E�X�ϕ����ƃ��[�v
            ! -----J�}�g���N�X�̌v�Z-----------------------------------------
            call GetAllShapeFunc(obj%FEMDomain%ShapeFunction, elem_id=i, nod_coord=obj%FEMDomain%Mesh%NodCoord, &
                                 elem_nod=obj%FEMDomain%Mesh%ElemNod, OptionalGpID=j)
            volume(i) = obj%FEMDomain%ShapeFunction%detJ
         end do
      end do

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

end module FiniteDeformationClass

