module StressClass
   use, intrinsic :: iso_fortran_env
   use MathClass
   use StrainClass
   implicit none

   type :: Stress_

      ! hyper
      real(real64), allocatable :: sigma(:, :)
      real(real64), allocatable :: sigma_n(:, :)
      real(real64), allocatable :: S(:, :)
      real(real64), allocatable :: P(:, :)

      ! derivatives
      real(real64), allocatable :: dSdC(:, :, :, :)

      ! hypo
      real(real64), allocatable :: sigma_dot(:, :)
      real(real64), allocatable :: sigma_j(:, :)
      real(real64), allocatable :: sigma_o(:, :)
      real(real64), allocatable :: sigma_t(:, :)

      ! derivatives (dsigma/deps)
      real(real64), allocatable :: E(:, :, :, :)

      integer(int32) :: TheoryID

      character*40 :: StrainTheory

      ! Please input one of following keyphrases

      ! 1 : Finite_Elasticity
      ! 2 : Finite_ElastoPlasticity
      ! 3 : Infinitesimal_Elasticity
      ! 4 : Infinitesimal_ElastoPlasticity
      ! 5 : Small_strain

   contains
      procedure, public :: init => initStress
      procedure, public :: getRate => getStressRate
      procedure, public :: getStress => getStress
      procedure, public :: getDerivative => getStressDerivative
      procedure, public :: delete => deleteStress
   end type
contains

! ###############################
   subroutine initStress(obj, StrainTheory)
      class(Stress_), intent(inout) :: obj
      character(*), intent(in) :: StrainTheory
      real(real64) :: delta(3, 3)

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

      if (allocated(obj%sigma)) then
         call obj%delete()
      end if

      obj%StrainTheory = StrainTheory

      if (obj%StrainTheory == "Finite_Elasticity") then
         obj%theoryID = 1

         ! hyper
         allocate (obj%sigma(3, 3))
         allocate (obj%sigma_n(3, 3))
         allocate (obj%S(3, 3))
         allocate (obj%P(3, 3))

         ! hypo
         allocate (obj%sigma_dot(0, 0))
         allocate (obj%sigma_j(0, 0))
         allocate (obj%sigma_o(0, 0))
         allocate (obj%sigma_t(0, 0))

         obj%sigma(:, :) = 0.0d0
         obj%sigma_n(:, :) = 0.0d0
         obj%S(:, :) = 0.0d0
         obj%P(:, :) = 0.0d0

      elseif (obj%StrainTheory == "Finite_ElastoPlasticity") then
         obj%theoryID = 2

         ! hyper
         allocate (obj%sigma(3, 3))
         allocate (obj%sigma_n(3, 3))
         allocate (obj%S(3, 3))
         allocate (obj%P(3, 3))

         ! hypo
         allocate (obj%sigma_dot(0, 0))
         allocate (obj%sigma_j(0, 0))
         allocate (obj%sigma_o(0, 0))
         allocate (obj%sigma_t(0, 0))

         obj%sigma(:, :) = 0.0d0
         obj%sigma_n(:, :) = 0.0d0
         obj%S(:, :) = 0.0d0
         obj%P(:, :) = 0.0d0

      elseif (obj%StrainTheory == "Infinitesimal_Elasticity") then
         obj%theoryID = 3

         ! hyper
         allocate (obj%sigma(3, 3))
         allocate (obj%sigma_n(3, 3))
         allocate (obj%S(0, 0))
         allocate (obj%P(0, 0))

         ! hypo
         allocate (obj%sigma_dot(3, 3))
         allocate (obj%sigma_j(3, 3))
         allocate (obj%sigma_o(3, 3))
         allocate (obj%sigma_t(3, 3))

         obj%sigma(:, :) = 0.0d0
         obj%sigma_n(:, :) = 0.0d0
         obj%sigma_dot(:, :) = 0.0d0
         obj%sigma_j(:, :) = 0.0d0
         obj%sigma_o(:, :) = 0.0d0
         obj%sigma_t(:, :) = 0.0d0

      elseif (obj%StrainTheory == "Infinitesimal_ElastoPlasticity") then
         obj%theoryID = 4

         ! hyper
         allocate (obj%sigma(3, 3))
         allocate (obj%sigma_n(3, 3))
         allocate (obj%S(0, 0))
         allocate (obj%P(0, 0))

         ! hypo
         allocate (obj%sigma_dot(3, 3))
         allocate (obj%sigma_j(3, 3))
         allocate (obj%sigma_o(3, 3))
         allocate (obj%sigma_t(3, 3))

         obj%sigma(:, :) = 0.0d0
         obj%sigma_n(:, :) = 0.0d0
         obj%sigma_dot(:, :) = 0.0d0
         obj%sigma_j(:, :) = 0.0d0
         obj%sigma_o(:, :) = 0.0d0
         obj%sigma_t(:, :) = 0.0d0

      elseif (obj%StrainTheory == "Small_strain") then
         obj%theoryID = 5

         ! hyper
         allocate (obj%sigma(3, 3))
         allocate (obj%sigma_n(3, 3))
         allocate (obj%S(0, 0))
         allocate (obj%P(0, 0))

         ! hypo
         allocate (obj%sigma_dot(0, 0))
         allocate (obj%sigma_j(0, 0))
         allocate (obj%sigma_o(0, 0))
         allocate (obj%sigma_t(0, 0))

         obj%sigma(:, :) = 0.0d0
         obj%sigma_n(:, :) = 0.0d0

      else
         print *, StrainTheory
         print *, "Please input one of following keyphrases"
         print *, " ---->"
         print *, "Finite_Elasticity"
         print *, "Finite_ElastoPlasticity"
         print *, "Infinitesimal_Elasticity"
         print *, "Infinitesimal_ElastoPlasticity"
         print *, "Small_strain"
         print *, " <----"
      end if

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

! ###############################
   subroutine deleteStress(obj)
      class(Stress_), intent(inout) :: obj

      ! hyper
      deallocate (obj%sigma)
      deallocate (obj%sigma_n)
      deallocate (obj%S)
      deallocate (obj%P)

      ! hypo
      deallocate (obj%sigma_dot)
      deallocate (obj%sigma_j)
      deallocate (obj%sigma_o)
      deallocate (obj%sigma_t)

      obj%TheoryID = 0

      obj%StrainTheory = " "

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

! ###############################
   subroutine getStressRate(obj, Strain, Type)
      class(Stress_), intent(inout) :: obj
      class(Strain_), intent(inout) :: strain
      character(*), intent(in) :: Type

      if (Type == "Jaumann") then
         obj%sigma_j = obj%sigma_dot + matmul(obj%sigma, strain%w) - matmul(strain%w, obj%sigma)
      elseif (Type == "Oldroyd") then
         obj%sigma_o = obj%sigma_dot + matmul(strain%l, obj%sigma) + matmul(obj%sigma, transpose(strain%l))
      elseif (Type == "Truesdell") then
         obj%sigma_t = obj%sigma_dot - matmul(strain%l, obj%sigma) &
                       - matmul(obj%sigma, transpose(strain%l)) + trace(strain%l)*obj%sigma
      else
         print *, "ERROR :: getStressRate :: invalid stress ", Type
         return
      end if

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

! ###############################
   subroutine getStress(obj, Strain, ConstitutiveModel, lambda, mu, K, G, c, phi, TimeIntegral, StressRate, dt)
      class(Stress_), intent(inout) :: obj
      class(Strain_), intent(inout) :: strain
      character(*), optional, intent(in) :: ConstitutiveModel, TimeIntegral, StressRate
      real(real64), optional, intent(in) :: lambda, mu, K, G, c, phi, dt
      real(real64), allocatable :: F_inv(:, :), Cp_inv(:, :), delta(:, :), C_inv(:, :), M(:, :), E(:, :, :, :)
      real(real64) :: detC, detCp, f, J2_M, I1_M, theta_M, BI, fc, f_MC
      integer(int32) :: i, j, l

      allocate (delta(3, 3))
      delta(:, :) = 0.0d0
      delta(1, 1) = 1.0d0
      delta(2, 2) = 1.0d0
      delta(3, 3) = 1.0d0

      if (size(Strain%F, 1) == 3) then
         ! Finite Strain Theory
         allocate (F_inv(3, 3), Cp_inv(3, 3), C_inv(3, 3), M(3, 3))
         call inverse_rank_2(strain%F, F_inv)
         call inverse_rank_2(strain%Cp, Cp_inv)
         call inverse_rank_2(strain%C, C_inv)
         detC = strain%detF*strain%detF
         detCp = det_mat(strain%Cp, size(strain%Cp, 1))

         if (ConstitutiveModel == "StVenant") then
            ! Conventional St.Venant
            obj%S(:, :) = lambda*0.50d0*trace(strain%C - delta)*delta(:, :) + mu*(strain%C(:, :) - delta(:, :))
            obj%sigma(:, :) = 1.0d0/strain%detF*(matmul(strain%F, matmul(obj%S, transpose(strain%F))))
         elseif (ConstitutiveModel == "StVenant_Kirchhoff") then
            ! St.Venant-Kirchhoff (Wallin and Ristinmaa, 2005 )
            print *, "Finite Strain St.Venant-Kirchhoff will be implemented soon."
         elseif (ConstitutiveModel == "NeoHookean") then
            ! Conventional Neo-Hookean
            obj%S(:, :) = lambda*log(strain%detF)*C_inv(:, :) + mu*(delta(:, :) - C_inv(:, :))
            obj%sigma(:, :) = 1.0d0/strain%detF*(matmul(strain%F, matmul(obj%S, transpose(strain%F))))
         elseif (ConstitutiveModel == "MCDP") then
            ! elastic part is Modified Neo-Hookean (Vladimirov, 2008; 2010)
            obj%S(:, :) = lambda*0.50d0*(detC/detCp - 1.0d0)*C_inv(:, :) &
                          - mu*(Cp_inv(:, :) - C_inv(:, :))
            obj%sigma(:, :) = 1.0d0/strain%detF*(matmul(strain%F, matmul(obj%S, transpose(strain%F))))

            ! Check yield criterion (Mohr-Coulomb)
            M(:, :) = matmul(strain%C, obj%S)
            J2_M = invariant_J2(M)
            I1_M = invariant_I1(M)
            theta_M = invariant_theta(M)
            BI = (1.0d0 + sin(phi))/(1.0d0 - sin(phi))
            fc = (2.0d0*c*cos(phi))/(1.0d0 - sin(phi))
            f_MC = sqrt(J2_M) + ((BI - 1)*I1_M - 3.0d0*fc)/(3.0d0*(BI + 1.0d0)*cos(theta_M) + sqrt(BI - 1.0d0)*sin(theta_M)*fc)

            ! Return-mapping

         elseif (ConstitutiveModel == "CamClay") then
            print *, "FeFp Cam-clay will be implemented soon."
            !obj%S(:,:)=P0*( 1.0d0/ )
            !obj%sigma(:,:) = 1.0d0/obj%detF*( matmul(obj%F, matmul(obj%S, transpose(obj%F)  ) ) )
         else
            print *, "ERROR :: getStressFinitestrain :: invalid stress ", ConstitutiveModel
            return
         end if
      else
         if (size(obj%sigma_t, 1) == 3) then
            ! Infinitesimal strain theory
            if (ConstitutiveModel == "LinearElastic") then
               ! get stiffness matrix
               allocate (E(3, 3, 3, 3))
               E(:, :, :, :) = 0.0d0
               !do i=1,3
               !    do j=1,3
               !        do k=1,3
               !            do l=1,3
               !               E(i,j,k,l)=lambda*delta(i,j)
               !            enddo
               !        enddo
               !    enddo
               !enddo

               if (TimeIntegral == "ForwardEuler") then
                  if (StressRate == "Jaumann") then
                     obj%sigma_j(:, :) = lambda*trace(strain%d) + delta(:, :) + 2.0d0*mu*strain%d(:, :)
                     obj%sigma_dot = obj%sigma_j - matmul(obj%sigma, strain%w) + matmul(strain%w, obj%sigma)
                  elseif (StressRate == "Oldroyd") then
                     obj%sigma_o(:, :) = lambda*trace(strain%d) + delta(:, :) + 2.0d0*mu*strain%d(:, :)
                     obj%sigma_dot = obj%sigma_o - matmul(obj%sigma, strain%w) + matmul(strain%w, obj%sigma)
                  elseif (StressRate == "Truesdell") then
                     obj%sigma_t(:, :) = lambda*trace(strain%d) + delta(:, :) + 2.0d0*mu*strain%d(:, :)
                     obj%sigma_dot = obj%sigma_t - matmul(obj%sigma, strain%w) + matmul(strain%w, obj%sigma)
                  else
                     print *, "ERROR :: getStress :: invalid stress ", ConstitutiveModel
                     return
                  end if
                  obj%sigma(:, :) = obj%sigma_n(:, :) + dt*obj%sigma_dot(:, :)

               elseif (TimeIntegral == "BackwardEuler") then
                  if (StressRate == "Jaumann") then
                     !    obj%sigma_j(:,:)= lambda*trace(strain%d) + delta(:,:) + 2.0d0*mu*strain%d(:,:)
                     !    obj%sigma_dot = obj%sigma_j - matmul(obj%sigma, strain%w) + matmul(strain%w, obj%sigma)
                  elseif (StressRate == "Oldroyd") then
                     !    obj%sigma_o(:,:)= lambda*trace(strain%d) + delta(:,:) + 2.0d0*mu*strain%d(:,:)
                     !    obj%sigma_dot = obj%sigma_o - matmul(obj%sigma, strain%w) + matmul(strain%w, obj%sigma)
                  elseif (StressRate == "Truesdell") then
                     !    obj%sigma_t(:,:)= lambda*trace(strain%d) + delta(:,:) + 2.0d0*mu*strain%d(:,:)
                     !    obj%sigma_dot = obj%sigma_t - matmul(obj%sigma, strain%w) + matmul(strain%w, obj%sigma)
                  else
                     print *, "ERROR :: getStress :: invalid stress ", ConstitutiveModel
                     return
                  end if
                  !obj%sigma(:,:) = obj%sigma_n(:,:) + dt*  obj%sigma_dot(:,:)

               elseif (TimeIntegral == "ClankNicolson") then
                  if (StressRate == "Jaumann") then
                     !    obj%sigma_j(:,:)= lambda*trace(strain%d) + delta(:,:) + 2.0d0*mu*strain%d(:,:)
                     !    obj%sigma_dot = obj%sigma_j - matmul(obj%sigma, strain%w) + matmul(strain%w, obj%sigma)
                  elseif (StressRate == "Oldroyd") then
                     !    obj%sigma_o(:,:)= lambda*trace(strain%d) + delta(:,:) + 2.0d0*mu*strain%d(:,:)
                     !    obj%sigma_dot = obj%sigma_o - matmul(obj%sigma, strain%w) + matmul(strain%w, obj%sigma)
                  elseif (StressRate == "Truesdell") then
                     !    obj%sigma_t(:,:)= lambda*trace(strain%d) + delta(:,:) + 2.0d0*mu*strain%d(:,:)
                     !    obj%sigma_dot = obj%sigma_t - matmul(obj%sigma, strain%w) + matmul(strain%w, obj%sigma)
                  else
                     print *, "ERROR :: getStress :: invalid stress ", ConstitutiveModel
                     return
                  end if
                  !obj%sigma(:,:) = obj%sigma_n(:,:) + dt*  obj%sigma_dot(:,:)

               elseif (TimeIntegral == "SpaceTime") then

                  if (StressRate == "Jaumann") then
                     !    obj%sigma_j(:,:)= lambda*trace(strain%d) + delta(:,:) + 2.0d0*mu*strain%d(:,:)
                     !    obj%sigma_dot = obj%sigma_j - matmul(obj%sigma, strain%w) + matmul(strain%w, obj%sigma)
                  elseif (StressRate == "Oldroyd") then
                     !    obj%sigma_o(:,:)= lambda*trace(strain%d) + delta(:,:) + 2.0d0*mu*strain%d(:,:)
                     !    obj%sigma_dot = obj%sigma_o - matmul(obj%sigma, strain%w) + matmul(strain%w, obj%sigma)
                  elseif (StressRate == "Truesdell") then
                     !    obj%sigma_t(:,:)= lambda*trace(strain%d) + delta(:,:) + 2.0d0*mu*strain%d(:,:)
                     !    obj%sigma_dot = obj%sigma_t - matmul(obj%sigma, strain%w) + matmul(strain%w, obj%sigma)
                  else
                     print *, "ERROR :: getStress :: invalid stress ", ConstitutiveModel
                     return
                  end if
                  !obj%sigma(:,:) = obj%sigma_n(:,:) + dt*  obj%sigma_dot(:,:)

               else
                  print *, "ERROR :: stressClass :: Linear Elastic TimeIntegral = ", TimeIntegral, "is not implemented yet."
               end if
            elseif (ConstitutiveModel == "MCDP") then
               ! get stiffness matrix
               allocate (E(3, 3, 3, 3))
               E(:, :, :, :) = 0.0d0
               !do i=1,3
               !    do j=1,3
               !        do k=1,3
               !            do l=1,3
               !               E(i,j,k,l)=lambda*delta(i,j)
               !            enddo
               !        enddo
               !    enddo
               !enddo

               if (TimeIntegral == "ForwardEuler") then
                  ! Elastic part is linear elastic
                  if (StressRate == "Jaumann") then
                     obj%sigma_j(:, :) = lambda*trace(strain%de) + delta(:, :) + 2.0d0*mu*strain%de(:, :)
                     obj%sigma_dot = obj%sigma_j - matmul(obj%sigma, strain%w) + matmul(strain%w, obj%sigma)
                  elseif (StressRate == "Oldroyd") then
                     obj%sigma_o(:, :) = lambda*trace(strain%de) + delta(:, :) + 2.0d0*mu*strain%de(:, :)
                     obj%sigma_dot = obj%sigma_o - matmul(obj%sigma, strain%w) + matmul(strain%w, obj%sigma)
                  elseif (StressRate == "Truesdell") then
                     obj%sigma_t(:, :) = lambda*trace(strain%de) + delta(:, :) + 2.0d0*mu*strain%de(:, :)
                     obj%sigma_dot = obj%sigma_t - matmul(obj%sigma, strain%w) + matmul(strain%w, obj%sigma)
                  else
                     print *, "ERROR :: getStress :: invalid stress ", ConstitutiveModel
                     return
                  end if
                  obj%sigma(:, :) = obj%sigma_n(:, :) + dt*obj%sigma_dot(:, :)

                  ! Return-mapping

               elseif (TimeIntegral == "BackwardEuler") then
                  if (StressRate == "Jaumann") then
                     !    obj%sigma_j(:,:)= lambda*trace(strain%d) + delta(:,:) + 2.0d0*mu*strain%d(:,:)
                     !    obj%sigma_dot = obj%sigma_j - matmul(obj%sigma, strain%w) + matmul(strain%w, obj%sigma)
                  elseif (StressRate == "Oldroyd") then
                     !    obj%sigma_o(:,:)= lambda*trace(strain%d) + delta(:,:) + 2.0d0*mu*strain%d(:,:)
                     !    obj%sigma_dot = obj%sigma_o - matmul(obj%sigma, strain%w) + matmul(strain%w, obj%sigma)
                  elseif (StressRate == "Truesdell") then
                     !    obj%sigma_t(:,:)= lambda*trace(strain%d) + delta(:,:) + 2.0d0*mu*strain%d(:,:)
                     !    obj%sigma_dot = obj%sigma_t - matmul(obj%sigma, strain%w) + matmul(strain%w, obj%sigma)
                  else
                     print *, "ERROR :: getStress :: invalid stress ", ConstitutiveModel
                     return
                  end if
                  !obj%sigma(:,:) = obj%sigma_n(:,:) + dt*  obj%sigma_dot(:,:)

               elseif (TimeIntegral == "ClankNicolson") then
                  if (StressRate == "Jaumann") then
                     !    obj%sigma_j(:,:)= lambda*trace(strain%d) + delta(:,:) + 2.0d0*mu*strain%d(:,:)
                     !    obj%sigma_dot = obj%sigma_j - matmul(obj%sigma, strain%w) + matmul(strain%w, obj%sigma)
                  elseif (StressRate == "Oldroyd") then
                     !    obj%sigma_o(:,:)= lambda*trace(strain%d) + delta(:,:) + 2.0d0*mu*strain%d(:,:)
                     !    obj%sigma_dot = obj%sigma_o - matmul(obj%sigma, strain%w) + matmul(strain%w, obj%sigma)
                  elseif (StressRate == "Truesdell") then
                     !    obj%sigma_t(:,:)= lambda*trace(strain%d) + delta(:,:) + 2.0d0*mu*strain%d(:,:)
                     !    obj%sigma_dot = obj%sigma_t - matmul(obj%sigma, strain%w) + matmul(strain%w, obj%sigma)
                  else
                     print *, "ERROR :: getStress :: invalid stress ", ConstitutiveModel
                     return
                  end if
                  !obj%sigma(:,:) = obj%sigma_n(:,:) + dt*  obj%sigma_dot(:,:)

               elseif (TimeIntegral == "SpaceTime") then

                  if (StressRate == "Jaumann") then
                     !    obj%sigma_j(:,:)= lambda*trace(strain%d) + delta(:,:) + 2.0d0*mu*strain%d(:,:)
                     !    obj%sigma_dot = obj%sigma_j - matmul(obj%sigma, strain%w) + matmul(strain%w, obj%sigma)
                  elseif (StressRate == "Oldroyd") then
                     !    obj%sigma_o(:,:)= lambda*trace(strain%d) + delta(:,:) + 2.0d0*mu*strain%d(:,:)
                     !    obj%sigma_dot = obj%sigma_o - matmul(obj%sigma, strain%w) + matmul(strain%w, obj%sigma)
                  elseif (StressRate == "Truesdell") then
                     !    obj%sigma_t(:,:)= lambda*trace(strain%d) + delta(:,:) + 2.0d0*mu*strain%d(:,:)
                     !    obj%sigma_dot = obj%sigma_t - matmul(obj%sigma, strain%w) + matmul(strain%w, obj%sigma)
                  else
                     print *, "ERROR :: getStress :: invalid stress ", ConstitutiveModel
                     return
                  end if
                  !obj%sigma(:,:) = obj%sigma_n(:,:) + dt*  obj%sigma_dot(:,:)

               else
                  print *, "ERROR :: stressClass :: MCDP TimeIntegral = ", TimeIntegral, "is not implemented yet."
               end if
            elseif (ConstitutiveModel == "CamClay") then
               ! get stiffness matrix
               allocate (E(3, 3, 3, 3))
               E(:, :, :, :) = 0.0d0
               !do i=1,3
               !    do j=1,3
               !        do k=1,3
               !            do l=1,3
               !               E(i,j,k,l)=lambda*delta(i,j)
               !            enddo
               !        enddo
               !    enddo
               !enddo

               if (TimeIntegral == "ForwardEuler") then
                  ! Elastic part is linear elastic
                  if (StressRate == "Jaumann") then
                     !obj%sigma_j(:,:)= lambda*trace(strain%de) + delta(:,:) + 2.0d0*mu*strain%de(:,:)
                     !obj%sigma_dot = obj%sigma_j - matmul(obj%sigma, strain%w) + matmul(strain%w, obj%sigma)
                  elseif (StressRate == "Oldroyd") then
                     !obj%sigma_o(:,:)= lambda*trace(strain%de) + delta(:,:) + 2.0d0*mu*strain%de(:,:)
                     !obj%sigma_dot = obj%sigma_o - matmul(obj%sigma, strain%w) + matmul(strain%w, obj%sigma)
                  elseif (StressRate == "Truesdell") then
                     !obj%sigma_t(:,:)= lambda*trace(strain%de) + delta(:,:) + 2.0d0*mu*strain%de(:,:)
                     !obj%sigma_dot = obj%sigma_t - matmul(obj%sigma, strain%w) + matmul(strain%w, obj%sigma)
                  else
                     print *, "ERROR :: getStress :: invalid stress ", ConstitutiveModel
                     return
                  end if
                  !obj%sigma(:,:) = obj%sigma_n(:,:) + dt*  obj%sigma_dot(:,:)

                  ! Return-mapping

               elseif (TimeIntegral == "BackwardEuler") then
                  if (StressRate == "Jaumann") then
                     !    obj%sigma_j(:,:)= lambda*trace(strain%d) + delta(:,:) + 2.0d0*mu*strain%d(:,:)
                     !    obj%sigma_dot = obj%sigma_j - matmul(obj%sigma, strain%w) + matmul(strain%w, obj%sigma)
                  elseif (StressRate == "Oldroyd") then
                     !    obj%sigma_o(:,:)= lambda*trace(strain%d) + delta(:,:) + 2.0d0*mu*strain%d(:,:)
                     !    obj%sigma_dot = obj%sigma_o - matmul(obj%sigma, strain%w) + matmul(strain%w, obj%sigma)
                  elseif (StressRate == "Truesdell") then
                     !    obj%sigma_t(:,:)= lambda*trace(strain%d) + delta(:,:) + 2.0d0*mu*strain%d(:,:)
                     !    obj%sigma_dot = obj%sigma_t - matmul(obj%sigma, strain%w) + matmul(strain%w, obj%sigma)
                  else
                     print *, "ERROR :: getStress :: invalid stress ", ConstitutiveModel
                     return
                  end if
                  !obj%sigma(:,:) = obj%sigma_n(:,:) + dt*  obj%sigma_dot(:,:)

               elseif (TimeIntegral == "ClankNicolson") then
                  if (StressRate == "Jaumann") then
                     !    obj%sigma_j(:,:)= lambda*trace(strain%d) + delta(:,:) + 2.0d0*mu*strain%d(:,:)
                     !    obj%sigma_dot = obj%sigma_j - matmul(obj%sigma, strain%w) + matmul(strain%w, obj%sigma)
                  elseif (StressRate == "Oldroyd") then
                     !    obj%sigma_o(:,:)= lambda*trace(strain%d) + delta(:,:) + 2.0d0*mu*strain%d(:,:)
                     !    obj%sigma_dot = obj%sigma_o - matmul(obj%sigma, strain%w) + matmul(strain%w, obj%sigma)
                  elseif (StressRate == "Truesdell") then
                     !    obj%sigma_t(:,:)= lambda*trace(strain%d) + delta(:,:) + 2.0d0*mu*strain%d(:,:)
                     !    obj%sigma_dot = obj%sigma_t - matmul(obj%sigma, strain%w) + matmul(strain%w, obj%sigma)
                  else
                     print *, "ERROR :: getStress :: invalid stress ", ConstitutiveModel
                     return
                  end if
                  !obj%sigma(:,:) = obj%sigma_n(:,:) + dt*  obj%sigma_dot(:,:)

               elseif (TimeIntegral == "SpaceTime") then

                  if (StressRate == "Jaumann") then
                     !    obj%sigma_j(:,:)= lambda*trace(strain%d) + delta(:,:) + 2.0d0*mu*strain%d(:,:)
                     !    obj%sigma_dot = obj%sigma_j - matmul(obj%sigma, strain%w) + matmul(strain%w, obj%sigma)
                  elseif (StressRate == "Oldroyd") then
                     !    obj%sigma_o(:,:)= lambda*trace(strain%d) + delta(:,:) + 2.0d0*mu*strain%d(:,:)
                     !    obj%sigma_dot = obj%sigma_o - matmul(obj%sigma, strain%w) + matmul(strain%w, obj%sigma)
                  elseif (StressRate == "Truesdell") then
                     !    obj%sigma_t(:,:)= lambda*trace(strain%d) + delta(:,:) + 2.0d0*mu*strain%d(:,:)
                     !    obj%sigma_dot = obj%sigma_t - matmul(obj%sigma, strain%w) + matmul(strain%w, obj%sigma)
                  else
                     print *, "ERROR :: getStress :: invalid stress ", ConstitutiveModel
                     return
                  end if
                  !obj%sigma(:,:) = obj%sigma_n(:,:) + dt*  obj%sigma_dot(:,:)

               else
                  print *, "ERROR :: stressClass :: MCDP TimeIntegral = ", TimeIntegral, "is not implemented yet."
               end if
            else
               print *, "ERROR :: getStressInfinitesimal :: invalid stress ", ConstitutiveModel
               return
            end if
         else
            ! Small strain
            if (ConstitutiveModel == "LinearElastic") then
               obj%sigma(:, :) = lambda*trace(strain%eps)*delta(:, :) + 2*mu*strain%eps(:, :)
            elseif (ConstitutiveModel == "MCDP") then
               obj%sigma(:, :) = lambda*trace(strain%eps)*delta(:, :) + 2*mu*strain%eps(:, :)
               ! return-mapping
            elseif (ConstitutiveModel == "CamClay") then
               obj%sigma(:, :) = lambda*trace(strain%eps)*delta(:, :) + 2*mu*strain%eps(:, :)
               ! return-mapping
            else
               print *, "ERROR :: getStressSmallStrain :: invalid stress ", ConstitutiveModel
               return
            end if
         end if
      end if

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

! ###############################
   subroutine getStressDerivative(obj, Strain, ConstitutiveModel, lambda, mu)
      class(Stress_), intent(inout) :: obj
      class(Strain_), intent(inout) :: strain
      character(*), intent(in) :: ConstitutiveModel
      real(real64), optional, intent(in) :: lambda, mu
      real(real64), allocatable :: F_inv(:, :)

      if (size(Strain%F, 1) == 3) then
         ! Finite Strain Theory
         allocate (F_inv(3, 3))
         call inverse_rank_2(strain%F, F_inv)
         if (ConstitutiveModel == "StVenant") then

         elseif (ConstitutiveModel == "NeoHookean") then

         elseif (ConstitutiveModel == "MCDP") then

         elseif (ConstitutiveModel == "CamClay") then

         else
            print *, "ERROR :: getStressFinitestrain :: invalid stress ", ConstitutiveModel
            return
         end if
      else
         if (size(obj%sigma_t, 1) == 3) then
            ! Infinitesimal strain theory
            if (ConstitutiveModel == "LinearElastic") then

            elseif (ConstitutiveModel == "MCDP") then

            elseif (ConstitutiveModel == "CamClay") then

            else
               print *, "ERROR :: getStressInfinitesimal :: invalid stress ", ConstitutiveModel
               return
            end if
         else
            ! Small strain
            if (ConstitutiveModel == "LinearElastic") then

            elseif (ConstitutiveModel == "MCDP") then

            elseif (ConstitutiveModel == "CamClay") then

            else
               print *, "ERROR :: getStressSmallStrain :: invalid stress ", ConstitutiveModel
               return
            end if
         end if
      end if

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

end module
