#include "../MEF90/mef90.inc" module m_MEF90_DefMechSplitHD #include "petsc/finclude/petsc.h" use m_MEF90_DefMechSplit_class use m_MEF90_Materials use m_MEF90_HookesLaw implicit none(type) private public :: MEF90DefMechSplitHD type, extends(MEF90DefMechSplit) :: MEF90DefMechSplitHD PetscReal :: gamma = 1.0e-1 contains procedure, pass(self) :: setFromOptions => setFromOptionsHD procedure, pass(self) :: view_internal => viewHD procedure, pass(self) :: setup => setupHD procedure, pass(self) :: EED => EEDHD procedure, pass(self) :: DEED => DEEDHD procedure, pass(self) :: D2EED => D2EEDHD end type MEF90DefMechSplitHD contains #undef __FUNCT__ #define __FUNCT__ "setFromOptionsHD" !!! author: Blaise Bourdin (2020, bourdin@lsu.edu) !!! author: Blaise Bourdin (2026, bourdin@mcmaster.ca) !!! !!! setFromOptionsHD: the default constructor for a MEF90DefMechSplitHD !!! subroutine setFromOptionsHD(self, ierr) class(MEF90DefMechSplitHD), intent(inout) :: self PetscErrorCode, intent(inout) :: ierr PetscInt :: verbose = 0 ! self%damageOrder = 0 self%quadratureOrder = 2 self%type = 'MEF90DefMechSplitHD' PetscCall(PetscOptionsBegin(self%comm, trim(self%prefix) // "split_HydrostaticDeviatoric_", "Options for MEF90DefMechSplitHD_type", "mef90DefMech", ierr)) PetscCall(PetscOptionsReal('-gamma', 'gamma parameter', '[]', self%gamma, self%gamma, PETSC_NULL_BOOL, ierr)) PetscCall(PetscOptionsBool("-hybrid", "Use a hybrid split", "MEF90", PETSC_FALSE, self%isHybrid, PETSC_NULL_BOOL, ierr)) PetscCall(PetscOptionsEnd(ierr)) PetscCall(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, "-verbose", verbose, PETSC_NULL_BOOL, ierr)) if (verbose > 0) then call self%view(PETSC_VIEWER_STDOUT_WORLD,ierr) end if end subroutine setFromOptionsHD #undef __FUNCT__ #define __FUNCT__ "viewHD" !!! author: Blaise Bourdin (2026, bourdin@mcmaster.ca) !!! !!! viewHD: view a MEF90DefMechSplitHD !!! subroutine viewHD(self, viewer, ierr) class(MEF90DefMechSplitHD), intent(in) :: self type(tPetscViewer), intent(in) :: viewer PetscErrorCode,intent(inout) :: ierr character(len=MEF90MXSTRLEN, kind=c_char) :: IOBuffer character(len=MEF90MXSTRLEN, kind=c_char) :: viewerType PetscCall(PetscViewerGetType(viewer, viewerType, ierr)) if (viewerType == 'ascii') then write(IOBuffer, "(A,': Options for MEF90DefMechSplit\n')") trim(self%prefix) // "split" PetscCall(PetscViewerASCIIPrintf(viewer, IOBuffer, ierr)) write(IOBuffer, "(' type: HydrostaticDeviatoric\n')") PetscCall(PetscViewerASCIIPrintf(viewer, IOBuffer, ierr)) write(IOBuffer, "(' hybrid: ', L1, '\n')") self%isHybrid PetscCall(PetscViewerASCIIPrintf(viewer, IOBuffer, ierr)) write(IOBuffer, "(' gamma: ',ES12.5,' []\n')") self%gamma PetscCall(PetscViewerASCIIPrintf(viewer, IOBuffer, ierr)) end if end subroutine viewHD #undef __FUNCT__ #define __FUNCT__ "setupHD" !!! author: Blaise Bourdin (2020, bourdin@lsu.edu) !!! author: Blaise Bourdin (2026, bourdin@mcmaster.ca) !!! !!! setupHD: the setup routine for a MEF90DefMechSplitHD, which does nothing since there is no split !!! subroutine setupHD(self, Strain, ierr) use m_MEF90 implicit none(type, external) class(MEF90DefMechSplitHD), intent(inout) :: self class(mef90Mat), intent(IN) :: Strain PetscErrorCode, intent(inout) :: ierr self%strain = Strain end subroutine setupHD #undef __FUNCT__ #define __FUNCT__ "EEDHD" !!! author: Blaise Bourdin (2020, bourdin@lsu.edu) !!! author: Blaise Bourdin (2026, bourdin@mcmaster.ca) !!! !!! EEDHD: Compute the positive and negative part of the elastic energy density associated with a strain tensor !!! EEDMinus = SmoothPositiveSquare(-trace(e)) AI.I/2/N^2 !!! EEDPlus = Ae.e/2 - EEDMinus !!! with N = 3 for 3D and N = 2 for 2D, where e is the strain tensor, A is the Hookes law tensor !!! subroutine EEDHD(self, HookesLaw, phi, EEDPlus, EEDMinus, ierr) class(MEF90DefMechSplitHD), intent(IN) :: self class(MEF90HookesLaw), intent(IN) :: HookesLaw class(mef90Mat), intent(IN) :: phi PetscReal, intent(OUT) :: EEDPlus, EEDMinus PetscErrorCode, intent(inout) :: ierr PetscReal :: AIIN2 ! AI.I/N^2 PetscReal :: tre ! tr(e) AIIN2 = 0.0_Kr tre = 0.0_Kr select type(HookesLaw) type is (MEF90HookesLawIsotropic2D) AIIN2 = HookesLaw%BulkModulus type is (MEF90HookesLawIsotropic3D) AIIN2 = HookesLaw%BulkModulus class default select type(e => phi) type is (MatS2D) call HookesLaw%multmult(MEF90MatS2DIdentity, MEF90MatS2DIdentity, AIIN2, ierr) AIIN2 = AIIN2 / 4.0_Kr type is (MatS3D) call HookesLaw%multmult(MEF90MatS3DIdentity, MEF90MatS3DIdentity, AIIN2, ierr) AIIN2 = AIIN2 / 9.0_Kr end select ! self%strain end select ! HookesLaw select type(e => phi) type is (MatS2D) tre = trace(e) type is (MatS3D) tre = trace(e) end select ! phi EEDMinus = AIIN2 * MEF90DefMechSplit_SmoothPositiveSquare(-tre, self%gamma) / 2.0_Kr call HookesLaw%multmult(phi, phi, EEDPlus, ierr) EEDPlus = EEDPlus / 2.0_Kr - EEDMinus end subroutine EEDHD #undef __FUNCT__ #define __FUNCT__ "DEEDHD" !!! author: Blaise Bourdin (2020, bourdin@lsu.edu) !!! author: Blaise Bourdin (2026, bourdin@mcmaster.ca) !!! !!! DEEDHD: Compute the directional derivative of the positive and negative part of the elastic energy density in the direction phi at the point defined by the strain tensor !!! subroutine DEEDHD(self, HookesLaw, phi, DEEDPlus, DEEDMinus, ierr) class(MEF90DefMechSplitHD), intent(IN) :: self class(MEF90HookesLaw), intent(IN) :: HookesLaw class(mef90Mat), intent(IN) :: phi PetscReal, intent(OUT) :: DEEDPlus, DEEDMinus PetscErrorCode, intent(inout) :: ierr PetscReal :: AIIN2 ! AI.I/N^2 PetscReal :: trStrain ! tr(self%Strain) PetscReal :: trPhi ! tr(phi) AIIN2 = 0.0_Kr trStrain = 0.0_Kr trPhi = 0.0_Kr select type(HookesLaw) type is (MEF90HookesLawIsotropic2D) AIIN2 = HookesLaw%BulkModulus type is (MEF90HookesLawIsotropic3D) AIIN2 = HookesLaw%BulkModulus class default select type(e => self%strain) type is (MatS2D) call HookesLaw%multmult(MEF90MatS2DIdentity, MEF90MatS2DIdentity, AIIN2, ierr) AIIN2 = AIIN2 / 4.0_Kr type is (MatS3D) call HookesLaw%multmult(MEF90MatS3DIdentity, MEF90MatS3DIdentity, AIIN2, ierr) AIIN2 = AIIN2 / 9.0_Kr end select ! self%strain end select ! HookesLaw select type(e => self%strain) type is (MatS2D) trStrain = trace(e) type is (MatS3D) trStrain = trace(e) end select ! self%strain select type(p =>phi) type is (MatS2D) trPhi = trace(p) type is (MatS3D) trPhi = trace(p) end select DEEDMinus = -AIIN2 * MEF90DefMechSplit_DSmoothPositiveSquare(-trStrain, self%gamma) * trPhi / 2.0_Kr call HookesLaw%multmult(phi, self%strain, DEEDPlus, ierr) DEEDPlus = DEEDPlus - DEEDMinus end subroutine DEEDHD #undef __FUNCT__ #define __FUNCT__ "D2EEDHD" !!! author: Blaise Bourdin (2020, bourdin@lsu.edu) !!! !!! D2EEDHD: Compute the second derivative of the positive and negative part of the elastic energy density (positive and negative stress) !!! subroutine D2EEDHD(self, HookesLaw, phi, psi, D2EEDPlus, D2EEDMinus, ierr) class(MEF90DefMechSplitHD), intent(IN) :: self class(MEF90HookesLaw), intent(IN) :: HookesLaw class(mef90Mat), intent(IN) :: phi, psi PetscReal, intent(OUT) :: D2EEDPlus, D2EEDMinus PetscErrorCode, intent(inout) :: ierr PetscReal :: AIIN2 ! AI.I/N^2 PetscReal :: trStrain ! tr(self%Strain) PetscReal :: trPhi, trPsi ! tr(phi), tr(psi) AIIN2 = 0.0_Kr trStrain = 0.0_Kr trPhi = 0.0_Kr trPsi = 0.0_Kr select type(HookesLaw) type is (MEF90HookesLawIsotropic2D) AIIN2 = HookesLaw%BulkModulus type is (MEF90HookesLawIsotropic3D) AIIN2 = HookesLaw%BulkModulus class default select type(e => self%strain) type is (MatS2D) call HookesLaw%multmult(MEF90MatS2DIdentity, MEF90MatS2DIdentity, AIIN2, ierr) AIIN2 = AIIN2 / 4.0_Kr trStrain = trace(e) type is (MatS3D) call HookesLaw%multmult(MEF90MatS3DIdentity, MEF90MatS3DIdentity, AIIN2, ierr) AIIN2 = AIIN2 / 9.0_Kr trStrain = trace(e) end select ! self%strain end select ! HookesLaw select type(e => self%strain) type is (MatS2D) trStrain = trace(e) type is (MatS3D) trStrain = trace(e) end select ! self%strain select type(p =>phi) type is (MatS2D) trPhi = trace(p) type is (MatS3D) trPhi = trace(p) end select select type(p =>psi) type is (MatS2D) trPsi = trace(p) type is (MatS3D) trPsi = trace(p) end select D2EEDMinus = AIIN2 * MEF90DefMechSplit_D2SmoothPositiveSquare(-trStrain, self%gamma) * trPhi * trPsi / 2.0_Kr call HookesLaw%multmult(phi, psi, D2EEDPlus, ierr) D2EEDPlus = D2EEDPlus - D2EEDMinus end subroutine D2EEDHD end module m_MEF90_DefMechSplitHD