#include "../MEF90/mef90.inc"
module MEF90_APPEND(m_MEF90_DefMechPlasticityVonMises,MEF90_DIM)D
#include "petsc/finclude/petsc.h"

use MEF90_APPEND(m_MEF90_DefMechPlasticityCtx,MEF90_DIM)D
implicit none(type, external)
private
public :: FHG_VONMISES
public :: FHG_VONMISESPLANETHEORY
public :: FHG_VONMISES1D

contains
#undef __FUNCT__
#define __FUNCT__ "FHG_VONMISES"
!!! author: Erwan Tanne (2015, erwan.tanne@gmail.com)
!!!
!!!  fhg: VonMises
!!!

subroutine FHG_VONMISES(x, f, h, g, myctx) bind(c)
   use, intrinsic :: iso_c_binding

   real(kind=c_double)                       :: x(*)
   real(kind=c_double)                       :: f(*)
   real(kind=c_double)                       :: h(*)
   real(kind=c_double)                       :: g(*)
   real(Kind=Kr)                             :: StiffnessA, PlasticStrainCumulated
   real(Kind=Kr)                             :: StiffnessB, Stiffness
   type(c_ptr), intent(in), value            :: myctx

   type(MEF90DefMechPlasticityCtx), pointer  :: myctx_ptr
   type(MEF90_MATS)                          :: xMatS
   type(MEF90_MATS)                          :: Stress

   xMatS = x(1:SIZEOFMEF90_MATS)
      !! This is the fortran equivalent of casting ctx into a c_ptr
   call c_f_pointer(myctx, myctx_ptr)

      !! Select which softening young model
   if (myctx_ptr%CoefficientLinSoft == 0) then
      StiffnessA = (1.0_kr - myctx_ptr%Damage)**2 + myctx_ptr%residualStiffness
      StiffnessB = (1.0_kr - myctx_ptr%Damage)**myctx_ptr%DuctileCouplingPower + myctx_ptr%residualStiffness
   else
      StiffnessA = ((1.0_kr - myctx_ptr%Damage)**2 / (1.0_kr + (myctx_ptr%CoefficientLinSoft - 1.0_kr) * (1.0_kr - (1.0_kr - myctx_ptr%Damage)**2))) + myctx_ptr%residualStiffness
      StiffnessB = (1.0_kr - myctx_ptr%Damage)**myctx_ptr%DuctileCouplingPower + myctx_ptr%residualStiffness
   end if

   Stress = myctx_ptr%HookesLaw * (myctx_ptr%totalStrain - xMatS)

   if (myctx_ptr%isLinearIsotropicHardening .eqv. .true.) then
      PlasticStrainCumulated = (myctx_ptr%cumulatedDissipatedPlasticEnergy + (Stress.DotP. (xMatS - myctx_ptr%PlasticStrainOld))) / myctx_ptr%YieldStress
   else
      PlasticStrainCumulated = 0.0_kr
   end if

   Stiffness = StiffnessB * (1.0_kr + PlasticStrainCumulated)

   f(1) = ((myctx_ptr%HookesLaw * (xMatS - myctx_ptr%PlasticStrainOld)) .DotP. (xMatS - myctx_ptr%PlasticStrainOld))
   g(1) = StiffnessA * sqrt(MEF90_DIM / (MEF90_DIM - 1.0_kr) * (deviatoricPart(Stress) .DotP.deviatoricPart(Stress))) - ((1.0_kr - myctx_ptr%residualYieldStress) * Stiffness + myctx_ptr%residualYieldStress) * myctx_ptr%YieldStress
   g(2) = -PlasticStrainCumulated
   h(1) = Trace(xMatS)
end subroutine FHG_VONMISES

#undef __FUNCT__
#define __FUNCT__ "FHG_VONMISESPLANETHEORY"
!!! author: Erwan Tanne (2015, erwan.tanne@gmail.com)
!!!
!!!  fhg: VonMises
!!!

subroutine FHG_VONMISESPLANETHEORY(x, f, h, g, myctx) bind(c)
   use, intrinsic :: iso_c_binding

   real(kind=c_double)                       :: x(*)
   real(kind=c_double)                       :: f(*)
   real(kind=c_double)                       :: h(*)
   real(kind=c_double)                       :: g(*)
   real(Kind=Kr)                           :: StiffnessA         !! Stiffness = a(alpha)/b(alpha)
   real(Kind=Kr)                           :: StiffnessB
   type(c_ptr), intent(in), value              :: myctx

   type(MEF90DefMechPlasticityCtx), pointer   :: myctx_ptr
   type(MEF90_MATS)                          :: xMatS
   real(Kind=Kr)                           :: lambda, mu, E, nu
   type(MatS3D)                              :: Strain, PlasticStrainFlow, Stress

   xMatS = x(1:SIZEOFMEF90_MATS)
      !! This is the fortran equivalent of casting ctx into a c_ptr
   call c_f_pointer(myctx, myctx_ptr)

      !! Select which softening young model
   if (myctx_ptr%CoefficientLinSoft == 0) then
      StiffnessA = (1.0_kr - myctx_ptr%Damage)**2 + myctx_ptr%residualStiffness
      StiffnessB = (1.0_kr - myctx_ptr%Damage)**myctx_ptr%DuctileCouplingPower + myctx_ptr%residualStiffness
   else
      StiffnessA = ((1.0_kr - myctx_ptr%Damage)**2 / (1.0_kr + (myctx_ptr%CoefficientLinSoft - 1.0_kr) * (1.0_kr - (1.0_kr - myctx_ptr%Damage)**2))) + myctx_ptr%residualStiffness
      StiffnessB = (1.0_kr - myctx_ptr%Damage)**myctx_ptr%DuctileCouplingPower + myctx_ptr%residualStiffness
   end if

   E = myctx_ptr%HookesLaw%YoungsModulus
   nu = myctx_ptr%HookesLaw%PoissonRatio
   mu = E / (1.0_kr + nu)*.5_kr
   lambda = E * nu / (1.0_kr + nu) / (1 - 2.0_kr * nu)

   if (myctx_ptr%isPlaneStress .eqv. .true.) then
         !! If plane stress
      Strain = 0.0_kr
      Strain%XX = myctx_ptr%totalStrain%XX
      Strain%YY = myctx_ptr%totalStrain%YY
      Strain%XY = myctx_ptr%totalStrain%XY
      Strain%ZZ = (-lambda * Trace(myctx_ptr%totalStrain) - 2 * mu * Trace(myctx_ptr%plasticStrainPrevious)) / (lambda + 2 * mu)

      PlasticStrainFlow = 0.0_kr
      PlasticStrainFlow%XX = xMatS%XX - myctx_ptr%PlasticStrainOld%XX
      PlasticStrainFlow%YY = xMatS%YY - myctx_ptr%PlasticStrainOld%YY
      PlasticStrainFlow%XY = xMatS%XY - myctx_ptr%PlasticStrainOld%XY
      PlasticStrainFlow%ZZ = -(PlasticStrainFlow%XX + PlasticStrainFlow%YY)

      Stress = 0.0_kr
      Stress%XX = lambda * (Trace(Strain)) + 2 * mu * (Strain%XX - xMatS%XX)
      Stress%YY = lambda * (Trace(Strain)) + 2 * mu * (Strain%YY - xMatS%YY)
      Stress%XY = 2 * mu * (Strain%XY - xMatS%XY)
      Stress%ZZ = lambda * (Trace(Strain)) + 2 * mu * (Strain%ZZ + Trace(xMatS))
   else
         !! If plane strain
      Strain = 0.0_kr
      Strain%XX = myctx_ptr%totalStrain%XX
      Strain%YY = myctx_ptr%totalStrain%YY
      Strain%XY = myctx_ptr%totalStrain%XY

      PlasticStrainFlow = 0.0_kr
      PlasticStrainFlow%XX = xMatS%XX - myctx_ptr%PlasticStrainOld%XX
      PlasticStrainFlow%YY = xMatS%YY - myctx_ptr%PlasticStrainOld%YY
      PlasticStrainFlow%XY = xMatS%XY - myctx_ptr%PlasticStrainOld%XY
      PlasticStrainFlow%ZZ = -(PlasticStrainFlow%XX + PlasticStrainFlow%YY)

      Stress = 0.0_kr
      Stress%XX = lambda * (Trace(Strain)) + 2 * mu * (Strain%XX - xMatS%XX)
      Stress%YY = lambda * (Trace(Strain)) + 2 * mu * (Strain%YY - xMatS%YY)
      Stress%XY = 2 * mu * (Strain%XY - xMatS%XY)
      Stress%ZZ = lambda * (Trace(Strain)) + 2 * mu * (Strain%ZZ + Trace(xMatS))
   end if

   if (myctx_ptr%isNoPlCoupling .eqv. .true.) then
      f(1) = (PlasticStrainFlow.DotP.PlasticStrainFlow)
      g(1) = StiffnessA * sqrt((3.0 / 2.0) * (deviatoricPart(Stress) .dotP.deviatoricPart(Stress))) - myctx_ptr%YieldStress
   else
      f(1) = (PlasticStrainFlow.DotP.PlasticStrainFlow)
      g(1) = StiffnessA * sqrt((3.0 / 2.0) * (deviatoricPart(Stress) .dotP.deviatoricPart(Stress))) - ((1.0_kr - myctx_ptr%residualYieldStress) * StiffnessB + myctx_ptr%residualYieldStress) * myctx_ptr%YieldStress
   end if

end subroutine FHG_VONMISESPLANETHEORY

#undef __FUNCT__
#define __FUNCT__ "FHG_VONMISES1D"
!!! author: Erwan Tanne (2015, erwan.tanne@gmail.com)
!!!
!!!  fhg: VonMises
!!!

subroutine FHG_VONMISES1D(x, f, h, g, myctx) bind(c)
   use, intrinsic :: iso_c_binding

   real(kind=c_double)                       :: x(*)
   real(kind=c_double)                       :: f(*)
   real(kind=c_double)                       :: h(*)
   real(kind=c_double)                       :: g(*)
   real(Kind=Kr)                           :: StiffnessA
   real(Kind=Kr)                           :: StiffnessB
   real(Kind=Kr)                           :: Stiffness, PlasticStrainCumulated
   type(c_ptr), intent(in), value              :: myctx

   type(MEF90DefMechPlasticityCtx), pointer   :: myctx_ptr
   type(MEF90_MATS)                          :: xMatS, Stress

   xMatS = x(1:SIZEOFMEF90_MATS)
      !! This is the fortran equivalent of casting ctx into a c_ptr
   call c_f_pointer(myctx, myctx_ptr)

       !! Select which softening young model
   if (myctx_ptr%CoefficientLinSoft == 0) then
      StiffnessA = (1.0_kr - myctx_ptr%Damage)**2
      StiffnessB = (1.0_kr - myctx_ptr%Damage)**myctx_ptr%DuctileCouplingPower
   else
      StiffnessA = ((1.0_kr - myctx_ptr%Damage)**2 / (1.0_kr + (myctx_ptr%CoefficientLinSoft - 1.0_kr) * (1.0_kr - (1.0_kr - myctx_ptr%Damage)**2)))
      StiffnessB = (1.0_kr - myctx_ptr%Damage)**myctx_ptr%DuctileCouplingPower
   end if

   Stress = myctx_ptr%HookesLaw * (myctx_ptr%totalStrain - xMatS)

   if (myctx_ptr%isLinearIsotropicHardening .eqv. .true.) then
      PlasticStrainCumulated = (myctx_ptr%cumulatedDissipatedPlasticEnergy + (Stress.DotP. (xMatS - myctx_ptr%PlasticStrainOld))) / myctx_ptr%YieldStress
   else
      PlasticStrainCumulated = 0.0_kr
   end if

   Stiffness = StiffnessB * (1.0_kr + PlasticStrainCumulated)

   f(1) = ((myctx_ptr%HookesLaw * (xMatS - myctx_ptr%PlasticStrainOld)) .DotP. (xMatS - myctx_ptr%PlasticStrainOld))
   g(1) = StiffnessA * sqrt(MEF90_DIM / (MEF90_DIM - 1.0_kr) * (deviatoricPart(Stress) .DotP.deviatoricPart(Stress))) - myctx_ptr%YieldStress * Stiffness
   g(2) = -PlasticStrainCumulated
   h(1) = xMatS%YY
   h(2) = xMatS%XY
end subroutine FHG_VONMISES1D

end module MEF90_APPEND(m_MEF90_DefMechPlasticityVonMises,MEF90_DIM)D