StressClass.f90 Source File


Contents

Source Code


Source Code

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()
    endif

    obj%StrainTheory=StrainTheory

    if(    trim(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(trim(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(trim(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(trim(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(trim(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 *, trim(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 *, " <----"
    endif


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(trim(Type) == "Jaumann" )then
        obj%sigma_j = obj%sigma_dot + matmul(obj%sigma, strain%w) - matmul(strain%w, obj%sigma) 
    elseif(trim(Type) == "Oldroyd" )then
        obj%sigma_o = obj%sigma_dot + matmul(strain%l, obj%sigma) + matmul(obj%sigma, transpose(strain%l) )  
    elseif(trim(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 ",trim(Type)
        return
    endif

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(trim(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(trim(ConstitutiveModel) == "StVenant_Kirchhoff" )then
            ! St.Venant-Kirchhoff (Wallin and Ristinmaa, 2005 )
            print *, "Finite Strain St.Venant-Kirchhoff will be implemented soon."
        elseif(trim(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(trim(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(trim(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 ",trim(ConstitutiveModel)
            return
        endif
    else
        if(size(obj%sigma_t,1)==3 )then
            ! Infinitesimal strain theory
            if(trim(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 ",trim(ConstitutiveModel)
                        return
                    endif
                    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 ",trim(ConstitutiveModel)
                        return
                    endif
                    !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 ",trim(ConstitutiveModel)
                        return
                    endif
                    !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 ",trim(ConstitutiveModel)
                        return
                    endif
                    !obj%sigma(:,:) = obj%sigma_n(:,:) + dt*  obj%sigma_dot(:,:)

                else
                    print *, "ERROR :: stressClass :: Linear Elastic TimeIntegral = ",TimeIntegral,"is not implemented yet."
                endif
            elseif(trim(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 ",trim(ConstitutiveModel)
                        return
                    endif
                    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 ",trim(ConstitutiveModel)
                        return
                    endif
                    !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 ",trim(ConstitutiveModel)
                        return
                    endif
                    !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 ",trim(ConstitutiveModel)
                        return
                    endif
                    !obj%sigma(:,:) = obj%sigma_n(:,:) + dt*  obj%sigma_dot(:,:)

                else
                    print *, "ERROR :: stressClass :: MCDP TimeIntegral = ",TimeIntegral,"is not implemented yet."
                endif
            elseif(trim(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 ",trim(ConstitutiveModel)
                        return
                    endif
                    !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 ",trim(ConstitutiveModel)
                        return
                    endif
                    !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 ",trim(ConstitutiveModel)
                        return
                    endif
                    !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 ",trim(ConstitutiveModel)
                        return
                    endif
                    !obj%sigma(:,:) = obj%sigma_n(:,:) + dt*  obj%sigma_dot(:,:)
                
                else
                    print *, "ERROR :: stressClass :: MCDP TimeIntegral = ",TimeIntegral,"is not implemented yet."
                endif
            else
                print *, "ERROR :: getStressInfinitesimal :: invalid stress ",trim(ConstitutiveModel)
                return
            endif
        else
            ! Small strain
            if(trim(ConstitutiveModel) == "LinearElastic" )then
                obj%sigma(:,:) = lambda*trace(strain%eps)*delta(:,:) + 2*mu*strain%eps(:,:)
            elseif(trim(ConstitutiveModel) == "MCDP" )then
                obj%sigma(:,:) = lambda*trace(strain%eps)*delta(:,:) + 2*mu*strain%eps(:,:)
                ! return-mapping
            elseif(trim(ConstitutiveModel) == "CamClay" )then
                obj%sigma(:,:) = lambda*trace(strain%eps)*delta(:,:) + 2*mu*strain%eps(:,:)
                ! return-mapping
            else
                print *, "ERROR :: getStressSmallStrain :: invalid stress ",trim(ConstitutiveModel)
                return
            endif
        endif
    endif

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(trim(ConstitutiveModel) == "StVenant" )then

        elseif(trim(ConstitutiveModel) == "NeoHookean" )then
        
        elseif(trim(ConstitutiveModel) == "MCDP" )then
        
        elseif(trim(ConstitutiveModel) == "CamClay" )then

        else
            print *, "ERROR :: getStressFinitestrain :: invalid stress ",trim(ConstitutiveModel)
            return
        endif
    else
        if(size(obj%sigma_t,1)==3 )then
            ! Infinitesimal strain theory
            if(trim(ConstitutiveModel) == "LinearElastic" )then
            
            elseif(trim(ConstitutiveModel) == "MCDP" )then
            
            elseif(trim(ConstitutiveModel) == "CamClay" )then
    
            else
                print *, "ERROR :: getStressInfinitesimal :: invalid stress ",trim(ConstitutiveModel)
                return
            endif
        else
            ! Small strain
            if(trim(ConstitutiveModel) == "LinearElastic" )then
            
            elseif(trim(ConstitutiveModel) == "MCDP" )then
            
            elseif(trim(ConstitutiveModel) == "CamClay" )then
    
            else
                print *, "ERROR :: getStressSmallStrain :: invalid stress ",trim(ConstitutiveModel)
                return
            endif
        endif
    endif

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

end module