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