#include "../MEF90/mef90.inc"
module m_MEF90_DefMechSplit_class
#include "petsc/finclude/petsc.h"

use m_MEF90_LinAlg_class
use m_MEF90_Parameters
use m_MEF90_BaseClass
use m_MEF90_HookesLaw_Class
implicit none(type, external)
private
public :: MEF90DefMechSplit
public :: MEF90DefMechSplit_SmoothPositiveSquare
public :: MEF90DefMechSplit_DSmoothPositiveSquare
public :: MEF90DefMechSplit_D2SmoothPositiveSquare
!!! author: Blaise Bourdin (2020, bourdin@lsu.edu)
!!! author: Blaise Bourdin (2026, bourdin@mcmaster.ca)
!!!
!!!  MEF90DefMechSplit: The abstract class used to define an energy split for
!!!                           handling unilateral contact
!!!

type, extends(MEF90Object), abstract :: MEF90DefMechSplit
   character(len=MEF90MXSTRLEN)                          :: type = ''
   integer                                               :: quadratureOrder = 0
   PetscBool                                             :: isHybrid = PETSC_FALSE
   class(mef90Mat), allocatable                          :: strain
contains
   procedure(setupInterface), pass(self), deferred       :: setup
   procedure(EEDInterface), pass(self), deferred         :: EED
   procedure(DEEDInterface), pass(self), deferred        :: DEED
   procedure(D2EEDInterface), pass(self), deferred       :: D2EED
end type MEF90DefMechSplit

abstract interface

   subroutine setupInterface(self, Strain, ierr)
      use m_MEF90
      import :: MEF90DefMechSplit
      implicit none(type, external)

      class(MEF90DefMechSplit), intent(inout) :: self
      class(mef90Mat), intent(IN)             :: Strain
      PetscErrorCode, intent(inout)           :: ierr
   end subroutine setupInterface

   subroutine EEDInterface(self, HookesLaw, phi, EEDPlus, EEDMinus, ierr)
      use m_MEF90
      import :: MEF90DefMechSplit
      implicit none(type, external)

      class(MEF90DefMechSplit), intent(IN) :: self
      class(MEF90HookesLaw), intent(IN)    :: HookesLaw
      class(mef90Mat), intent(IN)          :: phi
      PetscReal, intent(OUT)               :: EEDPlus, EEDMinus
      PetscErrorCode, intent(inout)        :: ierr
   end subroutine EEDInterface

   subroutine DEEDInterface(self, HookesLaw, phi, DEEDPlus, DEEDMinus, ierr)
      use m_MEF90
      import :: MEF90DefMechSplit
      implicit none(type, external)

      class(MEF90DefMechSplit), intent(IN) :: self
      class(MEF90HookesLaw), intent(IN)    :: HookesLaw
      class(mef90Mat), intent(IN)          :: phi
      PetscReal, intent(OUT)               :: DEEDPlus, DEEDMinus
      PetscErrorCode, intent(inout)        :: ierr
   end subroutine DEEDInterface

   subroutine D2EEDInterface(self, HookesLaw, phi, psi, D2EEDPlus, D2EEDMinus, ierr)
      use m_MEF90
      import :: MEF90DefMechSplit
      implicit none(type, external)

      class(MEF90DefMechSplit), intent(IN) :: self
      class(MEF90HookesLaw), intent(IN)    :: HookesLaw
      class(mef90Mat), intent(IN)          :: phi, psi
      PetscReal, intent(OUT)               :: D2EEDPlus, D2EEDMinus
      PetscErrorCode, intent(inout)        :: ierr
   end subroutine D2EEDInterface
end interface

contains
#undef __FUNCT__
#define __FUNCT__ "MEF90DefMechSplit_SmoothPositiveSquare"
!!! author: Blaise Bourdin (2018-2020, bourdin@lsu.edu)
!!! author: Blaise Bourdin (2026, bourdin@mcmaster.ca)
!!!
!!!  MEF90DefMechSplit_SmoothPositiveSquare: a low order polynomial C^2 regularization of (max(x,0))**2, defined by
!!!       0                         if x \le -gamma/2
!!!       (x+\gamma/2)^3/3/\gamma   if -\gamma/2 < x \le \gamma/2
!!!       x^2+gamma^2/12            otherwise
!!!

function MEF90DefMechSplit_SmoothPositiveSquare(x, gamma)
   PetscReal, intent(IN) :: x
   PetscReal, intent(IN) :: gamma
   PetscReal             :: MEF90DefMechSplit_SmoothPositiveSquare

   PetscReal             :: gammaOver2

   gammaOver2 = gamma * 0.5_kr
   if (x <= -gammaOver2) then
      MEF90DefMechSplit_SmoothPositiveSquare = 0.0_kr
   else if (x <= gammaOver2) then
      MEF90DefMechSplit_SmoothPositiveSquare = (x + gammaOver2)**3 / 3.0_kr / gamma
   else
      MEF90DefMechSplit_SmoothPositiveSquare = x**2 + gammaOver2**2 / 3.0_kr
   end if
end function MEF90DefMechSplit_SmoothPositiveSquare

#undef __FUNCT__
#define __FUNCT__ "MEF90DefMechSplit_DSmoothPositiveSquare"
!!! author: Blaise Bourdin (2018-2020, bourdin@lsu.edu)
!!! author: Blaise Bourdin (2026, bourdin@mcmaster.ca)
!!!
!!!  MEF90DefMechSplit_DSmoothPositiveSquare: the first derivative of MEF90_DefMechSplitHD_PenaltyFunction
!!!

function MEF90DefMechSplit_DSmoothPositiveSquare(x, gamma)
   PetscReal, intent(IN) :: x
   PetscReal, intent(IN) :: gamma
   PetscReal             :: MEF90DefMechSplit_DSmoothPositiveSquare

   PetscReal             :: gammaOver2

   gammaOver2 = gamma * 0.5_kr
   if (x <= -gammaOver2) then
      MEF90DefMechSplit_DSmoothPositiveSquare = 0.0_kr
   else if (x <= gammaOver2) then
      MEF90DefMechSplit_DSmoothPositiveSquare = (x + gammaOver2)**2 / gamma
   else
      MEF90DefMechSplit_DSmoothPositiveSquare = 2.0_kr * x
   end if
end function MEF90DefMechSplit_DSmoothPositiveSquare

#undef __FUNCT__
#define __FUNCT__ "MEF90DefMechSplit_D2SmoothPositiveSquare"
!!! author: Blaise Bourdin (2018-2020, bourdin@lsu.edu)
!!! author: Blaise Bourdin (2026, bourdin@mcmaster.ca)
!!!
!!!  MEF90DefMechSplit_D2SmoothPositiveSquare: the second derivative of MEF90_DefMechSplitHD_PenaltyFunction
!!!

function MEF90DefMechSplit_D2SmoothPositiveSquare(x, gamma)
   PetscReal, intent(IN) :: x
   PetscReal, intent(IN) :: gamma
   PetscReal             :: MEF90DefMechSplit_D2SmoothPositiveSquare

   PetscReal             :: gammaOver2

   gammaOver2 = gamma * 0.5_kr
   if (x <= -gammaOver2) then
      MEF90DefMechSplit_D2SmoothPositiveSquare = 0.0_kr
   else if (x <= gammaOver2) then
      MEF90DefMechSplit_D2SmoothPositiveSquare = 1.0_kr + x / gammaOver2
   else
      MEF90DefMechSplit_D2SmoothPositiveSquare = 2.0_kr
   end if
end function MEF90DefMechSplit_D2SmoothPositiveSquare

end module m_MEF90_DefMechSplit_class

