#include "../MEF90/mef90.inc" module m_MEF90_DefMech_class #include "petsc/finclude/petsc.h" use petscsys use m_MEF90_BaseClass use m_MEF90_Ctx use m_MEF90_DMPlex use m_MEF90_LinAlg use, intrinsic :: iso_c_binding implicit none(type) !!! author: Blaise Bourdin (2012-20, bourdin@lsu.edu) !!! author: Blaise Bourdin (2026, bourdin@mcmaster.ca) !!! !!! MEF90DefMech_Type: The class holding the state and the options of a defect mechanics problem !!! enum, bind(c) enumerator :: MEF90DefMech_SolverTypeAltMin = 0, & MEF90DefMech_SolverTypeQuasiNewton1, & MEF90DefMech_SolverTypeQuasiNewton2 end enum character(len=MEF90MXSTRLEN), dimension(6), protected :: MEF90DefMech_SolverTypeList = [character(len=MEF90MXSTRLEN) :: & 'AltMin', & 'QuasiNewton1', & 'QuasiNewton2', & 'MEF90DefMech_SolverType', & '_MEF90DefMech_SolverType', & ''] enum, bind(c) enumerator :: MEF90DefMech_TimeSteppingTypeNULL = 0, & MEF90DefMech_TimeSteppingTypeQuasiStatic end enum character(len=MEF90MXSTRLEN), dimension(5), protected :: MEF90DefMech_TimeSteppingTypeList = [character(len=MEF90MXSTRLEN) :: & 'Null', & 'QuasiStatic', & 'MEF90DefMech_TimeSteppingType', & '_MEF90DefMech_TimeSteppingType', & ''] enum, bind(c) enumerator :: MEF90DefMech_DamageSolverTypeSNES = 0, & MEF90DefMech_DamageSolverTypeTAO end enum character(len=MEF90MXSTRLEN), dimension(5), protected :: MEF90DefMech_DamageSolverTypeList = [character(len=MEF90MXSTRLEN) :: & 'SNES', & 'Tao', & 'MEF90DefMech_DamageSolverType', & '_MEF90DefMech_DamageSolverType', & ''] enum, bind(c) enumerator :: MEF90DefMech_BTTypeNULL = 0, & MEF90DefMech_BTTypeBackward, & MEF90DefMech_BTTypeForward end enum character(len=MEF90MXSTRLEN), dimension(6), protected :: MEF90DefMech_BTTypeList = [character(len=MEF90MXSTRLEN) :: & 'Null', & 'Backward', & 'Forward', & 'MEF90DefMech_BTType', & '_MEF90DefMech_BTType', & ''] enum, bind(c) enumerator :: MEF90DefMech_plasticityTypeNone = 0, & MEF90DefMech_plasticityTypeTresca, & MEF90DefMech_plasticityTypeVonMises, & MEF90DefMech_plasticityTypeVonMisesPlaneTheory, & MEF90DefMech_plasticityTypeCapModel, & MEF90DefMech_plasticityTypeDruckerPragerCapModel, & MEF90DefMech_plasticityTypeVonMises1D, & MEF90DefMech_plasticityTypeHillPlaneTheory, & MEF90DefMech_PlasticityTypeGreen, & MEF90DefMech_PlasticityTypeGurson end enum character(len=MEF90MXSTRLEN), dimension(13), protected :: MEF90DefMech_plasticityTypeList = [character(len=MEF90MXSTRLEN) :: & 'None', & 'Tresca', & 'VonMises', & 'VonMisesPlaneTheory', & 'CapModel', & 'DruckerPragerCapModel', & 'VonMises1D', & 'HillPlaneTheory', & 'Green', & 'Gurson', & 'MEF90DefMech_plasticityType', & '_MEF90DefMech_plasticityType', & ''] enum, bind(c) enumerator :: MEF90DefMech_unilateralContactTypeNone = 0, & MEF90DefMech_unilateralContactTypeHydrostaticDeviatoric, & MEF90DefMech_unilateralContactTypeHydrostatic, & MEF90DefMech_unilateralContactTypeDeviatoric end enum character(len=MEF90MXSTRLEN), dimension(7), protected :: MEF90DefMech_unilateralContactTypeList = [character(len=MEF90MXSTRLEN) :: & 'None', & 'HydrostaticDeviatoric', & 'Hydrostatic', & 'Deviatoric', & 'MEF90DefMech_unilateralContactTypeList', & '_MEF90DefMech_unilateralContactTypeList', & ''] !!! !!! MEF90DefMechGlobalOptions_Type: the problem-wide options of a MEF90DefMech_Type. !!! The values given here are the defaults used by setFromOptions !!! type, extends(MEF90Object) :: MEF90DefMechGlobalOptions_Type PetscEnum :: timeSteppingType = MEF90DefMech_TimeSteppingTypeQuasiStatic PetscEnum :: solverType = MEF90DefMech_SolverTypeAltMin PetscEnum :: damageSolverType = MEF90DefMech_DamageSolverTypeSNES !!! scaling = time (step) scaling law currently CST, Linear, or File PetscEnum :: boundaryDisplacementScaling = MEF90Scaling_Linear PetscEnum :: displacementLowerBoundScaling = MEF90Scaling_CST PetscEnum :: displacementUpperBoundScaling = MEF90Scaling_CST PetscEnum :: cohesiveDisplacementScaling = MEF90Scaling_CST PetscEnum :: boundaryDamageScaling = MEF90Scaling_CST PetscEnum :: bodyForceScaling = MEF90Scaling_Linear PetscEnum :: boundaryForceScaling = MEF90Scaling_Linear PetscEnum :: pressureForceScaling = MEF90Scaling_Linear PetscEnum :: CrackPressureScaling = MEF90Scaling_Linear PetscReal :: damageATol = 1.0e-3_kr PetscInt :: damageMaxIt = 1000_ki PetscInt :: PCLag = 10_ki PetscReal :: SOROmega = 1.0_kr PetscReal :: irrevthres = 0.0_kr PetscBool :: multiPhaseField = PETSC_FALSE PetscEnum :: BTType = MEF90DefMech_BTTypeNULL PetscInt :: BTInterval = -1_ki PetscInt :: BTScope = -1_ki PetscReal :: BTTol = 1.0e-2_kr PetscReal :: plasticStrainATol = 1.0e-4_kr PetscReal :: InjectedVolumeATol = 1.0e-3_kr PetscReal :: dampingCoefficientDisplacement = 0.0_kr PetscReal :: dampingCoefficientDamage = 0.0_kr PetscBool :: temperatureExport = PETSC_FALSE PetscBool :: displacementExport = PETSC_TRUE PetscBool :: damageExport = PETSC_TRUE PetscBool :: stressExport = PETSC_TRUE PetscBool :: plasticStrainExport = PETSC_FALSE PetscBool :: cumulatedPlasticDissipationExport = PETSC_FALSE contains procedure, pass(self) :: view_internal => MEF90DefMechGlobalOptionsView end type MEF90DefMechGlobalOptions_Type type, extends(MEF90Object) :: MEF90DefMechCellSetOptions_Type PetscReal, dimension(3) :: bodyforce = [0.0_kr, 0.0_kr, 0.0_kr] PetscReal :: crackPressure = 0.0_kr PetscEnum :: plasticityType = MEF90DefMech_plasticityTypeNone PetscReal, dimension(3) :: cohesiveDisplacement = [0.0_kr, 0.0_kr, 0.0_kr] PetscBool, dimension(3) :: Has_displacementBC = [PETSC_FALSE, PETSC_FALSE, PETSC_FALSE] PetscReal, dimension(3) :: boundaryDisplacement = [0.0_kr, 0.0_kr, 0.0_kr] PetscReal, dimension(3) :: displacementLowerBound = [MEF90NINFINITY, MEF90NINFINITY, MEF90NINFINITY] PetscReal, dimension(3) :: displacementUpperBound = [MEF90INFINITY, MEF90INFINITY, MEF90INFINITY] PetscBool :: Has_damageBC = PETSC_FALSE PetscReal :: boundaryDamage = 0.0_kr PetscBool :: CrackVolumeControlled = PETSC_FALSE PetscBool :: WorkControlled = PETSC_FALSE !!! Material properties. linearThermalExpansion is a symmetric matrix, so it is dimension !!! dependent and allocated by MEF90DefMechCellSetOptionsSetFromOptions class(mef90Mat), allocatable :: linearThermalExpansion PetscReal :: cohesiveStiffness = 0.0_kr contains procedure, pass(self) :: view_internal => MEF90DefMechCellSetOptionsView end type MEF90DefMechCellSetOptions_Type type, extends(MEF90Object) :: MEF90DefMechFaceSetOptions_Type PetscReal, dimension(3) :: boundaryforce = [0.0_kr, 0.0_kr, 0.0_kr] PetscReal :: pressureForce = 0.0_kr PetscBool, dimension(3) :: Has_displacementBC = [PETSC_FALSE, PETSC_FALSE, PETSC_FALSE] PetscReal, dimension(3) :: boundaryDisplacement = [0.0_kr, 0.0_kr, 0.0_kr] PetscReal, dimension(3) :: displacementLowerBound = [MEF90NINFINITY, MEF90NINFINITY, MEF90NINFINITY] PetscReal, dimension(3) :: displacementUpperBound = [MEF90INFINITY, MEF90INFINITY, MEF90INFINITY] PetscBool :: Has_damageBC = PETSC_FALSE PetscReal :: boundaryDamage = 0.0_kr contains procedure, pass(self) :: view_internal => MEF90DefMechFaceSetOptionsView end type MEF90DefMechFaceSetOptions_Type type, extends(MEF90Object) :: MEF90DefMechVertexSetOptions_Type PetscBool, dimension(3) :: Has_displacementBC = [PETSC_FALSE, PETSC_FALSE, PETSC_FALSE] PetscReal, dimension(3) :: boundaryDisplacement = [0.0_kr, 0.0_kr, 0.0_kr] PetscReal, dimension(3) :: displacementLowerBound = [MEF90NINFINITY, MEF90NINFINITY, MEF90NINFINITY] PetscReal, dimension(3) :: displacementUpperBound = [MEF90INFINITY, MEF90INFINITY, MEF90INFINITY] PetscBool :: Has_damageBC = PETSC_FALSE PetscReal :: boundaryDamage = 0.0_kr contains procedure, pass(self) :: view_internal => MEF90DefMechVertexSetOptionsView end type MEF90DefMechVertexSetOptions_Type type, extends(MEF90Object) :: MEF90DefMech_Type type(MEF90Ctx_Type), pointer :: MEF90Ctx => null() type(tDM) :: megaDM PetscInt :: dim PetscReal :: analysisTime, timeStep !!! vertex based vec type(tVec), pointer :: displacementLocal => null() type(tVec), pointer :: displacementPreviousStepLocal => null() type(tVec), pointer :: damageLocal => null() type(tVec), dimension(:), pointer :: partialDamageLocal => null() type(tVec), pointer :: displacementLowerBoundLocal => null() type(tVec), pointer :: displacementUpperBoundLocal => null() type(tVec), pointer :: temperatureLocal => null() !!! cell based vec type(tVec), pointer :: bodyForce => null() type(tVec), pointer :: boundaryForce => null() type(tVec), pointer :: pressureForce => null() !Type(tVec),pointer :: crackPressure => null() type(tVec), pointer :: cohesiveDisplacement => null() type(tVec), pointer :: plasticStrain => null() type(tVec), pointer :: cumulatedPlasticDissipation => null() type(tVec), pointer :: stress => null() type(tPetscSF) :: displacementToIOSF, IOToDisplacementSF type(tPetscSF) :: cohesiveDisplacementToIOSF, IOToCohesiveDisplacementSF type(tPetscSF) :: displacementConstraintsSF type(tPetscSF) :: damageToIOSF, IOToDamageSF type(tPetscSF) :: damageConstraintsSF type(tPetscSF) :: temperatureToIOSF, IOToTemperatureSF type(tPetscSF) :: bodyForceToIOSF, IOToBodyForceSF type(tPetscSF) :: boundaryForceToIOSF, IOToBoundaryForceSF type(tPetscSF) :: pressureForceToIOSF, IOToPressureForceSF type(tPetscSF) :: stressToIOSF, IOToStressSF type(tPetscSF) :: plasticStrainToIOSF, IOToPlasticStrainSF type(tPetscSF) :: cumulatedPlasticDissToIOSF, IOToCumulatedPlasticDissSF !!! Neither the problem-wide nor the per-set options are stored: they are read from the options !!! database where they are needed, with MEF90DefMechGlobalOptionsSetFromOptions and !!! MEF90DefMech[Cell,Face,Vertex]SetOptionsSetFromOptions, the same way the AT model, the !!! energy split, and the Hooke's law are obtained. The number of sets is obtained from the !!! megaDM with MEF90DMGetNumSets. type(tPetscViewer) :: globalEnergyViewer type(tPetscViewer), dimension(:), pointer :: setEnergyViewer => null() PetscBool :: hasDisplacementBounds = PETSC_FALSE PetscBool :: hasUnilateralContact = PETSC_FALSE PetscInt :: currentSet = 0_Ki ! used to pass the PF number in multiPhaseField !!! Handle on self, set once in MEF90DefMechCreate with PETScCtx = c_loc(DefMech), and handed to !!! PETSc wherever an application context is expected: SNESSetFunction, SNESSetJacobian, !!! TAOSetObjective, ... The callbacks in m_MEF90_DefMech recover the context with !!! call c_f_pointer(PETScCtx, MEF90DefMechCtx) !!! !!! MEF90DefMech_Type cannot be passed to PETSc directly, the way MEF90DefMechCtx_Type was up to !!! mef90 0.5.2, because extending MEF90Object gives it type-bound procedures and PETSc declares !!! its context arguments as assumed-type. This has to be a component rather than a local or an !!! inline c_loc(), for lifetime reasons that are easy to get wrong: see the note at the top of !!! MEF90/m_MEF90_BaseClass.F90 before changing any of this. type(c_ptr) :: PETScCtx = C_NULL_PTR contains procedure, pass(self) :: setFromOptions => MEF90DefMechSetFromOptions procedure, pass(self) :: view_internal => MEF90DefMechView end type MEF90DefMech_Type contains #undef __FUNCT__ #define __FUNCT__ "MEF90DefMechCreate" !!! author: Blaise Bourdin (2012-20, bourdin@lsu.edu) !!! author: Blaise Bourdin (2026, bourdin@mcmaster.ca) !!! !!! MEF90DefMechCreate: allocates all the fields of a MEF90DefMech_Type !!! subroutine MEF90DefMechCreate(DefMechCtx, dm, MEF90Ctx, prefix, ierr) !!! DefMech has the target attribute so that DefMech%PETScCtx can be handed over to PETSc as an !!! application context. The actual argument MUST therefore also have the target attribute. type(MEF90DefMech_Type), target, intent(OUT) :: DefMechCtx type(tDM), target, intent(IN) :: dm type(MEF90Ctx_Type), target, intent(IN) :: MEF90Ctx character(len=*), intent(IN) :: prefix PetscErrorCode, intent(INOUT) :: ierr type(MEF90CtxGlobalOptions_Type) :: MEF90CtxGlobalOptions type(MEF90DefMechGlobalOptions_Type) :: DefMechGlobalOptions type(tIS) :: setIS PetscInt :: set PetscInt, dimension(:), pointer :: setID character(len=MEF90MXSTRLEN) :: filename, IOBuffer character(len=MEF90MXSTRLEN) :: vecName type(tDM), dimension(:), pointer :: dmList type(tPetscSF) :: dummySF DefMechCtx%MEF90Ctx => MEF90Ctx DefMechCtx%comm = MEF90Ctx%comm DefMechCtx%prefix = prefix DefMechCtx%name = trim(prefix)//"DefMech" DefMechCtx%PETScCtx = c_loc(DefMechCtx) !! !! Create energy viewers !! filename = trim(MEF90FilePrefix(MEF90Ctx%resultFile))//'.ener' PetscCall(PetscViewerASCIIOpen(MEF90Ctx%comm, filename, DefMechCtx%globalEnergyViewer, ierr)) PetscCall(PetscViewerASCIIPrintf(DefMechCtx%globalEnergyViewer, "# step load elastic energy work cohesive energy surface energy total energy plastic dissipation \n", ierr)) PetscCall(PetscViewerFlush(DefMechCtx%globalEnergyViewer, ierr)) PetscCall(DMGetLabelIdIS(dm, MEF90CellSetLabelName, setIS, ierr)) PetscCall(MEF90ISAllGatherMerge(MEF90Ctx%comm, setIS, ierr)) PetscCall(ISGetIndices(setIS, setID, ierr)) allocate (DefMechCtx%setEnergyViewer(size(setID)), stat=ierr) do set = 1, size(setID) write (filename, 101) trim(MEF90FilePrefix(MEF90Ctx%resultFile)), setID(set) PetscCall(PetscViewerASCIIOpen(MEF90Ctx%comm, filename, DefMechCtx%setEnergyViewer(set), ierr)) write (IOBuffer, 102) setID(set) PetscCall(PetscViewerASCIIPrintf(DefMechCtx%setEnergyViewer(set), IOBuffer, ierr)) PetscCall(PetscViewerASCIIPrintf(DefMechCtx%setEnergyViewer(set), "# step load elastic energy work cohesive energy surface energy total energy plastic dissipation\n", ierr)) PetscCall(PetscViewerFlush(DefMechCtx%setEnergyViewer(set), ierr)) end do 101 format(A, '-', I4.4, '.enerblk') 102 format("# cell set ", I4, "\n") DefMechCtx%analysisTime = 0.0_kr DefMechCtx%timeStep = 0.0_kr !! Create Vecs and SF PetscCall(DMGetDimension(dm, DefMechCtx%dim, ierr)) PetscCall(MEF90CtxGlobalOptionsSetFromOptions(DefMechCtx%MEF90Ctx%comm, trim(DefMechCtx%MEF90Ctx%prefix), MEF90CtxGlobalOptions, ierr)) PetscCall(MEF90DefMechGlobalOptionsSetFromOptions(DefMechCtx%comm, trim(DefMechCtx%prefix), DefMechGlobalOptions, ierr)) vecName = "Displacement" allocate (DefMechCtx%displacementLocal, stat=ierr) PetscCall(MEF90CreateLocalVector(dm, MEF90CtxGlobalOptions%elementFamily, MEF90CtxGlobalOptions%elementOrder, DefMechCtx%dim, vecName, DefMechCtx%displacementLocal, ierr)) allocate (DefMechCtx%displacementPreviousStepLocal, stat=ierr) PetscCall(VecDuplicate(DefMechCtx%displacementLocal, DefMechCtx%displacementPreviousStepLocal, ierr)) PetscCall(PetscObjectSetName(DefMechCtx%displacementPreviousStepLocal, "DisplacementPreviousStep", ierr)) allocate (DefMechCtx%displacementLowerBoundLocal, stat=ierr) PetscCall(VecDuplicate(DefMechCtx%displacementLocal, DefMechCtx%displacementLowerBoundLocal, ierr)) PetscCall(PetscObjectSetName(DefMechCtx%displacementLowerBoundLocal, "DisplacementLowerBound", ierr)) allocate (DefMechCtx%displacementUpperBoundLocal, stat=ierr) PetscCall(VecDuplicate(DefMechCtx%displacementLocal, DefMechCtx%displacementUpperBoundLocal, ierr)) PetscCall(PetscObjectSetName(DefMechCtx%displacementUpperBoundLocal, "DisplacementUpperBound", ierr)) allocate (DefMechCtx%damageLocal, stat=ierr) vecName = "Damage" PetscCall(MEF90CreateLocalVector(dm, MEF90CtxGlobalOptions%elementFamily, MEF90CtxGlobalOptions%elementOrder, 1_ki, vecName, DefMechCtx%damageLocal, ierr)) if (DefMechGlobalOptions%multiPhaseField) then allocate (DefMechCtx%partialDamageLocal(size(setID)), stat=ierr) do set = 1, size(DefMechCtx%partialDamageLocal) !! I need to cheat here: !! I create the Vec with name "Damage" so that it inherits the proper damage BC, then change its name to the proper value vecName = "Damage" PetscCall(MEF90CreateLocalVector(dm, MEF90CtxGlobalOptions%elementFamily, MEF90CtxGlobalOptions%elementOrder, 1_ki, vecName, DefMechCtx%partialDamageLocal(set), ierr)) write(Vecname,'("partialDamage-", I4.4)') set PetscCall(PetscObjectSetName(DefMechCtx%partialDamageLocal(set), vecName, ierr)) end do end if allocate (DefMechCtx%TemperatureLocal, stat=ierr) PetscCall(VecDuplicate(DefMechCtx%damageLocal, DefMechCtx%TemperatureLocal, ierr)) PetscCall(PetscObjectSetName(DefMechCtx%TemperatureLocal, "Temperature", ierr)) vecName = "cohesiveDisplacement" allocate (DefMechCtx%cohesiveDisplacement, stat=ierr) PetscCall(MEF90CreateLocalVector(dm, MEF90CtxGlobalOptions%elementFamily, MEF90CtxGlobalOptions%elementOrder, DefMechCtx%dim, vecName, DefMechCtx%cohesiveDisplacement, ierr)) vecName = "bodyForce" allocate (DefMechCtx%bodyForce, stat=ierr) PetscCall(MEF90CreateCellVector(dm, DefMechCtx%dim, vecName, DefMechCtx%bodyForce, ierr)) vecName = "boundaryForce" allocate (DefMechCtx%boundaryForce, stat=ierr) PetscCall(MEF90CreateBoundaryCellVector(dm, DefMechCtx%dim, vecName, DefMechCtx%boundaryForce, ierr)) vecName = "pressureForce" allocate (DefMechCtx%pressureForce, stat=ierr) PetscCall(MEF90CreateBoundaryCellVector(dm, 1_ki, vecName, DefMechCtx%pressureForce, ierr)) vecName = "plasticStrain" allocate (DefMechCtx%plasticStrain, stat=ierr) PetscCall(MEF90CreateCellVector(dm, (DefMechCtx%dim * (DefMechCtx%dim + 1_ki)) / 2_ki, vecName, DefMechCtx%plasticStrain, ierr)) allocate (DefMechCtx%cumulatedPlasticDissipation, stat=ierr) PetscCall(VecDuplicate(DefMechCtx%plasticStrain, DefMechCtx%cumulatedPlasticDissipation, ierr)) PetscCall(PetscObjectSetName(DefMechCtx%cumulatedPlasticDissipation, "cumulatedPlasticDissipation", ierr)) allocate (DefMechCtx%stress, stat=ierr) PetscCall(VecDuplicate(DefMechCtx%plasticStrain, DefMechCtx%stress, ierr)) PetscCall(PetscObjectSetName(DefMechCtx%stress, "Stress", ierr)) !! Create megaDM !! This needs to be modified to add the individual damage fields if needed allocate (dmList(7)) PetscCall(VecGetDM(DefMechCtx%displacementLocal, dmList(1), ierr)) PetscCall(VecGetDM(DefMechCtx%damageLocal, dmList(2), ierr)) PetscCall(VecGetDM(DefMechCtx%temperatureLocal, dmList(3), ierr)) PetscCall(VecGetDM(DefMechCtx%bodyForce, dmList(4), ierr)) PetscCall(VecGetDM(DefMechCtx%boundaryForce, dmList(5), ierr)) PetscCall(VecGetDM(DefMechCtx%pressureForce, dmList(6), ierr)) PetscCall(VecGetDM(DefMechCtx%plasticStrain, dmList(7), ierr)) PetscCall(DMCreateSuperDM(dmList, 7_ki, PETSC_NULL_IS_POINTER, DefMechCtx%megaDM, ierr)) deallocate (dmList) !! Create the IO SF for all fields PetscCall(MEF90IOSFCreate(MEF90Ctx, DefMechCtx%displacementLocal, DefMechCtx%displacementToIOSF, DefMechCtx%IOTodisplacementSF, ierr)) PetscCall(MEF90IOSFCreate(MEF90Ctx, DefMechCtx%damageLocal, DefMechCtx%damageToIOSF, DefMechCtx%IOTodamageSF, ierr)) PetscCall(MEF90IOSFCreate(MEF90Ctx, DefMechCtx%cohesiveDisplacement, DefMechCtx%cohesiveDisplacementToIOSF, DefMechCtx%IOToCohesiveDisplacementSF, ierr)) PetscCall(MEF90IOSFCreate(MEF90Ctx, DefMechCtx%temperatureLocal, DefMechCtx%temperatureToIOSF, DefMechCtx%IOTotemperatureSF, ierr)) PetscCall(MEF90IOSFCreate(MEF90Ctx, DefMechCtx%bodyForce, DefMechCtx%bodyForceToIOSF, DefMechCtx%IOTobodyForceSF, ierr)) ! PetscCall(MEF90IOSFCreate(MEF90Ctx,DefMech%boundaryForce,DefMech%boundaryForceToIOSF,DefMech%IOToboundaryForceSF,ierr)) ! PetscCall(MEF90IOSFCreate(MEF90Ctx,DefMech%pressureForce,DefMech%pressureForceToIOSF,DefMech%IOTopressureForceSF,ierr)) PetscCall(MEF90IOSFCreate(MEF90Ctx, DefMechCtx%stress, DefMechCtx%stressToIOSF, DefMechCtx%IOToStressSF, ierr)) PetscCall(MEF90IOSFCreate(MEF90Ctx, DefMechCtx%plasticStrain, DefMechCtx%plasticStrainToIOSF, DefMechCtx%IOToplasticStrainSF, ierr)) PetscCall(MEF90IOSFCreate(MEF90Ctx, DefMechCtx%cumulatedPlasticDissipation, DefMechCtx%cumulatedPlasticDissToIOSF, DefMechCtx%IOToCumulatedPlasticDissSF, ierr)) !! Create the SF to exchange boundary values of the displacement and damage. PetscCall(MEF90ConstraintSFCreate(DefMechCtx%MEF90Ctx, DefMechCtx%displacementLocal, DefMechCtx%displacementLocal, DefMechCtx%displacementConstraintsSF, dummySF, ierr)) PetscCall(PetscSFDestroy(dummySF, ierr)) PetscCall(MEF90ConstraintSFCreate(DefMechCtx%MEF90Ctx, DefMechCtx%damageLocal, DefMechCtx%damageLocal, DefMechCtx%damageConstraintsSF, dummySF, ierr)) PetscCall(PetscSFDestroy(dummySF, ierr)) PetscCall(ISRestoreIndices(setIS, setID, ierr)) PetscCall(ISDestroy(setIS, ierr)) end subroutine MEF90DefMechCreate #undef __FUNCT__ #define __FUNCT__ "MEF90DefMechDestroy" !!! author: Blaise Bourdin (2012-14, bourdin@lsu.edu) !!! author: Blaise Bourdin (2026, bourdin@mcmaster.ca) !!! !!! MEF90DefMechDestroy: destroys a MEF90DefMech_Type !!! subroutine MEF90DefMechDestroy(DefMech, ierr) type(MEF90DefMech_Type), intent(INOUT) :: DefMech PetscErrorCode, intent(INOUT) :: ierr PetscInt :: set DefMech%PETScCtx = C_NULL_PTR !! !! Close energy viewers !! PetscCall(PetscViewerDestroy(DefMech%globalEnergyViewer, ierr)) do set = 1, size(DefMech%setEnergyViewer) PetscCall(PetscViewerDestroy(DefMech%setEnergyViewer(set), ierr)) end do deallocate (DefMech%setEnergyViewer) !! Destroy Vecs and SF and deAllocate them if (associated(DefMech%displacementLocal)) then PetscCall(VecDestroy(DefMech%displacementLocal, ierr)) deallocate (DefMech%displacementLocal) nullify (DefMech%displacementLocal) end if if (associated(DefMech%displacementPreviousStepLocal)) then PetscCall(VecDestroy(DefMech%displacementPreviousStepLocal, ierr)) deallocate (DefMech%displacementPreviousStepLocal) nullify (DefMech%displacementPreviousStepLocal) end if if (associated(DefMech%displacementLowerBoundLocal)) then PetscCall(VecDestroy(DefMech%displacementLowerBoundLocal, ierr)) deallocate (DefMech%displacementLowerBoundLocal) nullify (DefMech%displacementLowerBoundLocal) end if if (associated(DefMech%displacementUpperBoundLocal)) then PetscCall(VecDestroy(DefMech%displacementUpperBoundLocal, ierr)) deallocate (DefMech%displacementUpperBoundLocal) nullify (DefMech%displacementUpperBoundLocal) end if if (associated(DefMech%damageLocal)) then PetscCall(VecDestroy(DefMech%damageLocal, ierr)) deallocate (DefMech%damageLocal) nullify (DefMech%damageLocal) end if if (associated(DefMech%partialDamageLocal)) then do set = 1, size(DefMech%partialDamageLocal) PetscCall(VecDestroy(DefMech%partialDamageLocal(set), ierr)) end do deallocate (DefMech%partialDamageLocal) nullify (DefMech%partialDamageLocal) end if if (associated(DefMech%temperatureLocal)) then PetscCall(VecDestroy(DefMech%temperatureLocal, ierr)) deallocate (DefMech%temperatureLocal) nullify (DefMech%temperatureLocal) end if if (associated(DefMech%cohesiveDisplacement)) then PetscCall(VecDestroy(DefMech%cohesiveDisplacement, ierr)) deallocate (DefMech%cohesiveDisplacement) nullify (DefMech%cohesiveDisplacement) end if if (associated(DefMech%bodyForce)) then PetscCall(VecDestroy(DefMech%bodyForce, ierr)) deallocate (DefMech%bodyForce) nullify (DefMech%bodyForce) end if if (associated(DefMech%boundaryForce)) then PetscCall(VecDestroy(DefMech%boundaryForce, ierr)) deallocate (DefMech%boundaryForce) nullify (DefMech%boundaryForce) end if if (associated(DefMech%pressureForce)) then PetscCall(VecDestroy(DefMech%pressureForce, ierr)) deallocate (DefMech%pressureForce) nullify (DefMech%pressureForce) end if if (associated(DefMech%plasticStrain)) then PetscCall(VecDestroy(DefMech%plasticStrain, ierr)) deallocate (DefMech%plasticStrain) nullify (DefMech%plasticStrain) end if if (associated(DefMech%cumulatedPlasticDissipation)) then PetscCall(VecDestroy(DefMech%cumulatedPlasticDissipation, ierr)) deallocate (DefMech%cumulatedPlasticDissipation) nullify (DefMech%cumulatedPlasticDissipation) end if if (associated(DefMech%Stress)) then PetscCall(VecDestroy(DefMech%Stress, ierr)) deallocate (DefMech%Stress) nullify (DefMech%Stress) end if !! Destroy all PetscSF PetscCall(PetscSFDestroy(DefMech%displacementToIOSF, ierr)) PetscCall(PetscSFDestroy(DefMech%IOTodisplacementSF, ierr)) PetscCall(PetscSFDestroy(DefMech%cohesiveDisplacementToIOSF, ierr)) PetscCall(PetscSFDestroy(DefMech%IOToCohesiveDisplacementSF, ierr)) PetscCall(PetscSFDestroy(DefMech%damageToIOSF, ierr)) PetscCall(PetscSFDestroy(DefMech%IOTodamageSF, ierr)) PetscCall(PetscSFDestroy(DefMech%temperatureToIOSF, ierr)) PetscCall(PetscSFDestroy(DefMech%IOTotemperatureSF, ierr)) PetscCall(PetscSFDestroy(DefMech%bodyForceToIOSF, ierr)) PetscCall(PetscSFDestroy(DefMech%IOTobodyForceSF, ierr)) PetscCall(PetscSFDestroy(DefMech%boundaryForceToIOSF, ierr)) PetscCall(PetscSFDestroy(DefMech%IOToboundaryForceSF, ierr)) PetscCall(PetscSFDestroy(DefMech%pressureForceToIOSF, ierr)) PetscCall(PetscSFDestroy(DefMech%IOTopressureForceSF, ierr)) PetscCall(PetscSFDestroy(DefMech%stressToIOSF, ierr)) PetscCall(PetscSFDestroy(DefMech%IOToStressSF, ierr)) PetscCall(PetscSFDestroy(DefMech%plasticStrainToIOSF, ierr)) PetscCall(PetscSFDestroy(DefMech%IOToplasticStrainSF, ierr)) PetscCall(PetscSFDestroy(DefMech%cumulatedPlasticDissToIOSF, ierr)) PetscCall(PetscSFDestroy(DefMech%IOToCumulatedPlasticDissSF, ierr)) !! Destroy the SF to exchange boundary values of the displacement and damage. PetscCall(PetscSFDestroy(DefMech%displacementConstraintsSF, ierr)) PetscCall(PetscSFDestroy(DefMech%damageConstraintsSF, ierr)) !! Destroy the megaDM PetscCall(DMDestroy(DefMech%megaDM, ierr)) end subroutine MEF90DefMechDestroy #undef __FUNCT__ #define __FUNCT__ "MEF90DefMechGlobalOptionsSetFromOptions" !!! author: Blaise Bourdin (2026, bourdin@mcmaster.ca) !!! !!! MEF90DefMechGlobalOptionsSetFromOptions: reads the problem-wide options of a MEF90DefMech_Type !!! subroutine MEF90DefMechGlobalOptionsSetFromOptions(comm, prefix, options, ierr) MPIU_Comm, intent(IN) :: comm character(len=*), intent(IN) :: prefix type(MEF90DefMechGlobalOptions_Type), intent(OUT) :: options PetscErrorCode, intent(INOUT) :: ierr PetscInt :: verbose = 0 options%comm = comm options%prefix = prefix options%name = trim(prefix)//"DefMechGlobalOptions" PetscCall(PetscOptionsBegin(comm, prefix, "Options for MEF90DefMech_Type", "mef90DefMech", ierr)) PetscCall(PetscOptionsEnum('-DefMech_TimeStepping_Type', 'Type of defect mechanics Time steping', 'mef90DefMech', MEF90DefMech_TimeSteppingTypeList, options%timeSteppingType, options%timeSteppingType, PETSC_NULL_BOOL, ierr)) PetscCall(PetscOptionsEnum('-DefMech_solver_Type', 'Type of defect mechanics solver', 'mef90DefMech', MEF90DefMech_SolverTypeList, options%solverType, options%solverType, PETSC_NULL_BOOL, ierr)) PetscCall(PetscOptionsEnum('-DefMech_damageSolver_Type', 'Type of defect mechanics damage solver', 'mef90DefMech', MEF90DefMech_DamageSolverTypeList, options%damageSolverType, options%damageSolverType, PETSC_NULL_BOOL, ierr)) PetscCall(PetscOptionsBool('-temperature_export', 'Export temperature in result file', 'mef90DefMech', options%temperatureExport, options%temperatureExport, PETSC_NULL_BOOL, ierr)) PetscCall(PetscOptionsBool('-displacement_export', 'Export displacement in result file', 'mef90DefMech', options%displacementExport, options%displacementExport, PETSC_NULL_BOOL, ierr)) PetscCall(PetscOptionsBool('-damage_export', 'Export damage in result file', 'mef90DefMech', options%damageExport, options%damageExport, PETSC_NULL_BOOL, ierr)) PetscCall(PetscOptionsBool('-stress_export', 'Export stress in result file', 'mef90DefMech', options%stressExport, options%stressExport, PETSC_NULL_BOOL, ierr)) PetscCall(PetscOptionsBool('-plasticstrain_export', 'Export plastic strain in result file', 'mef90DefMech', options%plasticStrainExport, options%plasticStrainExport, PETSC_NULL_BOOL, ierr)) PetscCall(PetscOptionsBool('-cumulatedplasticdissipation_export', 'Export cumulated plastic dissipation in result file', 'mef90DefMech', options%cumulatedPlasticDissipationExport, options%cumulatedPlasticDissipationExport, PETSC_NULL_BOOL, ierr)) PetscCall(PetscOptionsEnum('-boundaryDisplacement_scaling', 'Boundary displacement scaling', 'mef90DefMech', MEF90ScalingList, options%boundaryDisplacementScaling, options%boundaryDisplacementScaling, PETSC_NULL_BOOL, ierr)) PetscCall(PetscOptionsEnum('-displacementlowerbound_scaling', 'Displacement lower bound scaling', 'mef90DefMech', MEF90ScalingList, options%displacementLowerBoundScaling, options%displacementLowerBoundScaling, PETSC_NULL_BOOL, ierr)) PetscCall(PetscOptionsEnum('-displacementupperbound_scaling', 'Displacement upper bound scaling', 'mef90DefMech', MEF90ScalingList, options%displacementUpperBoundScaling, options%displacementUpperBoundScaling, PETSC_NULL_BOOL, ierr)) PetscCall(PetscOptionsEnum('-cohesiveDisplacement_scaling', 'Cohesive displacement scaling', 'mef90DefMech', MEF90ScalingList, options%cohesiveDisplacementScaling, options%cohesiveDisplacementScaling, PETSC_NULL_BOOL, ierr)) PetscCall(PetscOptionsEnum('-boundaryDamage_scaling', 'Boundary damage scaling', 'mef90DefMech', MEF90ScalingList, options%boundaryDamageScaling, options%boundaryDamageScaling, PETSC_NULL_BOOL, ierr)) PetscCall(PetscOptionsEnum('-bodyforce_scaling', 'Body force scaling', 'mef90DefMech', MEF90ScalingList, options%bodyForceScaling, options%bodyForceScaling, PETSC_NULL_BOOL, ierr)) PetscCall(PetscOptionsEnum('-boundaryforce_scaling', 'Boundary force scaling', 'mef90DefMech', MEF90ScalingList, options%boundaryForceScaling, options%boundaryForceScaling, PETSC_NULL_BOOL, ierr)) PetscCall(PetscOptionsEnum('-pressureForce_scaling', 'Pressure force scaling', 'mef90DefMech', MEF90ScalingList, options%pressureForceScaling, options%pressureForceScaling, PETSC_NULL_BOOL, ierr)) PetscCall(PetscOptionsEnum('-crackPressure_scaling', 'Crack Pressure scaling', 'mef90DefMech', MEF90ScalingList, options%CrackPressureScaling, options%CrackPressureScaling, PETSC_NULL_BOOL, ierr)) PetscCall(PetscOptionsReal('-defmech_damage_atol', 'Absolute tolerance on damage error', 'mef90DefMech', options%damageATol, options%damageATol, PETSC_NULL_BOOL, ierr)) PetscCall(PetscOptionsInt('-defmech_damage_maxit', 'Maximum number of alternate minimizations for damage', 'mef90DefMech', options%damageMaxIt, options%damageMaxIt, PETSC_NULL_BOOL, ierr)) PetscCall(PetscOptionsInt('-defmech_pclag', 'Interval at which the PC is recomputed during alternate minimization', 'mef90DefMech', options%PCLag, options%PCLag, PETSC_NULL_BOOL, ierr)) PetscCall(PetscOptionsReal('-defmech_SOR_Omega', 'Alterate Minimization over relaxation factor (>0 for limited, <0 for projected)', 'mef90DefMech', options%SOROmega, options%SOROmega, PETSC_NULL_BOOL, ierr)) PetscCall(PetscOptionsReal('-defmech_irrevThres', 'Threshold above which irreversibility is enforced (0 for monotonicity, .99 for equality)', 'mef90DefMech', options%irrevthres, options%irrevthres, PETSC_NULL_BOOL, ierr)) PetscCall(PetscOptionsBool('-multiPhaseField', 'Use one damage variable per cell set', 'mef90DefMech', options%multiPhaseField, options%multiPhaseField, PETSC_NULL_BOOL, ierr)) PetscCall(PetscOptionsEnum('-BT_Type', 'Backtracking type', 'mef90DefMech', MEF90DefMech_BTTypeList, options%BTType, options%BTType, PETSC_NULL_BOOL, ierr)) PetscCall(PetscOptionsInt('-BT_Interval', 'Interval at which Backtracking is run in inner loop (0 for outer loop)', 'mef90DefMech', options%BTInterval, options%BTInterval, PETSC_NULL_BOOL, ierr)) PetscCall(PetscOptionsInt('-BT_Scope', 'Backtracking scope (0 for unlimited)', 'mef90DefMech', options%BTScope, options%BTScope, PETSC_NULL_BOOL, ierr)) PetscCall(PetscOptionsReal('-BT_Tol', 'Backtracking relative tolerance', 'mef90DefMech', options%BTTol, options%BTTol, PETSC_NULL_BOOL, ierr)) PetscCall(PetscOptionsReal('-defmech_plasticstrain_atol', 'Absolute tolerance on plastic error', 'mef90DefMech', options%plasticStrainATol, options%plasticStrainATol, PETSC_NULL_BOOL, ierr)) PetscCall(PetscOptionsReal('-defmech_InjectedVolume_atol', 'Absolute tolerance on injected volume error', 'mef90DefMech', options%InjectedVolumeATol, options%InjectedVolumeATol, PETSC_NULL_BOOL, ierr)) PetscCall(PetscOptionsReal('-defmech_dampingCoefficient_displacement', 'Damping coefficient on displacement field (0 for minimization, 1 for semi-implicit gradient flow)', 'mef90DefMech', options%dampingCoefficientDisplacement, options%dampingCoefficientDisplacement, PETSC_NULL_BOOL, ierr)) PetscCall(PetscOptionsReal('-defmech_dampingCoefficient_damage', 'Damping coefficient on damage field (0 for minimization, 1 for semi-implicit gradient flow)', 'mef90DefMech', options%dampingCoefficientDamage, options%dampingCoefficientDamage, 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 !! PETSC_VIEWER_STDOUT_WORLD and not PetscViewerASCIIGetStdout: the latter is collective, and the !! options are read where they are needed, which is not always on all ranks. view itself only !! prints the first time a given set is read, so it is collective only then. call options%view(PETSC_VIEWER_STDOUT_WORLD, ierr) end if end subroutine MEF90DefMechGlobalOptionsSetFromOptions #undef __FUNCT__ #define __FUNCT__ "MEF90DefMechGlobalOptionsView" !!! author: Blaise Bourdin (2026, bourdin@mcmaster.ca) !!! !!! MEF90DefMechGlobalOptionsView: the default viewer for a MEF90DefMechGlobalOptions_Type !!! subroutine MEF90DefMechGlobalOptionsView(self, viewer, ierr) class(MEF90DefMechGlobalOptions_Type), 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') return write (IOBuffer, "(A,': Options for MEF90DefMechGlobalOptions_Type\n')") trim(self%prefix)//"defmech" PetscCall(PetscViewerASCIIPrintf(viewer, IOBuffer, ierr)) write (IOBuffer, "(' time stepping type: ',A,'\n')") trim(MEF90DefMech_TimeSteppingTypeList(self%timeSteppingType + 1)) PetscCall(PetscViewerASCIIPrintf(viewer, IOBuffer, ierr)) write (IOBuffer, "(' solver type: ',A,'\n')") trim(MEF90DefMech_SolverTypeList(self%solverType + 1)) PetscCall(PetscViewerASCIIPrintf(viewer, IOBuffer, ierr)) write (IOBuffer, "(' damage solver type: ',A,'\n')") trim(MEF90DefMech_DamageSolverTypeList(self%damageSolverType + 1)) PetscCall(PetscViewerASCIIPrintf(viewer, IOBuffer, ierr)) write (IOBuffer, "(' multi phase field: ',L1,'\n')") self%multiPhaseField PetscCall(PetscViewerASCIIPrintf(viewer, IOBuffer, ierr)) write (IOBuffer, "(' damage atol / maxit: ',ES12.5,' / ',I6,'\n')") self%damageATol, self%damageMaxIt PetscCall(PetscViewerASCIIPrintf(viewer, IOBuffer, ierr)) write (IOBuffer, "(' PC lag: ',I6,' SOR omega: ',ES12.5,' irreversibility threshold: ',ES12.5,'\n')") & self%PCLag, self%SOROmega, self%irrevthres PetscCall(PetscViewerASCIIPrintf(viewer, IOBuffer, ierr)) write (IOBuffer, "(' backtracking type: ',A,' interval: ',I6,' scope: ',I6,' tol: ',ES12.5,'\n')") & trim(MEF90DefMech_BTTypeList(self%BTType + 1)), self%BTInterval, self%BTScope, self%BTTol PetscCall(PetscViewerASCIIPrintf(viewer, IOBuffer, ierr)) write (IOBuffer, "(' plastic strain atol / injected volume atol: ',ES12.5,' / ',ES12.5,'\n')") & self%plasticStrainATol, self%InjectedVolumeATol PetscCall(PetscViewerASCIIPrintf(viewer, IOBuffer, ierr)) write (IOBuffer, "(' damping coefficient displacement / damage: ',ES12.5,' / ',ES12.5,'\n')") & self%dampingCoefficientDisplacement, self%dampingCoefficientDamage PetscCall(PetscViewerASCIIPrintf(viewer, IOBuffer, ierr)) write (IOBuffer, "(' scaling boundary displacement / lower bound / upper bound: ',A,' / ',A,' / ',A,'\n')") & trim(MEF90ScalingList(self%boundaryDisplacementScaling + 1)), & trim(MEF90ScalingList(self%displacementLowerBoundScaling + 1)), & trim(MEF90ScalingList(self%displacementUpperBoundScaling + 1)) PetscCall(PetscViewerASCIIPrintf(viewer, IOBuffer, ierr)) write (IOBuffer, "(' scaling cohesive displacement / boundary damage: ',A,' / ',A,'\n')") & trim(MEF90ScalingList(self%cohesiveDisplacementScaling + 1)), trim(MEF90ScalingList(self%boundaryDamageScaling + 1)) PetscCall(PetscViewerASCIIPrintf(viewer, IOBuffer, ierr)) write (IOBuffer, "(' scaling body force / boundary force / pressure force / crack pressure: ',A,' / ',A,' / ',A,' / ',A,'\n')") & trim(MEF90ScalingList(self%bodyForceScaling + 1)), trim(MEF90ScalingList(self%boundaryForceScaling + 1)), & trim(MEF90ScalingList(self%pressureForceScaling + 1)), trim(MEF90ScalingList(self%CrackPressureScaling + 1)) PetscCall(PetscViewerASCIIPrintf(viewer, IOBuffer, ierr)) write (IOBuffer, "(' export displacement / damage / stress: ',3(L1,' '),'\n')") & self%displacementExport, self%damageExport, self%stressExport PetscCall(PetscViewerASCIIPrintf(viewer, IOBuffer, ierr)) write (IOBuffer, "(' export temperature / plastic strain / cumulated plastic dissipation: ',3(L1,' '),'\n')") & self%temperatureExport, self%plasticStrainExport, self%cumulatedPlasticDissipationExport PetscCall(PetscViewerASCIIPrintf(viewer, IOBuffer, ierr)) end subroutine MEF90DefMechGlobalOptionsView #undef __FUNCT__ #define __FUNCT__ "MEF90DefMechCellSetOptionsSetFromOptions" !!! author: Blaise Bourdin (2026, bourdin@mcmaster.ca) !!! !!! MEF90DefMechCellSetOptionsSetFromOptions: reads the options of a single cell set !!! subroutine MEF90DefMechCellSetOptionsSetFromOptions(comm, prefix, dim, options, ierr) MPIU_Comm, intent(IN) :: comm character(len=*), intent(IN) :: prefix PetscInt, intent(IN) :: dim type(MEF90DefMechCellSetOptions_Type), intent(OUT) :: options PetscErrorCode, intent(INOUT) :: ierr PetscInt :: nOpt PetscReal, dimension(:), allocatable :: tmpArray PetscInt :: verbose = 0 options%comm = comm options%prefix = prefix options%name = trim(prefix)//"DefMechCellSetOptions" select case (dim) case (2) options%linearThermalExpansion = MEF90MatS2DIdentity case (3) options%linearThermalExpansion = MEF90MatS3DIdentity case default SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "dim must be 2 or 3 in " // __FUNCT__) end select PetscCall(PetscOptionsBegin(comm, prefix, "Options for a MEF90DefMech_Type cell set", "mef90DefMech", ierr)) nOpt = 3 PetscCall(PetscOptionsRealArray('-bodyForce', '[N.m^(-3) / N.m^(-2)] (f): body force', 'mef90DefMech', options%bodyForce, nOpt, PETSC_NULL_BOOL, ierr)) PetscCall(PetscOptionsReal('-CrackPressure', 'without unit: internal crack pressure', 'mef90DefMech', options%CrackPressure, options%CrackPressure, PETSC_NULL_BOOL, ierr)) PetscCall(PetscOptionsEnum('-plasticity_type', 'Type of plasticity law', 'mef90DefMech', MEF90DefMech_plasticityTypeList, options%plasticityType, options%plasticityType, PETSC_NULL_BOOL, ierr)) nOpt = 3 PetscCall(PetscOptionsRealArray('-cohesiveDisplacement', '[m] (U): Cohesive displacement value', 'mef90DefMech', options%cohesiveDisplacement, nOpt, PETSC_NULL_BOOL, ierr)) nOpt = 3 PetscCall(PetscOptionsBoolArray('-DisplacementBC', 'Displacement has Dirichlet boundary Condition (Y/N)', 'mef90DefMech', options%Has_displacementBC, nOpt, PETSC_NULL_BOOL, ierr)) nOpt = 3 PetscCall(PetscOptionsRealArray('-boundaryDisplacement', '[m] (U): Displacement boundary value', 'mef90DefMech', options%boundaryDisplacement, nOpt, PETSC_NULL_BOOL, ierr)) nOpt = 3 PetscCall(PetscOptionsRealArray('-displacementLowerBound', '[m] (U): Displacement lower bound', 'mef90DefMech', options%displacementLowerBound, nOpt, PETSC_NULL_BOOL, ierr)) nOpt = 3 PetscCall(PetscOptionsRealArray('-displacementUpperBound', '[m] (U): Displacement upper bound', 'mef90DefMech', options%displacementUpperBound, nOpt, PETSC_NULL_BOOL, ierr)) PetscCall(PetscOptionsBool('-DamageBC', 'Damage has Dirichlet boundary Condition (Y/N)', 'mef90DefMech', options%Has_DamageBC, options%Has_DamageBC, PETSC_NULL_BOOL, ierr)) PetscCall(PetscOptionsBool('-CrackVolumeControlled', 'Crack Pressure controlled by the crack volume in this block (Y/N)', 'mef90DefMech', options%CrackVolumeControlled, options%CrackVolumeControlled, PETSC_NULL_BOOL, ierr)) PetscCall(PetscOptionsBool('-WorkControlled', 'Force magnitude controlled by its work in this block (Y/N)', 'mef90DefMech', options%WorkControlled, options%WorkControlled, PETSC_NULL_BOOL, ierr)) PetscCall(PetscOptionsReal('-boundaryDamage', '[unit-less] (alpha): Damage boundary value', 'mef90DefMech', options%boundaryDamage, options%boundaryDamage, PETSC_NULL_BOOL, ierr)) PetscCall(PetscOptionsReal('-cohesiveStiffness', '[N.m^(-4)] (k) cohesive stiffness in Winkler-type models', 'mef90DefMech', options%cohesiveStiffness, options%cohesiveStiffness, PETSC_NULL_BOOL, ierr)) select type (alpha => options%linearThermalExpansion) type is (MatS2D) nOpt = 3 allocate (tmpArray(nOpt)) tmpArray = alpha PetscCall(PetscOptionsRealArray('-LinearThermalExpansion', '[K^(-1)] (alpha) Linear thermal expansion matrix', 'mef90DefMech', tmpArray, nOpt, PETSC_NULL_BOOL, ierr)) alpha = tmpArray deallocate (tmpArray) type is (MatS3D) nOpt = 6 allocate (tmpArray(nOpt)) tmpArray = alpha PetscCall(PetscOptionsRealArray('-LinearThermalExpansion', '[K^(-1)] (alpha) Linear thermal expansion matrix', 'mef90DefMech', tmpArray, nOpt, PETSC_NULL_BOOL, ierr)) alpha = tmpArray deallocate (tmpArray) end select PetscCall(PetscOptionsEnd(ierr)) PetscCall(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, "-verbose", verbose, PETSC_NULL_BOOL, ierr)) if (verbose > 0) then !! PETSC_VIEWER_STDOUT_WORLD and not PetscViewerASCIIGetStdout: the latter is collective, and the !! options are read where they are needed, which is not always on all ranks. view itself only !! prints the first time a given set is read, so it is collective only then. call options%view(PETSC_VIEWER_STDOUT_WORLD, ierr) end if end subroutine MEF90DefMechCellSetOptionsSetFromOptions #undef __FUNCT__ #define __FUNCT__ "MEF90DefMechCellSetOptionsView" !!! author: Blaise Bourdin (2026, bourdin@mcmaster.ca) !!! !!! MEF90DefMechCellSetOptionsView: the default viewer for a MEF90DefMechCellSetOptions_Type !!! subroutine MEF90DefMechCellSetOptionsView(self, viewer, ierr) class(MEF90DefMechCellSetOptions_Type), 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') return write (IOBuffer, "(A,': Options for MEF90DefMechCellSetOptions_Type\n')") trim(self%prefix) PetscCall(PetscViewerASCIIPrintf(viewer, IOBuffer, ierr)) write (IOBuffer, "(' body force: ',3(ES12.5,' '),'crack pressure: ',ES12.5,'\n')") & self%bodyForce, self%crackPressure PetscCall(PetscViewerASCIIPrintf(viewer, IOBuffer, ierr)) write (IOBuffer, "(' plasticity type: ',A,'\n')") trim(MEF90DefMech_plasticityTypeList(self%plasticityType + 1)) PetscCall(PetscViewerASCIIPrintf(viewer, IOBuffer, ierr)) write (IOBuffer, "(' cohesive displacement: ',3(ES12.5,' '),'stiffness: ',ES12.5,'\n')") & self%cohesiveDisplacement, self%cohesiveStiffness PetscCall(PetscViewerASCIIPrintf(viewer, IOBuffer, ierr)) write (IOBuffer, "(' displacement BC: ',3(L1,' '),'value: ',3(ES12.5,' '),'\n')") & self%Has_displacementBC, self%boundaryDisplacement PetscCall(PetscViewerASCIIPrintf(viewer, IOBuffer, ierr)) write (IOBuffer, "(' displacement lower / upper bound: ',3(ES12.5,' '),' / ',3(ES12.5,' '),'\n')") & self%displacementLowerBound, self%displacementUpperBound PetscCall(PetscViewerASCIIPrintf(viewer, IOBuffer, ierr)) write (IOBuffer, "(' damage BC: ',L1,' value: ',ES12.5,'\n')") & self%Has_damageBC, self%boundaryDamage PetscCall(PetscViewerASCIIPrintf(viewer, IOBuffer, ierr)) write (IOBuffer, "(' crack volume controlled / work controlled: ',L1,' / ',L1,'\n')") & self%CrackVolumeControlled, self%WorkControlled PetscCall(PetscViewerASCIIPrintf(viewer, IOBuffer, ierr)) if (allocated(self%linearThermalExpansion)) then select type (alpha => self%linearThermalExpansion) type is (MatS2D) write (IOBuffer, "(' linear thermal expansion (alpha): ',2(ES12.5,', '),ES12.5,' [K^(-1)]\n')") alpha type is (MatS3D) write (IOBuffer, "(' linear thermal expansion (alpha): ',5(ES12.5,', '),ES12.5,' [K^(-1)]\n')") alpha class default write (IOBuffer, *) 'somehow linear thermal expansion is neither MatS2D nor MatS3D. This is wrong\n' end select PetscCall(PetscViewerASCIIPrintf(viewer, IOBuffer, ierr)) end if end subroutine MEF90DefMechCellSetOptionsView #undef __FUNCT__ #define __FUNCT__ "MEF90DefMechFaceSetOptionsSetFromOptions" !!! author: Blaise Bourdin (2026, bourdin@mcmaster.ca) !!! !!! MEF90DefMechFaceSetOptionsSetFromOptions: reads the options of a single face set !!! subroutine MEF90DefMechFaceSetOptionsSetFromOptions(comm, prefix, options, ierr) MPIU_Comm, intent(IN) :: comm character(len=*), intent(IN) :: prefix type(MEF90DefMechFaceSetOptions_Type), intent(OUT) :: options PetscErrorCode, intent(INOUT) :: ierr PetscInt :: nOpt PetscInt :: verbose = 0 options%comm = comm options%prefix = prefix options%name = trim(prefix)//"DefMechFaceSetOptions" PetscCall(PetscOptionsBegin(comm, prefix, "Options for a MEF90DefMech_Type face set", "mef90DefMech", ierr)) nOpt = 3 PetscCall(PetscOptionsRealArray('-boundaryForce', '[N.m^(-2) / N.m^(-1)] (f): boundary force', 'mef90DefMech', options%boundaryForce, nOpt, PETSC_NULL_BOOL, ierr)) PetscCall(PetscOptionsReal('-pressureForce', 'without unit: internal crack pressure', 'mef90DefMech', options%pressureForce, options%pressureForce, PETSC_NULL_BOOL, ierr)) nOpt = 3 PetscCall(PetscOptionsBoolArray('-DisplacementBC', 'Displacement has Dirichlet boundary Condition (Y/N)', 'mef90DefMech', options%Has_displacementBC, nOpt, PETSC_NULL_BOOL, ierr)) nOpt = 3 PetscCall(PetscOptionsRealArray('-boundaryDisplacement', '[m] (U): Displacement boundary value', 'mef90DefMech', options%boundaryDisplacement, nOpt, PETSC_NULL_BOOL, ierr)) nOpt = 3 PetscCall(PetscOptionsRealArray('-displacementLowerBound', '[m] (U): Displacement lower bound', 'mef90DefMech', options%displacementLowerBound, nOpt, PETSC_NULL_BOOL, ierr)) nOpt = 3 PetscCall(PetscOptionsRealArray('-displacementUpperBound', '[m] (U): Displacement upper bound', 'mef90DefMech', options%displacementUpperBound, nOpt, PETSC_NULL_BOOL, ierr)) PetscCall(PetscOptionsBool('-DamageBC', 'Damage has Dirichlet boundary Condition (Y/N)', 'mef90DefMech', options%Has_DamageBC, options%Has_DamageBC, PETSC_NULL_BOOL, ierr)) PetscCall(PetscOptionsReal('-boundaryDamage', '[unit-less] (alpha): Damage boundary value', 'mef90DefMech', options%boundaryDamage, options%boundaryDamage, 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 !! PETSC_VIEWER_STDOUT_WORLD and not PetscViewerASCIIGetStdout: the latter is collective, and the !! options are read where they are needed, which is not always on all ranks. view itself only !! prints the first time a given set is read, so it is collective only then. call options%view(PETSC_VIEWER_STDOUT_WORLD, ierr) end if end subroutine MEF90DefMechFaceSetOptionsSetFromOptions #undef __FUNCT__ #define __FUNCT__ "MEF90DefMechFaceSetOptionsView" !!! author: Blaise Bourdin (2026, bourdin@mcmaster.ca) !!! !!! MEF90DefMechFaceSetOptionsView: the default viewer for a MEF90DefMechFaceSetOptions_Type !!! subroutine MEF90DefMechFaceSetOptionsView(self, viewer, ierr) class(MEF90DefMechFaceSetOptions_Type), 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') return write (IOBuffer, "(A,': Options for MEF90DefMechFaceSetOptions_Type\n')") trim(self%prefix) PetscCall(PetscViewerASCIIPrintf(viewer, IOBuffer, ierr)) write (IOBuffer, "(' boundary force: ',3(ES12.5,' '),'pressure force: ',ES12.5,'\n')") & self%boundaryForce, self%pressureForce PetscCall(PetscViewerASCIIPrintf(viewer, IOBuffer, ierr)) write (IOBuffer, "(' displacement BC: ',3(L1,' '),'value: ',3(ES12.5,' '),'\n')") & self%Has_displacementBC, self%boundaryDisplacement PetscCall(PetscViewerASCIIPrintf(viewer, IOBuffer, ierr)) write (IOBuffer, "(' displacement lower / upper bound: ',3(ES12.5,' '),' / ',3(ES12.5,' '),'\n')") & self%displacementLowerBound, self%displacementUpperBound PetscCall(PetscViewerASCIIPrintf(viewer, IOBuffer, ierr)) write (IOBuffer, "(' damage BC: ',L1,' value: ',ES12.5,'\n')") & self%Has_damageBC, self%boundaryDamage PetscCall(PetscViewerASCIIPrintf(viewer, IOBuffer, ierr)) end subroutine MEF90DefMechFaceSetOptionsView #undef __FUNCT__ #define __FUNCT__ "MEF90DefMechVertexSetOptionsSetFromOptions" !!! author: Blaise Bourdin (2026, bourdin@mcmaster.ca) !!! !!! MEF90DefMechVertexSetOptionsSetFromOptions: reads the options of a single vertex set !!! subroutine MEF90DefMechVertexSetOptionsSetFromOptions(comm, prefix, options, ierr) MPIU_Comm, intent(IN) :: comm character(len=*), intent(IN) :: prefix type(MEF90DefMechVertexSetOptions_Type), intent(OUT) :: options PetscErrorCode, intent(INOUT) :: ierr PetscInt :: nOpt PetscInt :: verbose = 0 options%comm = comm options%prefix = prefix options%name = trim(prefix)//"DefMechVertexSetOptions" PetscCall(PetscOptionsBegin(comm, prefix, "Options for a MEF90DefMech_Type vertex set", "mef90DefMech", ierr)) nOpt = 3 PetscCall(PetscOptionsBoolArray('-DisplacementBC', 'Displacement has Dirichlet boundary Condition (Y/N)', 'mef90DefMech', options%Has_displacementBC, nOpt, PETSC_NULL_BOOL, ierr)) nOpt = 3 PetscCall(PetscOptionsRealArray('-boundaryDisplacement', '[m] (U): Displacement boundary value', 'mef90DefMech', options%boundaryDisplacement, nOpt, PETSC_NULL_BOOL, ierr)) nOpt = 3 PetscCall(PetscOptionsRealArray('-displacementLowerBound', '[m] (U): Displacement lower bound', 'mef90DefMech', options%displacementLowerBound, nOpt, PETSC_NULL_BOOL, ierr)) nOpt = 3 PetscCall(PetscOptionsRealArray('-displacementUpperBound', '[m] (U): Displacement upper bound', 'mef90DefMech', options%displacementUpperBound, nOpt, PETSC_NULL_BOOL, ierr)) PetscCall(PetscOptionsBool('-DamageBC', 'Damage has Dirichlet boundary Condition (Y/N)', 'mef90DefMech', options%Has_DamageBC, options%Has_DamageBC, PETSC_NULL_BOOL, ierr)) PetscCall(PetscOptionsReal('-boundaryDamage', '[unit-less] (alpha): boundaryDamage', 'mef90DefMech', options%boundaryDamage, options%boundaryDamage, 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 !! PETSC_VIEWER_STDOUT_WORLD and not PetscViewerASCIIGetStdout: the latter is collective, and the !! options are read where they are needed, which is not always on all ranks. view itself only !! prints the first time a given set is read, so it is collective only then. call options%view(PETSC_VIEWER_STDOUT_WORLD, ierr) end if end subroutine MEF90DefMechVertexSetOptionsSetFromOptions #undef __FUNCT__ #define __FUNCT__ "MEF90DefMechVertexSetOptionsView" !!! author: Blaise Bourdin (2026, bourdin@mcmaster.ca) !!! !!! MEF90DefMechVertexSetOptionsView: the default viewer for a MEF90DefMechVertexSetOptions_Type !!! subroutine MEF90DefMechVertexSetOptionsView(self, viewer, ierr) class(MEF90DefMechVertexSetOptions_Type), 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') return write (IOBuffer, "(A,': Options for MEF90DefMechVertexSetOptions_Type\n')") trim(self%prefix) PetscCall(PetscViewerASCIIPrintf(viewer, IOBuffer, ierr)) write (IOBuffer, "(' displacement BC: ',3(L1,' '),'value: ',3(ES12.5,' '),'\n')") & self%Has_displacementBC, self%boundaryDisplacement PetscCall(PetscViewerASCIIPrintf(viewer, IOBuffer, ierr)) write (IOBuffer, "(' displacement lower / upper bound: ',3(ES12.5,' '),' / ',3(ES12.5,' '),'\n')") & self%displacementLowerBound, self%displacementUpperBound PetscCall(PetscViewerASCIIPrintf(viewer, IOBuffer, ierr)) write (IOBuffer, "(' damage BC: ',L1,' value: ',ES12.5,'\n')") & self%Has_damageBC, self%boundaryDamage PetscCall(PetscViewerASCIIPrintf(viewer, IOBuffer, ierr)) end subroutine MEF90DefMechVertexSetOptionsView #undef __FUNCT__ #define __FUNCT__ "MEF90DefMechSetFromOptions" !!! author: Blaise Bourdin (2012-14, bourdin@lsu.edu) !!! author: Blaise Bourdin (2026, bourdin@mcmaster.ca) !!! !!! MEF90DefMechSetFromOptions: initializes a MEF90DefMech_Type from options !!! subroutine MEF90DefMechSetFromOptions(self, ierr) class(MEF90DefMech_Type), intent(INOUT) :: self PetscErrorCode, intent(INOUT) :: ierr type(tIS) :: setIS PetscInt, dimension(:), pointer :: setID PetscInt :: set character(len=MEF90MXSTRLEN) :: setPrefix type(MEF90DefMechGlobalOptions_Type) :: globalOptions type(MEF90DefMechCellSetOptions_Type) :: cellSetOptions type(MEF90DefMechFaceSetOptions_Type) :: faceSetOptions PetscBool :: printHelp !! !! Problem-wide options !! PetscCall(MEF90DefMechGlobalOptionsSetFromOptions(self%comm, trim(self%prefix), globalOptions, ierr)) !! !! The per-set options are read where they are needed, but they have to be visited once here so !! that they are registered with the options database, and so that the displacement bounds of all !! sets are known before the solvers are set up !! self%hasDisplacementBounds = PETSC_FALSE PetscCall(DMGetLabelIdIS(self%megaDM, MEF90CellSetLabelName, setIS, ierr)) PetscCall(MEF90ISAllGatherMerge(self%comm, setIS, ierr)) PetscCall(ISGetIndices(setIS, setID, ierr)) do set = 1, size(setID) write (setPrefix, "(A,'cs',I4.4,'_')") trim(self%prefix), setID(set) PetscCall(MEF90DefMechCellSetOptionsSetFromOptions(self%comm, trim(setPrefix), self%dim, cellSetOptions, ierr)) self%hasDisplacementBounds = any(cellSetOptions%displacementLowerBound /= MEF90NINFINITY) .or. self%hasDisplacementBounds self%hasDisplacementBounds = any(cellSetOptions%displacementUpperBound /= MEF90INFINITY) .or. self%hasDisplacementBounds end do PetscCall(ISRestoreIndices(setIS, setID, ierr)) PetscCall(ISDestroy(setIS, ierr)) PetscCall(DMGetLabelIdIS(self%megaDM, MEF90FaceSetLabelName, setIS, ierr)) PetscCall(MEF90ISAllGatherMerge(self%comm, setIS, ierr)) PetscCall(ISGetIndices(setIS, setID, ierr)) do set = 1, size(setID) write (setPrefix, "(A,'fs',I4.4,'_')") trim(self%prefix), setID(set) PetscCall(MEF90DefMechFaceSetOptionsSetFromOptions(self%comm, trim(setPrefix), faceSetOptions, ierr)) self%hasDisplacementBounds = any(faceSetOptions%displacementLowerBound /= MEF90NINFINITY) .or. self%hasDisplacementBounds self%hasDisplacementBounds = any(faceSetOptions%displacementUpperBound /= MEF90INFINITY) .or. self%hasDisplacementBounds end do PetscCall(ISRestoreIndices(setIS, setID, ierr)) PetscCall(ISDestroy(setIS, ierr)) PetscCall(PetscOptionsGetBool(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, "-verbose", printHelp, PETSC_NULL_BOOL, ierr)) if (printHelp) then call self%view(PETSC_VIEWER_STDOUT_WORLD, ierr) end if end subroutine MEF90DefMechSetFromOptions #undef __FUNCT__ #define __FUNCT__ "MEF90DefMechView" !!! author: Blaise Bourdin (2026, bourdin@mcmaster.ca) !!! !!! MEF90DefMechView: the default viewer for a MEF90DefMech_Type !!! subroutine MEF90DefMechView(self, viewer, ierr) class(MEF90DefMech_Type), intent(IN) :: self type(tPetscViewer), intent(IN) :: viewer PetscErrorCode, intent(INOUT) :: ierr character(len=MEF90MXSTRLEN, kind=c_char) :: viewerType type(tIS) :: setIS PetscInt, dimension(:), pointer :: setID PetscInt :: set character(len=MEF90MXSTRLEN) :: setPrefix type(MEF90DefMechGlobalOptions_Type) :: globalOptions type(MEF90DefMechCellSetOptions_Type) :: cellSetOptions type(MEF90DefMechFaceSetOptions_Type) :: faceSetOptions type(MEF90DefMechVertexSetOptions_Type) :: vertexSetOptions PetscCall(PetscViewerGetType(viewer, viewerType, ierr)) if (viewerType /= 'ascii') return !! The problem-wide and per-set options print themselves: each of the option types below has its !! own view, and MEF90Object%view makes sure that a given set is only ever printed once. PetscCall(MEF90DefMechGlobalOptionsSetFromOptions(self%comm, trim(self%prefix), globalOptions, ierr)) call globalOptions%view(viewer, ierr) PetscCall(DMGetLabelIdIS(self%megaDM, MEF90CellSetLabelName, setIS, ierr)) PetscCall(MEF90ISAllGatherMerge(self%comm, setIS, ierr)) PetscCall(ISGetIndices(setIS, setID, ierr)) do set = 1, size(setID) write (setPrefix, "(A,'cs',I4.4,'_')") trim(self%prefix), setID(set) PetscCall(MEF90DefMechCellSetOptionsSetFromOptions(self%comm, trim(setPrefix), self%dim, cellSetOptions, ierr)) call cellSetOptions%view(viewer, ierr) end do PetscCall(ISRestoreIndices(setIS, setID, ierr)) PetscCall(ISDestroy(setIS, ierr)) PetscCall(DMGetLabelIdIS(self%megaDM, MEF90FaceSetLabelName, setIS, ierr)) PetscCall(MEF90ISAllGatherMerge(self%comm, setIS, ierr)) PetscCall(ISGetIndices(setIS, setID, ierr)) do set = 1, size(setID) write (setPrefix, "(A,'fs',I4.4,'_')") trim(self%prefix), setID(set) PetscCall(MEF90DefMechFaceSetOptionsSetFromOptions(self%comm, trim(setPrefix), faceSetOptions, ierr)) call faceSetOptions%view(viewer, ierr) end do PetscCall(ISRestoreIndices(setIS, setID, ierr)) PetscCall(ISDestroy(setIS, ierr)) PetscCall(DMGetLabelIdIS(self%megaDM, MEF90VertexSetLabelName, setIS, ierr)) PetscCall(MEF90ISAllGatherMerge(self%comm, setIS, ierr)) PetscCall(ISGetIndices(setIS, setID, ierr)) do set = 1, size(setID) write (setPrefix, "(A,'vs',I4.4,'_')") trim(self%prefix), setID(set) PetscCall(MEF90DefMechVertexSetOptionsSetFromOptions(self%comm, trim(setPrefix), vertexSetOptions, ierr)) call vertexSetOptions%view(viewer, ierr) end do PetscCall(ISRestoreIndices(setIS, setID, ierr)) PetscCall(ISDestroy(setIS, ierr)) end subroutine MEF90DefMechView end module m_MEF90_DefMech_class