#include "../MEF90/mef90.inc" module MEF90_APPEND(m_MEF90_DefMechAssembly,MEF90_DIM)D #include "petsc/finclude/petsc.h" #define MEF90_HDRegularization 0.01_Kr use petscsnes use petsctao use m_MEF90_DefMech_class use m_MEF90_Materials use m_MEF90_HookesLaw use m_MEF90_DefMechSplit use m_MEF90_DefMechAT implicit none(type, external) ! Private public MEF90DefMechOperatorDisplacement, & MEF90DefMechBilinearFormDisplacement, & MEF90DefMechWork, & MEF90DefMechCohesiveEnergy, & MEF90DefMechPlasticDissipation, & MEF90DefMechElasticEnergy, & MEF90DefMechOperatorDamage, & MEF90DefMechTAOGradientDamage, & MEF90DefMechBilinearFormDamage, & MEF90DefMechTAOHessianDamage, & MEF90DefMechSurfaceEnergy, & MEF90DefMechTAOObjectiveDamage, & MEF90DefMechCrackVolume, & MEF90DefMechStress contains #undef __FUNCT__ #define __FUNCT__ "MEF90DefMechOperatorDisplacement" !!! author: Blaise Bourdin (2012-20, bourdin@lsu.edu) !!! author: Erwan Tanne (2012-20, erwan.tanne@gmail.com) !!! author: Blaise Bourdin (2022-26, bourdin@mcmaster.ca) !!! author: Alexis Marboeuf (2022, marboeua@mcmaster.ca) !!! !!! MEF90DefMechOperatorDisplacement: Build the operator. When called in SNES, the solution time should always match the target time, !!! so there is no need for interpolation of the forces, external, and boundary values !!! subroutine MEF90DefMechOperatorDisplacement(snesDisplacement, displacement, residual, MEF90DefMechCtx, ierr) type(tSNES), intent(IN) :: snesDisplacement type(tVec), intent(IN) :: displacement type(tVec), intent(INOUT) :: residual type(MEF90DefMech_Type), intent(IN) :: MEF90DefMechCtx PetscErrorCode, intent(INOUT) :: ierr type(tDM) :: dmDisplacement, dmDamage, dmTemperature, dmCohesiveDisplacement, dmBodyForce, dmBoundaryForce, dmPressureForce, dmPlasticStrain type(tPetscSection) :: sectionBodyForce, sectionBoundaryForce, sectionPressureForce, sectionPlasticStrain PetscReal, dimension(:), pointer :: cohesiveDisplacementDof, displacementDof, damageDof, temperatureDof PetscReal, dimension(:), pointer :: bodyForceArray, boundaryForceArray, pressureForceArray, plasticStrainArray type(MEF90_MATS) :: plasticStrainCell type(tIS) :: setIS, setPointIS PetscInt, dimension(:), pointer :: setID, setPointID PetscInt :: set, QuadratureOrder, cell, dim, vecOffset type(MEF90DefMechCellSetOptions_Type) :: cellSetOptions type(MEF90_MATS) :: linearThermalExpansion type(MEF90DefMechFaceSetOptions_Type) :: faceSetOptions type(MEF90_ELEMENT_ELAST), dimension(:), pointer :: elemVect type(MEF90_ELEMENT_SCAL), dimension(:), pointer :: elemScal type(eDMPolytopeType) :: cellGeometry type(MEF90ElementType) :: elemVectType, elemScalType PetscReal, dimension(:), pointer :: residualDof PetscReal :: residualDoFPlus, residualDoFMinus type(MEF90CtxGlobalOptions_Type) :: MEF90CtxGlobalOptions type(MEF90DefMechGlobalOptions_Type) :: MEF90DefMechGlobalOptions type(tVec) :: residualLocal class(MEF90DefMechAT_Type), allocatable :: ATModel class(MEF90DefMechSplit), allocatable :: Split type(MEF90_MATS) :: totalStrainGauss type(MEF90_VECT) :: U0Gauss, bodyForce, boundaryForce, pressureForce PetscReal :: damageGauss, temperatureGauss PetscInt :: iDof, iGauss, numDofDisplacement, numDofDamage, numGauss class(MEF90HookesLaw), allocatable :: HookesLaw character(len=MEF90MXSTRLEN) :: prefix type(tVec), pointer :: damageLocal PetscCall(MEF90CtxGlobalOptionsSetFromOptions(MEF90DefMechCtx%MEF90Ctx%comm, trim(MEF90DefMechCtx%MEF90Ctx%prefix), MEF90CtxGlobalOptions, ierr)) PetscCall(MEF90DefMechGlobalOptionsSetFromOptions(MEF90DefMechCtx%comm, trim(MEF90DefMechCtx%prefix), MEF90DefMechGlobalOptions, ierr)) PetscCall(VecGetDM(MEF90DefMechCtx%displacementLocal, dmDisplacement, ierr)) PetscCall(VecGetDM(MEF90DefMechCtx%cohesiveDisplacement, dmCohesiveDisplacement, ierr)) PetscCall(VecGetDM(MEF90DefMechCtx%damageLocal, dmDamage, ierr)) PetscCall(VecGetDM(MEF90DefMechCtx%TemperatureLocal, dmTemperature, ierr)) PetscCall(VecGetDM(MEF90DefMechCtx%bodyForce, dmBodyForce, ierr)) PetscCall(DMGetLocalSection(dmBodyForce, sectionBodyForce, ierr)) PetscCall(VecGetDM(MEF90DefMechCtx%plasticStrain, dmPlasticStrain, ierr)) PetscCall(DMGetLocalSection(dmPlasticStrain, sectionPlasticStrain, ierr)) PetscCall(VecGetDM(MEF90DefMechCtx%boundaryForce, dmBoundaryForce, ierr)) PetscCall(DMGetLocalSection(dmBoundaryForce, sectionBoundaryForce, ierr)) PetscCall(VecGetDM(MEF90DefMechCtx%pressureForce, dmPressureForce, ierr)) PetscCall(DMGetLocalSection(dmPressureForce, sectionPressureForce, ierr)) PetscCall(DMGetDimension(dmDisplacement, dim, ierr)) PetscCall(DMGetLocalVector(dmDisplacement, residualLocal, ierr)) PetscCall(DMGlobalToLocal(dmDisplacement, displacement, INSERT_VALUES, MEF90DefMechCtx%displacementLocal, ierr)) !! Something subtle is going on here: !! I _have_ to use MEF90DefMechCtx%displacementLocal because I need the constrained values, which would not be initialized if I were !! to create a new local Vec, or duplicate MEF90DefMechCtx%displacementLocal !! This DOES create a side effect. PetscCall(VecSet(residual, 0.0_kr, ierr)) PetscCall(VecSet(residualLocal, 0.0_kr, ierr)) !! get IS for cell sets PetscCall(DMGetLabelIdIS(dmDisplacement, MEF90CellSetLabelName, setIS, ierr)) PetscCall(MEF90ISAllGatherMerge(MEF90DefMechCtx%comm, setIS, ierr)) if (.not. PetscObjectIsNull(setIS)) then damageLocal => MEF90DefMechCtx%damageLocal PetscCall(VecGetArray(MEF90DefMechCtx%plasticStrain, plasticStrainArray, ierr)) PetscCall(VecGetArray(MEF90DefMechCtx%bodyForce, bodyForceArray, ierr)) PetscCall(ISGetIndices(setIS, setID, ierr)) do set = 1, size(setID) if (MEF90DefMechGlobalOptions%multiPhaseField) then damageLocal => MEF90DefMechCtx%partialDamageLocal(set) end if PetscCall(DMGetStratumIS(dmDisplacement, MEF90CellSetLabelName, setID(set), setPointIS, ierr)) if (.not. PetscObjectIsNull(setPointIS)) then write (prefix, '(A,"cs",I4.4,"_")') trim(MEF90DefMechCtx%prefix), setID(set) PetscCall(MEF90DefMechCellSetOptionsSetFromOptions(MEF90DefMechCtx%comm, prefix, MEF90_DIM, cellSetOptions, ierr)) select type (linearThermalExpansionMat => cellSetOptions%linearThermalExpansion) type is (MEF90_MATS) linearThermalExpansion = linearThermalExpansionMat end select PetscCall(ISGetIndices(setPointIS, setPointID, ierr)) PetscCall(DMPlexGetCellType(dmDisplacement, setPointID(1), cellGeometry, ierr)) PetscCall(MEF90ElementGetType(MEF90CtxGlobalOptions%elementFamily, MEF90CtxGlobalOptions%elementOrder, cellGeometry, elemVectType, ierr)) PetscCall(MEF90ElementGetType(MEF90CtxGlobalOptions%elementFamily, MEF90CtxGlobalOptions%elementOrder, cellGeometry, elemScalType, ierr)) !! get the ATModel and split objects PetscCall(MEF90DefMechGetATModel(MEF90DefMechCtx%comm, prefix, dim, ATModel, ierr)) PetscCall(ATModel%setFromOptions(ierr)) PetscCall(MEF90DefMechGetSplit(MEF90DefMechCtx%comm, prefix, Split, ierr)) PetscCall(Split%setFromOptions(ierr)) PetscCall(MEF90GetHookesLaw(MEF90DefMechCtx%comm, prefix, MEF90_DIM, HookesLaw, ierr)) PetscCall(HookesLaw%setFromOptions(ierr)) !! Allocate elements QuadratureOrder = max(2 * elemVectType%order, Split%quadratureOrder) PetscCall(MEF90ElementCreate(dmDisplacement, setPointIS, elemVect, QuadratureOrder, elemVectType, ierr)) PetscCall(MEF90ElementCreate(dmDamage, setPointIS, elemScal, QuadratureOrder, elemScalType, ierr)) numDofDisplacement = size(elemVect(1)%BF(:, 1)) numDofDamage = size(elemScal(1)%BF(:, 1)) numGauss = size(elemVect(1)%Gauss_C) allocate (residualDof(numDofDisplacement)) do cell = 1, size(setPointID) residualDof = 0.0_kr do iGauss = 1, numGauss !! Main term: [a(\alpha) sigma^+(u) + sigma^-(u)] . e(v) damageGauss = 0.0_kr if (.not. ATModel%isElastic) then PetscCall(DMPlexVecGetClosure(dmDamage, PETSC_NULL_SECTION, damageLocal, setPointID(cell), PETSC_NULL_INTEGER, damageDof, ierr)) do iDof = 1, numDofDamage damageGauss = damageGauss + damageDof(iDof) * elemScal(cell)%BF(iDof, iGauss) end do ! iDof numDofDamage PetscCall(DMPlexVecRestoreClosure(dmDamage, PETSC_NULL_SECTION, damageLocal, setPointID(cell), PETSC_NULL_INTEGER, damageDof, ierr)) end if totalStrainGauss = 0.0_kr PetscCall(DMPlexVecGetClosure(dmDisplacement, PETSC_NULL_SECTION, MEF90DefMechCtx%displacementLocal, setPointID(cell), PETSC_NULL_INTEGER, displacementDof, ierr)) do iDof = 1, numDofDisplacement totalStrainGauss = totalStrainGauss + displacementDof(iDof) * elemVect(cell)%GradS_BF(iDof, iGauss) end do ! iDof numDofDisplacement PetscCall(DMPlexVecRestoreClosure(dmDisplacement, PETSC_NULL_SECTION, MEF90DefMechCtx%displacementLocal, setPointID(cell), PETSC_NULL_INTEGER, displacementDof, ierr)) temperatureGauss = 0.0_kr PetscCall(DMPlexVecGetClosure(dmTemperature, PETSC_NULL_SECTION, MEF90DefMechCtx%TemperatureLocal, setPointID(cell), PETSC_NULL_INTEGER, temperatureDof, ierr)) do iDof = 1, numDofDamage temperatureGauss = temperatureGauss + temperatureDof(iDof) * elemScal(cell)%BF(iDof, iGauss) end do ! iDof numDofDamage PetscCall(DMPlexVecRestoreClosure(dmTemperature, PETSC_NULL_SECTION, MEF90DefMechCtx%TemperatureLocal, setPointID(cell), PETSC_NULL_INTEGER, temperatureDof, ierr)) totalStrainGauss = totalStrainGauss - (temperatureGauss * linearThermalExpansion) ! #if MEF90_DIM == 2 !! We need something along these lines !! Adding terms in planestrain for plasticity with tr(p) = 0 ! If (.NOT. matProp%HookesLaw%isPlaneStress) Then ! stressGauss = stressGauss + stiffness * ( matProp%HookesLaw%lambda*trace(plasticStrainCell)*MEF90MatS2DIdentity ) ! End If ! #endif PetscCall(PetscSectionGetOffset(sectionPlasticStrain, setPointID(cell), vecOffset, ierr)) plasticStrainCell = plasticStrainArray(vecOffset + 1:vecOffset + 1 + SIZEOFMEF90_MATS) call Split%setup(totalStrainGauss - plasticStrainCell, ierr) do iDof = 1, numDofDisplacement call Split%DEED(HookesLaw, elemVect(cell)%GradS_BF(iDof, iGauss), residualDoFPlus, residualDoFminus, ierr) if (ATModel%isElastic) then residualDof(iDof) = residualDof(iDof) + elemVect(cell)%Gauss_C(iGauss) * (residualDoFPlus + residualDoFminus) else if (split%isHybrid) then residualDof(iDof) = residualDof(iDof) + elemVect(cell)%Gauss_C(iGauss) * ATModel%a(damageGauss) * (residualDoFPlus + residualDoFminus) else residualDof(iDof) = residualDof(iDof) + elemVect(cell)%Gauss_C(iGauss) * (ATModel%a(damageGauss) * residualDoFPlus + residualDoFminus) end if ! isHybrid end if! isElastic end do ! iDof numDofDisplacement end do ! iGauss PetscCall(DMPlexVecSetClosure(dmDisplacement, PETSC_NULL_SECTION, residualLocal, setPointID(cell), residualDof, ADD_VALUES, ierr)) end do ! cell if (norm2(cellSetOptions%bodyForce) /= 0.0_kr) then do cell = 1, size(setPointID) residualDof = 0.0_kr !! This could break if TemperatureLocal had no dof in any point in the closure of setPointID(set) !! If this happens, we will need to protect this loop PetscCall(PetscSectionGetOffset(sectionBodyForce, setPointID(cell), vecOffset, ierr)) bodyForce = bodyForceArray(vecOffset + 1:vecOffset + 1 + SIZEOFMEF90_VECT) do iGauss = 1, numGauss do iDof = 1, numDofDisplacement residualDof(iDof) = residualDof(iDof) - elemVect(cell)%Gauss_C(iGauss) * (bodyForce .DotP. elemVect(cell)%BF(iDof, iGauss)) end do ! iDof numDofDisplacement end do ! iGauss PetscCall(DMPlexVecSetClosure(dmDisplacement, PETSC_NULL_SECTION, residualLocal, setPointID(cell), residualDof, ADD_VALUES, ierr)) end do ! cell end if ! cellSetOptions%bodyForce ! !!! Damping ! If (MEF90DefMechGlobalOptions%dampingCoefficientDisplacement * MEF90DefMechCtx%timeStep /= 0.0_Kr) Then ! Do cell = 1,size(setPointID) ! residualDof = 0.0_Kr ! Do iGauss = 1,numGauss ! displacementDampingGauss = 0.0_Kr ! PetscCall(DMPlexVecGetClosure(dmDisplacement,PETSC_NULL_SECTION,locDisplacement,setPointID(cell),PETSC_NULL_INTEGER,displacementDof,ierr)) ! PetscCall(DMPlexVecGetClosure(dmDisplacement,PETSC_NULL_SECTION,MEF90DefMechCtx%displacementPreviousStepLocal,setPointID(cell),PETSC_NULL_INTEGER,displacementPreviousStepDof,ierr)) ! Do iDof = 1,numDofDisplacement ! displacementDampingGauss = displacementDampingGauss + ((displacementDof(iDof) - displacementPreviousStepDof(iDof)) * elemVect(cell)%BF(iDof,iGauss)) ! End Do ! iDof numDofDisplacement ! PetscCall(DMPlexVecRestoreClosure(dmDisplacement,PETSC_NULL_SECTION,locDisplacement,setPointID(cell),PETSC_NULL_INTEGER,displacementDof,ierr)) ! PetscCall(DMPlexVecRestoreClosure(dmDisplacement,PETSC_NULL_SECTION,MEF90DefMechCtx%displacementPreviousStepLocal,setPointID(cell),PETSC_NULL_INTEGER,displacementPreviousStepDof,ierr)) ! displacementDampingGauss = displacementDampingGauss * MEF90DefMechGlobalOptions%dampingCoefficientDisplacement / MEF90DefMechCtx%timeStep ! Do iDof = 1,numDofDisplacement ! residualDof(iDof) = residualDof(iDof) + elemVect(cell)%Gauss_C(iGauss) * & ! ( displacementDampingGauss .DotP. elemVect(cell)%BF(iDof,iGauss) ) ! End Do ! iDof numDofDisplacement ! End Do ! iGauss ! PetscCall(DMPlexVecSetClosure(dmDisplacement,PETSC_NULL_SECTION,locResidual,setPointID(cell),residualDof,ADD_VALUES,ierr)) ! End Do ! cell ! End If ! damping !! Cohesive force if ((norm2(cellSetOptions%cohesiveDisplacement) /= 0.0_kr) .and. (cellSetOptions%cohesiveStiffness /= 0.0_kr)) then do cell = 1, size(setPointID) residualDof = 0.0_kr U0Gauss = 0.0_kr PetscCall(DMPlexVecGetClosure(dmDisplacement, PETSC_NULL_SECTION, MEF90DefMechCtx%displacementLocal, setPointID(cell), PETSC_NULL_INTEGER, displacementDof, ierr)) PetscCall(DMPlexVecGetClosure(dmCohesiveDisplacement, PETSC_NULL_SECTION, MEF90DefMechCtx%cohesiveDisplacement, setPointID(cell), PETSC_NULL_INTEGER, cohesiveDisplacementDof, ierr)) do iGauss = 1, numGauss do iDof = 1, numDofDisplacement U0Gauss = U0Gauss + (displacementDof(iDof) - cohesiveDisplacementDof(iDof)) * elemVect(cell)%BF(iDof, iGauss) end do ! iDof numDofDisplacement U0Gauss = U0Gauss * cellSetOptions%cohesiveStiffness do iDof = 1, numDofDisplacement residualDof(iDof) = residualDof(iDof) + elemVect(cell)%Gauss_C(iGauss) * (U0Gauss .DotP. elemVect(cell)%BF(iDof, iGauss)) end do ! iDof numDofDisplacement end do ! iGauss PetscCall(DMPlexVecRestoreClosure(dmDisplacement, PETSC_NULL_SECTION, MEF90DefMechCtx%displacementLocal, setPointID(cell), PETSC_NULL_INTEGER, displacementDof, ierr)) PetscCall(DMPlexVecRestoreClosure(dmCohesiveDisplacement, PETSC_NULL_SECTION, MEF90DefMechCtx%cohesiveDisplacement, setPointID(cell), PETSC_NULL_INTEGER, cohesiveDisplacementDof, ierr)) PetscCall(DMPlexVecSetClosure(dmDisplacement, PETSC_NULL_SECTION, residualLocal, setPointID(cell), residualDof, ADD_VALUES, ierr)) end do ! cell end if ! cohesiveDisplacement ! !!! crack pressure ! If (cellSetOptions%crackPressure /= 0.0_Kr)) Then ! Do cell = 1,size(setPointID) ! residualDof = 0.0_Kr ! Do iGauss = 1,numGauss ! gradientDamageGauss = 0.0_Kr ! PetscCall(DMPlexVecGetClosure(dmDamage,PETSC_NULL_SECTION,MEF90DefMechCtx%damageLocal(1),setPointID(cell),PETSC_NULL_INTEGER,damageDof,ierr)) ! Do iDof = 1,numDofDamage ! GradientDamageGauss = GradientDamageGauss + damageDof(iDof) * elemScal(cell)%Grad_BF(iDof,iGauss) ! End Do ! iDof numDofDamage ! PetscCall(DMPlexVecRestoreClosure(dmDamage,PETSC_NULL_SECTION,MEF90DefMechCtx%damageLocal(1),setPointID(cell),PETSC_NULL_INTEGER,damageDof,ierr)) ! CrackPressureCell = 0.0_Kr ! PetscCall(DMPlexVecGetClosure(dmCrackPressure,PETSC_NULL_SECTION,MEF90DefMechCtx%crackPressure,setPointID(cell),PETSC_NULL_INTEGER,crackPressureDof,ierr)) ! Do iDof = 1,numDofDamage ! CrackPressureCell = CrackPressureCell + crackPressureDof(iDof) * elemScal(cell)%Grad_BF(iDof,iGauss) ! End Do ! iDof numDofDamage ! PetscCall(DMPlexVecRestoreClosure(dmCrackPressure,PETSC_NULL_SECTION,MEF90DefMechCtx%crackPressure,setPointID(cell),PETSC_NULL_INTEGER,crackPressureDof,ierr)) ! Do iDof = 1,numDofDisplacement ! residualDof(iDof) = residualDof(iDof) + elemVect(cell)%Gauss_C(iGauss) * & ! CrackPressureCell * (GradientDamageGauss .DotP. elemVect(cell)%BF(iDof,iGauss)) ! End Do ! iDof numDofDisplacement ! End Do ! iGauss ! PetscCall(DMPlexVecSetClosure(dmDisplacement,PETSC_NULL_SECTION,locResidual,setPointID(cell),residualDof,ADD_VALUES,ierr)) ! End Do ! cell ! End If ! crack Pressure deallocate (residualDof) PetscCall(ISRestoreIndices(setPointIS, setPointID, ierr)) PetscCall(MEF90ElementDestroy(elemVect, ierr)) PetscCall(MEF90ElementDestroy(elemScal, ierr)) PetscCall(ISDestroy(setPointIS, ierr)) end if ! setPointIS end do ! set PetscCall(ISRestoreIndices(setIS, setID, ierr)) PetscCall(VecRestoreArray(MEF90DefMechCtx%plasticStrain, plasticStrainArray, ierr)) PetscCall(VecRestoreArray(MEF90DefMechCtx%bodyForce, bodyForceArray, ierr)) PetscCall(ISDestroy(setIS, ierr)) end if ! setIS !! face-based contributions PetscCall(DMGetLabelIdIS(dmDisplacement, MEF90FaceSetLabelName, setIS, ierr)) PetscCall(MEF90ISAllGatherMerge(MEF90DefMechCtx%comm, setIS, ierr)) if (.not. PetscObjectIsNull(setIS)) then PetscCall(ISGetIndices(setIS, setID, ierr)) PetscCall(VecGetArray(MEF90DefMechCtx%boundaryForce, boundaryForceArray, ierr)) PetscCall(VecGetArray(MEF90DefMechCtx%pressureForce, pressureForceArray, ierr)) do set = 1, size(setID) PetscCall(DMGetStratumIS(dmDisplacement, MEF90FaceSetLabelName, setID(set), setPointIS, ierr)) if (.not. PetscObjectIsNull(setPointIS)) then write (prefix, '(A,"fs",I4.4,"_")') trim(MEF90DefMechCtx%prefix), setID(set) PetscCall(MEF90DefMechFaceSetOptionsSetFromOptions(MEF90DefMechCtx%comm, prefix, faceSetOptions, ierr)) PetscCall(ISGetIndices(setPointIS, setPointID, ierr)) PetscCall(DMPlexGetCellType(dmDisplacement, setPointID(1), cellGeometry, ierr)) PetscCall(MEF90ElementGetTypeBoundary(MEF90CtxGlobalOptions%elementFamily, MEF90CtxGlobalOptions%elementOrder, cellGeometry, elemVectType, ierr)) PetscCall(MEF90ElementGetTypeBoundary(MEF90CtxGlobalOptions%elementFamily, MEF90CtxGlobalOptions%elementOrder, cellGeometry, elemScalType, ierr)) !! Allocate elements QuadratureOrder = 2 * elemVectType%order PetscCall(MEF90ElementCreate(dmDisplacement, setPointIS, elemVect, QuadratureOrder, elemVectType, ierr)) PetscCall(MEF90ElementCreate(dmDamage, setPointIS, elemScal, QuadratureOrder, elemScalType, ierr)) numDofDisplacement = size(elemVect(1)%BF(:, 1)) numDofDamage = size(elemScal(1)%BF(:, 1)) numGauss = size(elemVect(1)%Gauss_C) allocate (residualDof(numDofDisplacement)) if (norm2(faceSetOptions%boundaryForce) /= 0.0_kr) then do cell = 1, size(setPointID) residualDof = 0.0_kr PetscCall(PetscSectionGetOffset(sectionBoundaryForce, setPointID(cell), vecOffset, ierr)) boundaryForce = boundaryForceArray(vecOffset + 1:vecOffset + 1 + SIZEOFMEF90_VECT) do iGauss = 1, numGauss do iDof = 1, numDofDisplacement residualDof(iDof) = residualDof(iDof) - elemVect(cell)%Gauss_C(iGauss) * (boundaryForce .DotP. elemVect(cell)%BF(iDof, iGauss)) end do ! iDof numDofDisplacement end do ! iGauss PetscCall(DMPlexVecSetClosure(dmDisplacement, PETSC_NULL_SECTION, residualLocal, setPointID(cell), residualDof, ADD_VALUES, ierr)) end do ! cell end if ! faceSetOptions%boundaryForce if (faceSetOptions%pressureForce /= 0.0_kr) then do cell = 1, size(setPointID) residualDof = 0.0_kr PetscCall(PetscSectionGetOffset(sectionPressureForce, setPointID(cell), vecOffset, ierr)) pressureForce = pressureForceArray(vecOffset + 1) * elemVect(cell)%outerNormal do iGauss = 1, numGauss do iDof = 1, numDofDisplacement residualDof(iDof) = residualDof(iDof) - elemVect(cell)%Gauss_C(iGauss) * (pressureForce .DotP. elemVect(cell)%BF(iDof, iGauss)) end do ! iDof numDofDisplacement end do ! iGauss PetscCall(DMPlexVecSetClosure(dmDisplacement, PETSC_NULL_SECTION, residualLocal, setPointID(cell), residualDof, ADD_VALUES, ierr)) end do ! cell end if ! faceSetOptions%boundaryForce deallocate (residualDof) PetscCall(MEF90ElementDestroy(elemVect, ierr)) PetscCall(MEF90ElementDestroy(elemScal, ierr)) PetscCall(ISRestoreIndices(setPointIS, setPointID, ierr)) PetscCall(ISDestroy(setPointIS, ierr)) end if ! setPointIS end do ! set PetscCall(ISRestoreIndices(setIS, setID, ierr)) PetscCall(VecRestoreArray(MEF90DefMechCtx%boundaryForce, boundaryForceArray, ierr)) PetscCall(VecRestoreArray(MEF90DefMechCtx%pressureForce, pressureForceArray, ierr)) PetscCall(ISDestroy(setIS, ierr)) end if ! setIS PetscCall(DMLocalToGlobalBegin(dmDisplacement, residualLocal, ADD_VALUES, residual, ierr)) PetscCall(DMLocalToGlobalEnd(dmDisplacement, residualLocal, ADD_VALUES, residual, ierr)) PetscCall(DMRestoreLocalVector(dmDisplacement, residualLocal, ierr)) end subroutine MEF90DefMechOperatorDisplacement #undef __FUNCT__ #define __FUNCT__ "MEF90DefMechBilinearFormDisplacement" !!! author: Blaise Bourdin (2012-20, bourdin@lsu.edu) !!! author: Erwan Tanne (2012-20, erwan.tanne@gmail.com) !!! author: Blaise Bourdin (2022-26, bourdin@mcmaster.ca) !!! author: Alexis Marboeuf (2022, marboeua@mcmaster.ca) !!! !!! MEF90DefMechBilinearFormDisplacement: !!! subroutine MEF90DefMechBilinearFormDisplacement(snesDisplacement, displacement, A, M, MEF90DefMechCtx, ierr) type(tSNES), intent(IN) :: snesDisplacement type(tVec), intent(IN) :: displacement type(tMat), intent(INOUT) :: A, M type(MEF90DefMech_Type), intent(IN) :: MEF90DefMechCtx PetscErrorCode, intent(INOUT) :: ierr type(tDM) :: dmDisplacement, dmDamage, dmTemperature, dmPlasticStrain type(tPetscSection) :: sectionPlasticStrain type(tIS) :: setIS, setPointIS PetscInt :: numDofDisplacement, numDofDamage, numGauss, set, cell, iGauss, iDof, jDof, QuadratureOrder, dim, vecOffset PetscReal :: damageGauss, temperatureGauss PetscInt, dimension(:), pointer :: setID, setPointID PetscReal, dimension(:,:), pointer :: matDof PetscReal, dimension(:), pointer :: displacementDof, damageDof, temperatureDof, plasticStrainArray type(MEF90_MATS) :: totalStrainGauss, plasticStrainCell type(MEF90_VECT) :: U0Gauss type(MEF90DefMechCellSetOptions_Type) :: cellSetOptions type(MEF90_MATS) :: linearThermalExpansion type(MEF90_ELEMENT_ELAST), dimension(:), pointer :: elemVect type(MEF90_ELEMENT_SCAL), dimension(:), pointer :: elemScal type(eDMPolytopeType) :: cellGeometry type(MEF90ElementType) :: elemVectType, elemScalType type(MEF90CtxGlobalOptions_Type) :: MEF90CtxGlobalOptions type(MEF90DefMechGlobalOptions_Type) :: MEF90DefMechGlobalOptions class(MEF90DefMechAT_Type), allocatable :: ATModel class(MEF90DefMechSplit), allocatable :: Split character(len=MEF90MXSTRLEN) :: prefix class(MEF90HookesLaw), allocatable :: HookesLaw PetscReal :: D2EEDPlus, D2EEDMinus type(tVec), pointer :: damageLocal PetscCall(MEF90CtxGlobalOptionsSetFromOptions(MEF90DefMechCtx%MEF90Ctx%comm, trim(MEF90DefMechCtx%MEF90Ctx%prefix), MEF90CtxGlobalOptions, ierr)) PetscCall(MEF90DefMechGlobalOptionsSetFromOptions(MEF90DefMechCtx%comm, trim(MEF90DefMechCtx%prefix), MEF90DefMechGlobalOptions, ierr)) PetscCall(VecGetDM(MEF90DefMechCtx%displacementLocal, dmDisplacement, ierr)) PetscCall(VecGetDM(MEF90DefMechCtx%damageLocal, dmDamage, ierr)) PetscCall(VecGetDM(MEF90DefMechCtx%temperatureLocal, dmTemperature, ierr)) PetscCall(VecGetDM(MEF90DefMechCtx%plasticStrain, dmPlasticStrain, ierr)) PetscCall(DMGetLocalSection(dmPlasticStrain, sectionPlasticStrain, ierr)) PetscCall(DMGetDimension(dmDisplacement, dim, ierr)) PetscCall(DMGlobalToLocal(dmDisplacement, displacement, INSERT_VALUES, MEF90DefMechCtx%displacementLocal, ierr)) PetscCall(MatZeroEntries(A, ierr)) !! get IS for cell sets PetscCall(DMGetLabelIdIS(dmDisplacement, MEF90CellSetLabelName, setIS, ierr)) PetscCall(MEF90ISAllGatherMerge(MEF90DefMechCtx%comm, setIS, ierr)) if (.not. PetscObjectIsNull(setIS)) then damageLocal => MEF90DefMechCtx%damageLocal PetscCall(ISGetIndices(setIS, setID, ierr)) PetscCall(VecGetArray(MEF90DefMechCtx%plasticStrain, plasticStrainArray, ierr)) do set = 1, size(setID) if (MEF90DefMechGlobalOptions%multiPhaseField) then damageLocal => MEF90DefMechCtx%partialDamageLocal(set) end if PetscCall(DMGetStratumIS(dmDisplacement, MEF90CellSetLabelName, setID(set), setPointIS, ierr)) if (.not. PetscObjectIsNull(setPointIS)) then write (prefix, '(A,"cs",I4.4,"_")') trim(MEF90DefMechCtx%prefix), setID(set) PetscCall(MEF90DefMechCellSetOptionsSetFromOptions(MEF90DefMechCtx%comm, prefix, MEF90_DIM, cellSetOptions, ierr)) select type (linearThermalExpansionMat => cellSetOptions%linearThermalExpansion) type is (MEF90_MATS) linearThermalExpansion = linearThermalExpansionMat end select PetscCall(ISGetIndices(setPointIS, setPointID, ierr)) PetscCall(DMPlexGetCellType(dmDisplacement, setPointID(1), cellGeometry, ierr)) PetscCall(MEF90ElementGetType(MEF90CtxGlobalOptions%elementFamily, MEF90CtxGlobalOptions%elementOrder, cellGeometry, elemVectType, ierr)) PetscCall(MEF90ElementGetType(MEF90CtxGlobalOptions%elementFamily, MEF90CtxGlobalOptions%elementOrder, cellGeometry, elemScalType, ierr)) !! get the ATModel and split objects PetscCall(MEF90DefMechGetATModel(MEF90DefMechCtx%comm, prefix, dim, ATModel, ierr)) PetscCall(ATModel%setFromOptions(ierr)) PetscCall(MEF90DefMechGetSplit(MEF90DefMechCtx%comm, prefix, Split, ierr)) PetscCall(Split%setFromOptions(ierr)) PetscCall(MEF90GetHookesLaw(MEF90DefMechCtx%Comm, prefix, MEF90_DIM, HookesLaw, ierr)) PetscCall(HookesLaw%setFromOptions(ierr)) !! Allocate elements QuadratureOrder = max(2 * elemVectType%order, Split%quadratureOrder) PetscCall(MEF90ElementCreate(dmDisplacement, setPointIS, elemVect, QuadratureOrder, elemVectType, ierr)) PetscCall(MEF90ElementCreate(dmDamage, setPointIS, elemScal, QuadratureOrder, elemScalType, ierr)) numDofDisplacement = size(elemVect(1)%BF(:, 1)) numDofDamage = size(elemScal(1)%BF(:, 1)) numGauss = size(elemVect(1)%Gauss_C) allocate (matDof(numDofDisplacement, numDofDisplacement)) do cell = 1, size(setPointID) matDof = 0.0_kr do iGauss = 1, numGauss !! Main term: [a(\alpha) sigma^+(u) + sigma^-(u)] . e(v) damageGauss = 0.0_kr if (.not. ATModel%isElastic) then PetscCall(DMPlexVecGetClosure(dmDamage, PETSC_NULL_SECTION, damageLocal, setPointID(cell), PETSC_NULL_INTEGER, damageDof, ierr)) do iDof = 1, numDofDamage damageGauss = damageGauss + damageDof(iDof) * elemScal(cell)%BF(iDof, iGauss) end do ! iDof numDofDamage PetscCall(DMPlexVecRestoreClosure(dmDamage, PETSC_NULL_SECTION, damageLocal, setPointID(cell), PETSC_NULL_INTEGER, damageDof, ierr)) end if totalStrainGauss = 0.0_kr PetscCall(DMPlexVecGetClosure(dmDisplacement, PETSC_NULL_SECTION, MEF90DefMechCtx%displacementLocal, setPointID(cell), PETSC_NULL_INTEGER, displacementDof, ierr)) do iDof = 1, numDofDisplacement totalStrainGauss = totalStrainGauss + displacementDof(iDof) * elemVect(cell)%GradS_BF(iDof, iGauss) end do ! iDof numDofDisplacement PetscCall(DMPlexVecRestoreClosure(dmDisplacement, PETSC_NULL_SECTION, MEF90DefMechCtx%displacementLocal, setPointID(cell), PETSC_NULL_INTEGER, displacementDof, ierr)) temperatureGauss = 0.0_kr PetscCall(DMPlexVecGetClosure(dmTemperature, PETSC_NULL_SECTION, MEF90DefMechCtx%temperatureLocal, setPointID(cell), PETSC_NULL_INTEGER, temperatureDof, ierr)) do iDof = 1, numDofDamage temperatureGauss = temperatureGauss + temperatureDof(iDof) * elemScal(cell)%BF(iDof, iGauss) end do ! iDof numDofDamage PetscCall(DMPlexVecRestoreClosure(dmTemperature, PETSC_NULL_SECTION, MEF90DefMechCtx%temperatureLocal, setPointID(cell), PETSC_NULL_INTEGER, temperatureDof, ierr)) totalStrainGauss = totalStrainGauss - (temperatureGauss * linearThermalExpansion) ! #if MEF90_DIM == 2 !! We need something along these lines !! Adding terms in planestrain for plasticity with tr(p) = 0 ! If (.NOT. matProp%HookesLaw%isPlaneStress) Then ! stressGauss = stressGauss + stiffness * ( matProp%HookesLaw%lambda*trace(plasticStrainCell)*MEF90MatS2DIdentity ) ! End If ! #endif PetscCall(PetscSectionGetOffset(sectionPlasticStrain, setPointID(cell), vecOffset, ierr)) plasticStrainCell = plasticStrainArray(vecOffset + 1:vecOffset + 1 + SIZEOFMEF90_MATS) call Split%setup(totalStrainGauss - plasticStrainCell, ierr) do iDof = 1, numDofDisplacement do jDof = 1, numDofDisplacement call Split%D2EED(HookesLaw, elemVect(cell)%GradS_BF(iDof, iGauss), elemVect(cell)%GradS_BF(jDof, iGauss), D2EEDPlus, D2EEDMinus, ierr) if (ATModel%isElastic) then matDof(jDof, iDof) = matDof(jDof, iDof) + elemVect(cell)%Gauss_C(iGauss) * (D2EEDPlus + D2EEDMinus) else if (split%isHybrid) then matDof(jDof, iDof) = matDof(jDof, iDof) + elemVect(cell)%Gauss_C(iGauss) * ATModel%a(damageGauss) * (D2EEDPlus + D2EEDMinus) else matDof(jDof, iDof) = matDof(jDof, iDof) + elemVect(cell)%Gauss_C(iGauss) * (ATModel%a(damageGauss) * D2EEDPlus + D2EEDMinus) end if ! isHybrid end if ! isElastic end do ! jDof numDofDisplacement end do ! iDof numDofDisplacement end do ! iGauss PetscCall(DMPlexMatSetClosure(dmDisplacement, PETSC_NULL_SECTION, PETSC_NULL_SECTION, A, setPointID(cell), reshape(matDof, [numDofDisplacement*numDofDisplacement]), ADD_VALUES, ierr)) end do ! cell ! !!! Damping ! If (MEF90DefMechGlobalOptions%dampingCoefficientDisplacement * MEF90DefMechCtx%timeStep /= 0.0_Kr) Then ! Do cell = 1,size(setPointID) ! matDof = 0.0_Kr ! Do iGauss = 1,numGauss ! Do jDof = 0,numDofDisplacement-1 ! Do iDof = 1,numDofDisplacement ! U0Gauss = elemVect(cell)%BF(iDof,iGauss) * MEF90DefMechGlobalOptions%dampingCoefficientDisplacement / MEF90DefMechCtx%timeStep ! matDof(jDof*numDofDisplacement+iDof) = matDof(jDof*numDofDisplacement+iDof) + elemVect(cell)%Gauss_C(iGauss) * (U0Gauss .DotP. elemVect(cell)%BF(jDof+1,iGauss)) ! End Do ! iDof numDofDisplacement ! End Do ! jDof numDofDisplacement ! End Do ! iGauss ! PetscCall(DMPlexMatSetClosure(dmDisplacement,PETSC_NULL_SECTION,PETSC_NULL_SECTION,A,setPointID(cell),matDof,ADD_VALUES,ierr)) ! End Do ! cell ! End If ! damping !! Cohesive force if ((norm2(cellSetOptions%cohesiveDisplacement) /= 0.0_kr) .and. (cellSetOptions%cohesiveStiffness /= 0.0_kr)) then do cell = 1, size(setPointID) matDof = 0.0_kr do iGauss = 1, numGauss do iDof = 1, numDofDisplacement U0Gauss = cellSetOptions%cohesiveStiffness * elemVect(cell)%BF(iDof, iGauss) do jDof = 1, numDofDisplacement matDof(jDof, iDof) = matDof(jDof, iDof) + elemVect(cell)%Gauss_C(iGauss) * (U0Gauss .DotP. elemVect(cell)%BF(jDof, iGauss)) end do ! jDof numDofDisplacement end do ! iDof numDofDisplacement end do ! iGauss PetscCall(DMPlexMatSetClosure(dmDisplacement, PETSC_NULL_SECTION, PETSC_NULL_SECTION, A, setPointID(cell), reshape(matDof, [numDofDisplacement*numDofDisplacement]), ADD_VALUES, ierr)) end do ! cell end if ! cohesiveDisplacement deallocate (matDof) PetscCall(MEF90ElementDestroy(elemVect, ierr)) PetscCall(MEF90ElementDestroy(elemScal, ierr)) PetscCall(ISRestoreIndices(setPointIS, setPointID, ierr)) PetscCall(ISDestroy(setPointIS, ierr)) end if ! setPointIS end do ! set PetscCall(ISRestoreIndices(setIS, setID, ierr)) PetscCall(VecRestoreArray(MEF90DefMechCtx%plasticStrain, plasticStrainArray, ierr)) PetscCall(ISDestroy(setIS, ierr)) end if ! setIS PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr)) PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr)) PetscCall(MatCopy(A, M, SAME_NONZERO_PATTERN, ierr)) end subroutine MEF90DefMechBilinearFormDisplacement #undef __FUNCT__ #define __FUNCT__ "MEF90DefMechWork" !!! author: Blaise Bourdin (2014-2022, bourdin@lsu.edu) !!! author: Blaise Bourdin (2022-26, bourdin@mcmaster.ca) !!! author: Alexis Marboeuf (2022, marboeua@mcmaster.ca) !!! !!! MEF90DefMechWork: !!! subroutine MEF90DefMechWork(MEF90DefMechCtx, bodyForceWork, boundaryForceWork, ierr) type(MEF90DefMech_Type), intent(IN) :: MEF90DefMechCtx PetscReal, dimension(:), pointer :: bodyForceWork, boundaryForceWork PetscErrorCode, intent(INOUT) :: ierr type(tDM) :: dmDisplacement, dmBodyForce, dmBoundaryForce, dmPressureForce type(tPetscSection) :: sectionBodyForce, sectionBoundaryForce, sectionPressureForce PetscReal, dimension(:), pointer :: displacementDof, bodyForceArray, boundaryForceArray, pressureForceArray type(tIS) :: setIS, setPointIS PetscInt, dimension(:), pointer :: setID, setPointID PetscInt :: set, QuadratureOrder, cell, dim, vecOffset type(MEF90DefMechCellSetOptions_Type) :: cellSetOptions type(MEF90DefMechFaceSetOptions_Type) :: faceSetOptions type(MEF90_ELEMENT_ELAST), dimension(:), pointer :: elemVect type(eDMPolytopeType) :: cellGeometry type(MEF90ElementType) :: elemVectType character(len=MEF90MXSTRLEN) :: prefix type(MEF90CtxGlobalOptions_Type) :: MEF90CtxGlobalOptions type(MEF90DefMechGlobalOptions_Type) :: MEF90DefMechGlobalOptions type(MEF90_VECT) :: bodyForce, boundaryForce, pressureForce PetscReal :: myWork PetscInt :: iDof, iGauss, numDofDisplacement, numGauss PetscCall(MEF90CtxGlobalOptionsSetFromOptions(MEF90DefMechCtx%MEF90Ctx%comm, trim(MEF90DefMechCtx%MEF90Ctx%prefix), MEF90CtxGlobalOptions, ierr)) PetscCall(MEF90DefMechGlobalOptionsSetFromOptions(MEF90DefMechCtx%comm, trim(MEF90DefMechCtx%prefix), MEF90DefMechGlobalOptions, ierr)) PetscCall(VecGetDM(MEF90DefMechCtx%displacementLocal, dmDisplacement, ierr)) PetscCall(VecGetDM(MEF90DefMechCtx%bodyForce, dmBodyForce, ierr)) PetscCall(DMGetLocalSection(dmBodyForce, sectionBodyForce, ierr)) PetscCall(VecGetDM(MEF90DefMechCtx%boundaryForce, dmBoundaryForce, ierr)) PetscCall(DMGetLocalSection(dmBoundaryForce, sectionBoundaryForce, ierr)) PetscCall(VecGetDM(MEF90DefMechCtx%pressureForce, dmPressureForce, ierr)) PetscCall(DMGetLocalSection(dmPressureForce, sectionPressureForce, ierr)) PetscCall(DMGetDimension(dmDisplacement, dim, ierr)) !! get IS for cell sets PetscCall(DMGetLabelIdIS(dmDisplacement, MEF90CellSetLabelName, setIS, ierr)) PetscCall(MEF90ISAllGatherMerge(MEF90DefMechCtx%comm, setIS, ierr)) bodyForceWork = 0.0_kr if (.not. PetscObjectIsNull(setIS)) then PetscCall(ISGetIndices(setIS, setID, ierr)) PetscCall(VecGetArray(MEF90DefMechCtx%bodyForce, bodyForceArray, ierr)) do set = 1, size(setID) myWork = 0.0_kr PetscCall(DMGetStratumIS(dmDisplacement, MEF90CellSetLabelName, setID(set), setPointIS, ierr)) if (.not. PetscObjectIsNull(setPointIS)) then write (prefix, '(A,"cs",I4.4,"_")') trim(MEF90DefMechCtx%prefix), setID(set) PetscCall(MEF90DefMechCellSetOptionsSetFromOptions(MEF90DefMechCtx%comm, prefix, MEF90_DIM, cellSetOptions, ierr)) PetscCall(ISGetIndices(setPointIS, setPointID, ierr)) PetscCall(DMPlexGetCellType(dmDisplacement, setPointID(1), cellGeometry, ierr)) PetscCall(MEF90ElementGetType(MEF90CtxGlobalOptions%elementFamily, MEF90CtxGlobalOptions%elementOrder, cellGeometry, elemVectType, ierr)) !! Allocate elements QuadratureOrder = 2 * (elemVectType%order - 1) PetscCall(MEF90ElementCreate(dmDisplacement, setPointIS, elemVect, QuadratureOrder, elemVectType, ierr)) numDofDisplacement = size(elemVect(1)%BF(:, 1)) numGauss = size(elemVect(1)%Gauss_C) if (norm2(cellSetOptions%bodyForce) /= 0.0_kr) then do cell = 1, size(setPointID) PetscCall(PetscSectionGetOffset(sectionBodyForce, setPointID(cell), vecOffset, ierr)) bodyForce = bodyForceArray(vecOffset + 1:vecOffset + 1 + SIZEOFMEF90_VECT) PetscCall(DMPlexVecGetClosure(dmDisplacement, PETSC_NULL_SECTION, MEF90DefMechCtx%displacementLocal, setPointID(cell), PETSC_NULL_INTEGER, displacementDof, ierr)) do iGauss = 1, numGauss do iDof = 1, numDofDisplacement myWork = myWork + elemVect(cell)%Gauss_C(iGauss) * displacementDof(iDof) * (bodyForce .DotP. elemVect(cell)%BF(iDof, iGauss)) end do ! iDof numDofDisplacement end do ! iGauss PetscCall(DMPlexVecRestoreClosure(dmDisplacement, PETSC_NULL_SECTION, MEF90DefMechCtx%displacementLocal, setPointID(cell), PETSC_NULL_INTEGER, displacementDof, ierr)) end do ! cell end if ! cellSetOptions%bodyForce PetscCall(MEF90ElementDestroy(elemVect, ierr)) PetscCall(ISRestoreIndices(setPointIS, setPointID, ierr)) PetscCall(ISDestroy(setPointIS, ierr)) end if ! setPointIS PetscCallMPI(MPI_AllReduce(myWork, bodyForceWork(set), 1, MPIU_SCALAR, MPI_SUM, MEF90DefMechCtx%comm, ierr)) end do ! set PetscCall(ISRestoreIndices(setIS, setID, ierr)) PetscCall(VecRestoreArray(MEF90DefMechCtx%bodyForce, bodyForceArray, ierr)) PetscCall(ISDestroy(setIS, ierr)) end if ! setIS !! face-based contributions PetscCall(DMGetLabelIdIS(dmDisplacement, MEF90FaceSetLabelName, setIS, ierr)) PetscCall(MEF90ISAllGatherMerge(MEF90DefMechCtx%comm, setIS, ierr)) boundaryForceWork = 0.0_kr if (.not. PetscObjectIsNull(setIS)) then PetscCall(ISGetIndices(setIS, setID, ierr)) PetscCall(VecGetArray(MEF90DefMechCtx%boundaryForce, boundaryForceArray, ierr)) PetscCall(VecGetArray(MEF90DefMechCtx%pressureForce, pressureForceArray, ierr)) do set = 1, size(setID) myWork = 0.0_kr PetscCall(DMGetStratumIS(dmDisplacement, MEF90FaceSetLabelName, setID(set), setPointIS, ierr)) if (.not. PetscObjectIsNull(setPointIS)) then write (prefix, '(A,"fs",I4.4,"_")') trim(MEF90DefMechCtx%prefix), setID(set) PetscCall(MEF90DefMechFaceSetOptionsSetFromOptions(MEF90DefMechCtx%comm, prefix, faceSetOptions, ierr)) PetscCall(ISGetIndices(setPointIS, setPointID, ierr)) PetscCall(DMPlexGetCellType(dmDisplacement, setPointID(1), cellGeometry, ierr)) PetscCall(MEF90ElementGetTypeBoundary(MEF90CtxGlobalOptions%elementFamily, MEF90CtxGlobalOptions%elementOrder, cellGeometry, elemVectType, ierr)) !! Allocate elements QuadratureOrder = 2 * elemVectType%order PetscCall(MEF90ElementCreate(dmDisplacement, setPointIS, elemVect, QuadratureOrder, elemVectType, ierr)) numDofDisplacement = size(elemVect(1)%BF(:, 1)) numGauss = size(elemVect(1)%Gauss_C) if (norm2(faceSetOptions%boundaryForce) /= 0.0_kr) then do cell = 1, size(setPointID) PetscCall(PetscSectionGetOffset(sectionBoundaryForce, setPointID(cell), vecOffset, ierr)) boundaryForce = boundaryForceArray(vecOffset + 1:vecOffset + 1 + SIZEOFMEF90_VECT) PetscCall(DMPlexVecGetClosure(dmDisplacement, PETSC_NULL_SECTION, MEF90DefMechCtx%displacementLocal, setPointID(cell), PETSC_NULL_INTEGER, displacementDof, ierr)) do iGauss = 1, numGauss do iDof = 1, numDofDisplacement myWork = myWork + elemVect(cell)%Gauss_C(iGauss) * displacementDof(iDof) * (boundaryForce .DotP. elemVect(cell)%BF(iDof, iGauss)) end do ! iDof numDofDisplacement end do ! iGauss PetscCall(DMPlexVecRestoreClosure(dmDisplacement, PETSC_NULL_SECTION, MEF90DefMechCtx%displacementLocal, setPointID(cell), PETSC_NULL_INTEGER, displacementDof, ierr)) end do ! cell end if ! faceSetOptions%boundaryForce if (faceSetOptions%pressureForce /= 0.0_kr) then do cell = 1, size(setPointID) PetscCall(PetscSectionGetOffset(sectionPressureForce, setPointID(cell), vecOffset, ierr)) PetscCall(DMPlexVecGetClosure(dmDisplacement, PETSC_NULL_SECTION, MEF90DefMechCtx%displacementLocal, setPointID(cell), PETSC_NULL_INTEGER, displacementDof, ierr)) pressureForce = pressureForceArray(vecOffset + 1) * elemVect(cell)%outerNormal do iGauss = 1, numGauss do iDof = 1, numDofDisplacement myWork = myWork + elemVect(cell)%Gauss_C(iGauss) * displacementDof(iDof) * (pressureForce .DotP. elemVect(cell)%BF(iDof, iGauss)) end do ! iDof numDofDisplacement end do ! iGauss PetscCall(DMPlexVecRestoreClosure(dmDisplacement, PETSC_NULL_SECTION, MEF90DefMechCtx%displacementLocal, setPointID(cell), PETSC_NULL_INTEGER, displacementDof, ierr)) end do ! cell end if ! faceSetOptions%pressureForce PetscCall(MEF90ElementDestroy(elemVect, ierr)) PetscCall(ISRestoreIndices(setPointIS, setPointID, ierr)) PetscCall(ISDestroy(setPointIS, ierr)) end if ! setPointIS PetscCallMPI(MPI_AllReduce(myWork, boundaryForceWork(set), 1, MPIU_SCALAR, MPI_SUM, MEF90DefMechCtx%comm, ierr)) end do ! set PetscCall(VecRestoreArray(MEF90DefMechCtx%boundaryForce, boundaryForceArray, ierr)) PetscCall(VecRestoreArray(MEF90DefMechCtx%pressureForce, pressureForceArray, ierr)) PetscCall(ISDestroy(setIS, ierr)) end if ! setIS end subroutine MEF90DefMechWork #undef __FUNCT__ #define __FUNCT__ "MEF90DefMechCohesiveEnergy" !!! author: Blaise Bourdin (2014, bourdin@lsu.edu) !!! author: Blaise Bourdin (2022-26, bourdin@mcmaster.ca) !!! author: Alexis Marboeuf (2022, marboeua@mcmaster.ca) !!! !!! MEF90DefMechCohesiveEnergy: !!! subroutine MEF90DefMechCohesiveEnergy(MEF90DefMechCtx, cohesiveEnergy, ierr) type(MEF90DefMech_Type), intent(IN) :: MEF90DefMechCtx PetscReal, dimension(:), pointer :: cohesiveEnergy PetscErrorCode, intent(INOUT) :: ierr type(tDM) :: dmDisplacement, dmCohesiveDisplacement PetscReal, dimension(:), pointer :: displacementDof, cohesiveDisplacementDof type(tIS) :: setIS, setPointIS PetscInt, dimension(:), pointer :: setID, setPointID PetscInt :: set, QuadratureOrder, cell type(MEF90DefMechCellSetOptions_Type) :: cellSetOptions type(MEF90_ELEMENT_ELAST), dimension(:), pointer :: elemVect type(eDMPolytopeType) :: cellGeometry type(MEF90ElementType) :: elemVectType character(len=MEF90MXSTRLEN) :: prefix type(MEF90CtxGlobalOptions_Type) :: MEF90CtxGlobalOptions type(MEF90_VECT) :: U0Gauss PetscReal :: myCohesiveEnergy PetscInt :: iDof, iGauss, numDofDisplacement, numGauss PetscCall(MEF90CtxGlobalOptionsSetFromOptions(MEF90DefMechCtx%MEF90Ctx%comm, trim(MEF90DefMechCtx%MEF90Ctx%prefix), MEF90CtxGlobalOptions, ierr)) PetscCall(VecGetDM(MEF90DefMechCtx%displacementLocal, dmDisplacement, ierr)) PetscCall(VecGetDM(MEF90DefMechCtx%cohesiveDisplacement, dmCohesiveDisplacement, ierr)) !! get IS for cell sets PetscCall(DMGetLabelIdIS(dmDisplacement, MEF90CellSetLabelName, setIS, ierr)) PetscCall(MEF90ISAllGatherMerge(MEF90DefMechCtx%comm, setIS, ierr)) cohesiveEnergy = 0.0_kr if (.not. PetscObjectIsNull(setIS)) then PetscCall(ISGetIndices(setIS, setID, ierr)) do set = 1, size(setID) myCohesiveEnergy = 0.0_kr PetscCall(DMGetStratumIS(dmDisplacement, MEF90CellSetLabelName, setID(set), setPointIS, ierr)) if (.not. PetscObjectIsNull(setPointIS)) then write (prefix, '(A,"cs",I4.4,"_")') trim(MEF90DefMechCtx%prefix), setID(set) PetscCall(MEF90DefMechCellSetOptionsSetFromOptions(MEF90DefMechCtx%comm, prefix, MEF90_DIM, cellSetOptions, ierr)) PetscCall(ISGetIndices(setPointIS, setPointID, ierr)) PetscCall(DMPlexGetCellType(dmDisplacement, setPointID(1), cellGeometry, ierr)) PetscCall(MEF90ElementGetType(MEF90CtxGlobalOptions%elementFamily, MEF90CtxGlobalOptions%elementOrder, cellGeometry, elemVectType, ierr)) !! Allocate elements QuadratureOrder = 2 * (elemVectType%order - 1) PetscCall(MEF90ElementCreate(dmDisplacement, setPointIS, elemVect, QuadratureOrder, elemVectType, ierr)) numDofDisplacement = size(elemVect(1)%BF(:, 1)) numGauss = size(elemVect(1)%Gauss_C) if ((norm2(cellSetOptions%cohesiveDisplacement) /= 0.0_kr) .and. (cellSetOptions%cohesiveStiffness /= 0.0_kr)) then do cell = 1, size(setPointID) U0Gauss = 0.0_kr PetscCall(DMPlexVecGetClosure(dmDisplacement, PETSC_NULL_SECTION, MEF90DefMechCtx%displacementLocal, setPointID(cell), PETSC_NULL_INTEGER, displacementDof, ierr)) PetscCall(DMPlexVecGetClosure(dmCohesiveDisplacement, PETSC_NULL_SECTION, MEF90DefMechCtx%cohesiveDisplacement, setPointID(cell), PETSC_NULL_INTEGER, cohesiveDisplacementDof, ierr)) do iGauss = 1, numGauss do iDof = 1, numDofDisplacement U0Gauss = U0Gauss + (displacementDof(iDof) - cohesiveDisplacementDof(iDof)) * elemVect(cell)%BF(iDof, iGauss) end do ! iDof numDofDisplacement myCohesiveEnergy = myCohesiveEnergy + 0.5_kr * elemVect(cell)%Gauss_C(iGauss) * cellSetOptions%cohesiveStiffness * (U0Gauss .DotP. U0Gauss) end do ! iGauss PetscCall(DMPlexVecRestoreClosure(dmCohesiveDisplacement, PETSC_NULL_SECTION, MEF90DefMechCtx%cohesiveDisplacement, setPointID(cell), PETSC_NULL_INTEGER, cohesiveDisplacementDof, ierr)) PetscCall(DMPlexVecRestoreClosure(dmDisplacement, PETSC_NULL_SECTION, MEF90DefMechCtx%displacementLocal, setPointID(cell), PETSC_NULL_INTEGER, displacementDof, ierr)) end do ! cell end if ! cohesiveDisplacement PetscCall(MEF90ElementDestroy(elemVect, ierr)) PetscCall(ISRestoreIndices(setPointIS, setPointID, ierr)) PetscCall(ISDestroy(setPointIS, ierr)) end if ! setPointIS PetscCallMPI(MPI_AllReduce(myCohesiveEnergy, cohesiveEnergy(set), 1, MPIU_SCALAR, MPI_SUM, MEF90DefMechCtx%comm, ierr)) end do ! set PetscCall(ISRestoreIndices(setIS, setID, ierr)) PetscCall(ISDestroy(setIS, ierr)) end if ! setIS end subroutine MEF90DefMechCohesiveEnergy #undef __FUNCT__ #define __FUNCT__ "MEF90DefMechPlasticDissipation" !!! author: Erwan Tanne (2015, erwan.tanne@gmail.com) !!! !!! MEF90DefMechPlasticDissipation: !!! subroutine MEF90DefMechPlasticDissipation(MEF90DefMechCtx, plasticStrainOld, energy, ierr) type(tVec), intent(IN) :: plasticStrainOld type(MEF90DefMech_Type), intent(IN) :: MEF90DefMechCtx PetscReal, dimension(:), pointer :: energy PetscErrorCode, intent(INOUT) :: ierr end subroutine MEF90DefMechPlasticDissipation #undef __FUNCT__ #define __FUNCT__ "MEF90DefMechElasticEnergy" !!! author: Blaise Bourdin (2012-14, bourdin@lsu.edu) !!! author: Blaise Bourdin (2022-26, bourdin@mcmaster.ca) !!! author: Alexis Marboeuf (2022, marboeua@mcmaster.ca) !!! !!! MEF90DefMechElasticEnergy: !!! subroutine MEF90DefMechElasticEnergy(MEF90DefMechCtx, energy, ierr) type(MEF90DefMech_Type), intent(IN) :: MEF90DefMechCtx PetscReal, dimension(:), pointer :: energy PetscErrorCode, intent(INOUT) :: ierr type(tDM) :: dmDisplacement, dmDamage, dmTemperature, dmPlasticStrain type(tPetscSection) :: sectionPlasticStrain PetscReal, dimension(:), pointer :: displacementDof, damageDof, temperatureDof, plasticStrainArray PetscInt :: vecOffset type(MEF90_MATS) :: plasticStrainCell type(tIS) :: setIS, setPointIS PetscInt, dimension(:), pointer :: setID, setPointID PetscInt :: set, QuadratureOrder, cell, dim type(MEF90DefMechCellSetOptions_Type) :: cellSetOptions type(MEF90_MATS) :: linearThermalExpansion type(MEF90_ELEMENT_ELAST), dimension(:), pointer :: elemVect type(MEF90_ELEMENT_SCAL), dimension(:), pointer :: elemScal type(eDMPolytopeType) :: cellGeometryVect, cellGeometryScal type(MEF90ElementType) :: elemVectType, elemScalType type(MEF90CtxGlobalOptions_Type) :: MEF90CtxGlobalOptions type(MEF90DefMechGlobalOptions_Type) :: MEF90DefMechGlobalOptions class(MEF90DefMechAT_Type), allocatable :: ATModel class(MEF90DefMechSplit), allocatable :: Split type(MEF90_MATS) :: totalStrainGauss PetscReal :: damageGauss, temperatureGauss, myEnergy, EEDPlus, EEDMinus, elasticEnergyDensityGauss PetscInt :: iDof, iGauss, numDofDisplacement, numDofDamage, numGauss character(len=MEF90MXSTRLEN) :: prefix class(MEF90HookesLaw), allocatable :: HookesLaw type(tVec), pointer :: damageLocal PetscCall(MEF90CtxGlobalOptionsSetFromOptions(MEF90DefMechCtx%MEF90Ctx%comm, trim(MEF90DefMechCtx%MEF90Ctx%prefix), MEF90CtxGlobalOptions, ierr)) PetscCall(MEF90DefMechGlobalOptionsSetFromOptions(MEF90DefMechCtx%comm, trim(MEF90DefMechCtx%prefix), MEF90DefMechGlobalOptions, ierr)) PetscCall(VecGetDM(MEF90DefMechCtx%displacementLocal, dmDisplacement, ierr)) PetscCall(VecGetDM(MEF90DefMechCtx%damageLocal, dmDamage, ierr)) PetscCall(VecGetDM(MEF90DefMechCtx%TemperatureLocal, dmTemperature, ierr)) PetscCall(VecGetDM(MEF90DefMechCtx%plasticStrain, dmPlasticStrain, ierr)) PetscCall(DMGetLocalSection(dmPlasticStrain, sectionPlasticStrain, ierr)) PetscCall(DMGetDimension(dmDisplacement, dim, ierr)) !! get IS for cell sets PetscCall(DMGetLabelIdIS(dmDisplacement, MEF90CellSetLabelName, setIS, ierr)) PetscCall(MEF90ISAllGatherMerge(MEF90DefMechCtx%comm, setIS, ierr)) energy = 0.0_kr if (.not. PetscObjectIsNull(setIS)) then PetscCall(VecGetArray(MEF90DefMechCtx%plasticStrain, plasticStrainArray, ierr)) damageLocal => MEF90DefMechCtx%damageLocal PetscCall(ISGetIndices(setIS, setID, ierr)) do set = 1, size(setID) myEnergy = 0.0_kr if (MEF90DefMechGlobalOptions%multiPhaseField) then damageLocal => MEF90DefMechCtx%partialDamageLocal(set) end if PetscCall(DMGetStratumIS(dmDisplacement, MEF90CellSetLabelName, setID(set), setPointIS, ierr)) if (.not. PetscObjectIsNull(setPointIS)) then write (prefix, '(A,"cs",I4.4,"_")') trim(MEF90DefMechCtx%prefix), setID(set) PetscCall(MEF90DefMechCellSetOptionsSetFromOptions(MEF90DefMechCtx%comm, prefix, MEF90_DIM, cellSetOptions, ierr)) select type (linearThermalExpansionMat => cellSetOptions%linearThermalExpansion) type is (MEF90_MATS) linearThermalExpansion = linearThermalExpansionMat end select PetscCall(ISGetIndices(setPointIS, setPointID, ierr)) PetscCall(DMPlexGetCellType(dmDisplacement, setPointID(1), cellGeometryVect, ierr)) PetscCall(DMPlexGetCellType(dmDamage, setPointID(1), cellGeometryScal, ierr)) PetscCall(MEF90ElementGetType(MEF90CtxGlobalOptions%elementFamily, MEF90CtxGlobalOptions%elementOrder, cellGeometryVect, elemVectType, ierr)) PetscCall(MEF90ElementGetType(MEF90CtxGlobalOptions%elementFamily, MEF90CtxGlobalOptions%elementOrder, cellGeometryScal, elemScalType, ierr)) !! get the ATModel and split objects PetscCall(MEF90DefMechGetATModel(MEF90DefMechCtx%comm, prefix, dim, ATModel, ierr)) PetscCall(ATModel%setFromOptions(ierr)) PetscCall(MEF90DefMechGetSplit(MEF90DefMechCtx%comm, prefix, Split, ierr)) PetscCall(Split%setFromOptions(ierr)) PetscCall(MEF90GetHookesLaw(MEF90DefMechCtx%Comm, prefix, MEF90_DIM, HookesLaw, ierr)) PetscCall(HookesLaw%setFromOptions(ierr)) !! Allocate elements QuadratureOrder = max(2 * elemVectType%order, Split%quadratureOrder) PetscCall(MEF90ElementCreate(dmDisplacement, setPointIS, elemVect, QuadratureOrder, elemVectType, ierr)) PetscCall(MEF90ElementCreate(dmDamage, setPointIS, elemScal, QuadratureOrder, elemScalType, ierr)) numDofDisplacement = size(elemVect(1)%BF(:, 1)) numDofDamage = size(elemScal(1)%BF(:, 1)) numGauss = size(elemVect(1)%Gauss_C) do cell = 1, size(setPointID) PetscCall(DMPlexVecGetClosure(dmDamage, PETSC_NULL_SECTION, damageLocal, setPointID(cell), PETSC_NULL_INTEGER, damageDof, ierr)) PetscCall(DMPlexVecGetClosure(dmDisplacement, PETSC_NULL_SECTION, MEF90DefMechCtx%displacementLocal, setPointID(cell), PETSC_NULL_INTEGER, displacementDof, ierr)) PetscCall(DMPlexVecGetClosure(dmTemperature, PETSC_NULL_SECTION, MEF90DefMechCtx%TemperatureLocal, setPointID(cell), PETSC_NULL_INTEGER, temperatureDof, ierr)) PetscCall(PetscSectionGetOffset(sectionPlasticStrain, setPointID(cell), vecOffset, ierr)) plasticStrainCell = plasticStrainArray(vecOffset + 1:vecOffset + 1 + SIZEOFMEF90_MATS) do iGauss = 1, numGauss damageGauss = 0.0_kr if (.not. ATModel%isElastic) then do iDof = 1, numDofDamage damageGauss = damageGauss + damageDof(iDof) * elemScal(cell)%BF(iDof, iGauss) end do ! iDof numDofDamage end if totalStrainGauss = 0.0_kr do iDof = 1, numDofDisplacement totalStrainGauss = totalStrainGauss + displacementDof(iDof) * elemVect(cell)%GradS_BF(iDof, iGauss) end do ! iDof numDofDisplacement temperatureGauss = 0.0_kr do iDof = 1, numDofDamage temperatureGauss = temperatureGauss + temperatureDof(iDof) * elemScal(cell)%BF(iDof, iGauss) end do ! iDof numDofDamage totalStrainGauss = totalStrainGauss - (temperatureGauss * linearThermalExpansion) ! #if MEF90_DIM == 2 !! We need something along these lines !! Adding terms in planestrain for plasticity with tr(p) = 0 ! If (.NOT. matProp%HookesLaw%isPlaneStress) Then ! stressGauss = stressGauss + stiffness * ( matProp%HookesLaw%lambda*trace(plasticStrainCell)*MEF90MatS2DIdentity ) ! End If ! #endif !! Not actually sure if the split should be computed at totalStrainGauss or at totalStrainGauss - plasticStrainCell call Split%setup(totalStrainGauss - plasticStrainCell, ierr) call Split%EED(HookesLaw, totalStrainGauss - plasticStrainCell, EEDPlus, EEDMinus, ierr) if (ATModel%isElastic) then elasticEnergyDensityGauss = EEDPlus + EEDMinus else if (split%isHybrid) then elasticEnergyDensityGauss = ATModel%a(damageGauss) * (EEDPlus + EEDMinus) else elasticEnergyDensityGauss = ATModel%a(damageGauss) * EEDPlus + EEDMinus endif end if myEnergy = myEnergy + elemVect(cell)%Gauss_C(iGauss) * elasticEnergyDensityGauss end do ! iGauss PetscCall(DMPlexVecRestoreClosure(dmDamage, PETSC_NULL_SECTION, damageLocal, setPointID(cell), PETSC_NULL_INTEGER, damageDof, ierr)) PetscCall(DMPlexVecRestoreClosure(dmDisplacement, PETSC_NULL_SECTION, MEF90DefMechCtx%displacementLocal, setPointID(cell), PETSC_NULL_INTEGER, displacementDof, ierr)) PetscCall(DMPlexVecRestoreClosure(dmTemperature, PETSC_NULL_SECTION, MEF90DefMechCtx%TemperatureLocal, setPointID(cell), PETSC_NULL_INTEGER, temperatureDof, ierr)) end do ! cell PetscCall(MEF90ElementDestroy(elemVect, ierr)) PetscCall(MEF90ElementDestroy(elemScal, ierr)) PetscCall(ISRestoreIndices(setPointIS, setPointID, ierr)) PetscCall(ISDestroy(setPointIS, ierr)) end if ! setPointIS PetscCallMPI(MPI_AllReduce(myEnergy, energy(set), 1, MPIU_SCALAR, MPI_SUM, MEF90DefMechCtx%comm, ierr)) end do ! set PetscCall(ISRestoreIndices(setIS, setID, ierr)) PetscCall(VecRestoreArray(MEF90DefMechCtx%plasticStrain, plasticStrainArray, ierr)) PetscCall(ISDestroy(setIS, ierr)) end if ! setIS end subroutine MEF90DefMechElasticEnergy #undef __FUNCT__ #define __FUNCT__ "MEF90DefMechStress" !!! author: Blaise Bourdin (2012-20, bourdin@lsu.edu) !!! author: Blaise Bourdin (2022-26, bourdin@mcmaster.ca) !!! author: Alexis Marboeuf (2022, marboeua@mcmaster.ca) !!! !!! MEF90DefMechStress: !!! subroutine MEF90DefMechStress(MEF90DefMechCtx, stress, ierr) type(MEF90DefMech_Type), intent(IN) :: MEF90DefMechCtx PetscErrorCode, intent(INOUT) :: ierr type(tVec), intent(IN) :: stress type(tDM) :: dmDisplacement, dmDamage, dmTemperature, dmPlasticStrain, dmStress type(tPetscSection) :: sectionPlasticStrain PetscReal, dimension(:), pointer :: displacementDof, damageDof, temperatureDof, plasticStrainArray type(MEF90_MATS) :: plasticStrainCell type(tIS) :: setIS, setPointIS PetscInt, dimension(:), pointer :: setID, setPointID PetscInt :: set, QuadratureOrder, cell, dim, vecOffset type(MEF90DefMechCellSetOptions_Type) :: cellSetOptions type(MEF90_MATS) :: linearThermalExpansion type(MEF90_ELEMENT_ELAST), dimension(:), pointer :: elemVect type(MEF90_ELEMENT_SCAL), dimension(:), pointer :: elemScal type(eDMPolytopeType) :: cellGeometryVect, cellGeometryScal type(MEF90ElementType) :: elemVectType, elemScalType PetscReal, dimension(:), pointer :: stressDof type(MEF90CtxGlobalOptions_Type) :: MEF90CtxGlobalOptions type(MEF90DefMechGlobalOptions_Type) :: MEF90DefMechGlobalOptions class(MEF90DefMechAT_Type), allocatable :: ATModel class(MEF90DefMechSplit), allocatable :: Split PetscReal :: stressGaussPlus, stressGaussMinus type(MEF90_MATS) :: totalStrainGauss PetscReal :: damageGauss, temperatureGauss, cellSize PetscInt :: iDof, iGauss, numDofDisplacement, numDofDamage, numGauss, ij character(len=MEF90MXSTRLEN) :: prefix class(MEF90HookesLaw), allocatable :: HookesLaw type(tVec), pointer :: damageLocal #if MEF90_DIM == 2 type(MEF90_MATS), dimension(3) :: eij = [MatS2D(1.0_Kr, 0.0_Kr, 0.0_Kr), & MatS2D(0.0_Kr, 1.0_Kr, 0.0_Kr), & MatS2D(0.0_Kr, 0.0_Kr, 1.0_Kr)] #else type(MEF90_MATS), dimension(6) :: eij = [MatS3D(1.0_Kr, 0.0_Kr, 0.0_Kr, 0.0_Kr, 0.0_Kr, 0.0_Kr), & MatS3D(0.0_Kr, 1.0_Kr, 0.0_Kr, 0.0_Kr, 0.0_Kr, 0.0_Kr), & MatS3D(0.0_Kr, 0.0_Kr, 1.0_Kr, 0.0_Kr, 0.0_Kr, 0.0_Kr), & MatS3D(0.0_Kr, 0.0_Kr, 0.0_Kr, 1.0_Kr, 0.0_Kr, 0.0_Kr), & MatS3D(0.0_Kr, 0.0_Kr, 0.0_Kr, 0.0_Kr, 1.0_Kr, 0.0_Kr), & MatS3D(0.0_Kr, 0.0_Kr, 0.0_Kr, 0.0_Kr, 0.0_Kr, 1.0_Kr)] #endif PetscCall(MEF90CtxGlobalOptionsSetFromOptions(MEF90DefMechCtx%MEF90Ctx%comm, trim(MEF90DefMechCtx%MEF90Ctx%prefix), MEF90CtxGlobalOptions, ierr)) PetscCall(MEF90DefMechGlobalOptionsSetFromOptions(MEF90DefMechCtx%comm, trim(MEF90DefMechCtx%prefix), MEF90DefMechGlobalOptions, ierr)) PetscCall(VecGetDM(stress, dmStress, ierr)) PetscCall(VecGetDM(MEF90DefMechCtx%displacementLocal, dmDisplacement, ierr)) PetscCall(VecGetDM(MEF90DefMechCtx%damageLocal, dmDamage, ierr)) PetscCall(VecGetDM(MEF90DefMechCtx%TemperatureLocal, dmTemperature, ierr)) PetscCall(VecGetDM(MEF90DefMechCtx%plasticStrain, dmPlasticStrain, ierr)) PetscCall(DMGetLocalSection(dmPlasticStrain, sectionPlasticStrain, ierr)) PetscCall(DMGetDimension(dmDisplacement, dim, ierr)) PetscCall(VecSet(stress, 0.0_kr, ierr)) !! get IS for cell sets PetscCall(DMGetLabelIdIS(dmDisplacement, MEF90CellSetLabelName, setIS, ierr)) PetscCall(MEF90ISAllGatherMerge(MEF90DefMechCtx%comm, setIS, ierr)) if (.not. PetscObjectIsNull(setIS)) then damageLocal => MEF90DefMechCtx%damageLocal PetscCall(ISGetIndices(setIS, setID, ierr)) PetscCall(VecGetArray(MEF90DefMechCtx%plasticStrain, plasticStrainArray, ierr)) do set = 1, size(setID) if (MEF90DefMechGlobalOptions%multiPhaseField) then damageLocal => MEF90DefMechCtx%partialDamageLocal(set) end if PetscCall(DMGetStratumIS(dmDisplacement, MEF90CellSetLabelName, setID(set), setPointIS, ierr)) if (.not. PetscObjectIsNull(setPointIS)) then write (prefix, '(A,"cs",I4.4,"_")') trim(MEF90DefMechCtx%prefix), setID(set) PetscCall(MEF90DefMechCellSetOptionsSetFromOptions(MEF90DefMechCtx%comm, prefix, MEF90_DIM, cellSetOptions, ierr)) select type (linearThermalExpansionMat => cellSetOptions%linearThermalExpansion) type is (MEF90_MATS) linearThermalExpansion = linearThermalExpansionMat end select PetscCall(ISGetIndices(setPointIS, setPointID, ierr)) PetscCall(DMPlexGetCellType(dmDisplacement, setPointID(1), cellGeometryVect, ierr)) PetscCall(DMPlexGetCellType(dmDamage, setPointID(1), cellGeometryScal, ierr)) PetscCall(MEF90ElementGetType(MEF90CtxGlobalOptions%elementFamily, MEF90CtxGlobalOptions%elementOrder, cellGeometryVect, elemVectType, ierr)) PetscCall(MEF90ElementGetType(MEF90CtxGlobalOptions%elementFamily, MEF90CtxGlobalOptions%elementOrder, cellGeometryScal, elemScalType, ierr)) !! get the ATModel and split objects PetscCall(MEF90DefMechGetATModel(MEF90DefMechCtx%comm, prefix, dim, ATModel, ierr)) PetscCall(ATModel%setFromOptions(ierr)) PetscCall(MEF90DefMechGetSplit(MEF90DefMechCtx%comm, prefix, Split, ierr)) PetscCall(Split%setFromOptions(ierr)) PetscCall(MEF90GetHookesLaw(MEF90DefMechCtx%Comm, prefix, MEF90_DIM, HookesLaw, ierr)) PetscCall(HookesLaw%setFromOptions(ierr)) !! Allocate elements QuadratureOrder = max(2 * elemVectType%order, Split%quadratureOrder) PetscCall(MEF90ElementCreate(dmDisplacement, setPointIS, elemVect, QuadratureOrder, elemVectType, ierr)) PetscCall(MEF90ElementCreate(dmDamage, setPointIS, elemScal, QuadratureOrder, elemScalType, ierr)) numDofDisplacement = size(elemVect(1)%BF(:, 1)) numDofDamage = size(elemScal(1)%BF(:, 1)) numGauss = size(elemVect(1)%Gauss_C) allocate (stressDof(SIZEOFMEF90_MATS)) do cell = 1, size(setPointID) PetscCall(DMPlexVecGetClosure(dmDamage, PETSC_NULL_SECTION, damageLocal, setPointID(cell), PETSC_NULL_INTEGER, damageDof, ierr)) PetscCall(DMPlexVecGetClosure(dmDisplacement, PETSC_NULL_SECTION, MEF90DefMechCtx%displacementLocal, setPointID(cell), PETSC_NULL_INTEGER, displacementDof, ierr)) PetscCall(DMPlexVecGetClosure(dmTemperature, PETSC_NULL_SECTION, MEF90DefMechCtx%TemperatureLocal, setPointID(cell), PETSC_NULL_INTEGER, temperatureDof, ierr)) stressDof = 0.0_kr cellSize = 0.0_kr do iGauss = 1, numGauss damageGauss = 0.0_kr if (.not. ATModel%isElastic) then do iDof = 1, numDofDamage damageGauss = damageGauss + damageDof(iDof) * elemScal(cell)%BF(iDof, iGauss) end do ! iDof numDofDamage end if totalStrainGauss = 0.0_kr do iDof = 1, numDofDisplacement totalStrainGauss = totalStrainGauss + displacementDof(iDof) * elemVect(cell)%GradS_BF(iDof, iGauss) end do ! iDof numDofDisplacement temperatureGauss = 0.0_kr do iDof = 1, numDofDamage temperatureGauss = temperatureGauss + temperatureDof(iDof) * elemScal(cell)%BF(iDof, iGauss) end do ! iDof numDofDamage totalStrainGauss = totalStrainGauss - (temperatureGauss * linearThermalExpansion) ! #if MEF90_DIM == 2 !! We need something along these lines !! Adding terms in planestrain for plasticity with tr(p) = 0 ! If (.NOT. matProp%HookesLaw%isPlaneStress) Then ! stressGauss = stressGauss + stiffness * ( matProp%HookesLaw%lambda*trace(plasticStrainCell)*MEF90MatS2DIdentity ) ! End If ! #endif PetscCall(PetscSectionGetOffset(sectionPlasticStrain, setPointID(cell), vecOffset, ierr)) plasticStrainCell = plasticStrainArray(vecOffset + 1:vecOffset + 1 + SIZEOFMEF90_MATS) !! This is really silly but since I decided that EED, DEED, and D2EED would return directinal derivatives instead of linear operators, !! I have to call DEED for each component of the stress tensor (plus it only works for a linear law) do ij = 1, SIZEOFMEF90_MATS call Split%setup(totalStrainGauss - plasticStrainCell, ierr) call Split%DEED(HookesLaw, eij(ij), stressGaussPlus, stressGaussMinus, ierr) if (ATModel%isElastic) then stressDof(ij) = stressDof(ij) + elemVect(cell)%Gauss_C(iGauss) * (stressGaussPlus + stressGaussMinus) else if (split%isHybrid) then stressDof(ij) = stressDof(ij) + elemVect(cell)%Gauss_C(iGauss) * ATModel%a(damageGauss) * (stressGaussPlus + stressGaussMinus) else stressDof(ij) = stressDof(ij) + elemVect(cell)%Gauss_C(iGauss) * (ATModel%a(damageGauss) * stressGaussPlus + stressGaussMinus) end if ! isHybrid end if ! isElastic end do cellSize = cellSize + elemVect(cell)%Gauss_C(iGauss) end do ! iGauss stressDof = stressDof / cellSize PetscCall(DMPlexVecSetClosure(dmStress, PETSC_NULL_SECTION, stress, setPointID(cell), stressDof, INSERT_VALUES, ierr)) PetscCall(DMPlexVecRestoreClosure(dmDamage, PETSC_NULL_SECTION, MEF90DefMechCtx%damageLocal, setPointID(cell), PETSC_NULL_INTEGER, damageDof, ierr)) PetscCall(DMPlexVecRestoreClosure(dmDisplacement, PETSC_NULL_SECTION, MEF90DefMechCtx%displacementLocal, setPointID(cell), PETSC_NULL_INTEGER, displacementDof, ierr)) PetscCall(DMPlexVecRestoreClosure(dmTemperature, PETSC_NULL_SECTION, MEF90DefMechCtx%TemperatureLocal, setPointID(cell), PETSC_NULL_INTEGER, temperatureDof, ierr)) end do ! cell deallocate (stressDof) PetscCall(MEF90ElementDestroy(elemVect, ierr)) PetscCall(MEF90ElementDestroy(elemScal, ierr)) PetscCall(ISRestoreIndices(setPointIS, setPointID, ierr)) PetscCall(ISDestroy(setPointIS, ierr)) end if ! setPointIS end do ! set PetscCall(ISRestoreIndices(setIS, setID, ierr)) PetscCall(VecRestoreArray(MEF90DefMechCtx%plasticStrain, plasticStrainArray, ierr)) PetscCall(ISDestroy(setIS, ierr)) end if ! setIS end subroutine MEF90DefMechStress #undef __FUNCT__ #define __FUNCT__ "MEF90DefMechOperatorDamage" !!! author: Blaise Bourdin (2012-20, bourdin@lsu.edu) !!! author: Erwan Tanne (2012-20, erwan.tanne@gmail.com) !!! author: Blaise Bourdin (2022-26, bourdin@mcmaster.ca) !!! author: Alexis Marboeuf (2022, marboeua@mcmaster.ca) !!! !!! MEF90DefMechOperatorDamage: Build the operator. When called in SNES, the solution time should always match the target time, !!! so there is no need for interpolation of the forces, external, and boundary values !!! subroutine MEF90DefMechOperatorDamage(snesDamage, damage, residual, MEF90DefMechCtx, ierr) type(tSNES), intent(IN) :: snesDamage type(tVec), intent(IN) :: damage type(tVec), intent(INOUT) :: residual type(MEF90DefMech_Type), intent(IN) :: MEF90DefMechCtx PetscErrorCode, intent(INOUT) :: ierr type(tDM) :: dmDisplacement, dmDamage, dmTemperature, dmCohesiveDisplacement, dmPlasticStrain type(tPetscSection) :: sectionPlasticStrain PetscReal, dimension(:), pointer :: displacementDof, damageDof, temperatureDof, plasticStrainArray type(MEF90_MATS) :: plasticStrainCell type(tIS) :: setIS, setPointIS PetscInt, dimension(:), pointer :: setID, setPointID PetscInt :: set, QuadratureOrder, cell, dim, vecOffset type(MEF90DefMechCellSetOptions_Type) :: cellSetOptions type(MEF90_MATS) :: linearThermalExpansion type(MEF90_ELEMENT_ELAST), dimension(:), pointer :: elemVect type(MEF90_ELEMENT_SCAL), dimension(:), pointer :: elemScal type(eDMPolytopeType) :: cellGeometryVect, cellGeometryScal type(MEF90ElementType) :: elemVectType, elemScalType PetscReal, dimension(:), pointer :: residualDof type(MEF90CtxGlobalOptions_Type) :: MEF90CtxGlobalOptions type(MEF90DefMechGlobalOptions_Type) :: MEF90DefMechGlobalOptions type(tVec) :: residualLoc class(MEF90DefMechAT_Type), allocatable :: ATModel class(MEF90DefMechSplit), allocatable :: Split type(MEF90_MATS) :: totalStrainGauss, C2 type(MEF90_VECT) :: gradDamageGauss PetscReal :: damageGauss, temperatureGauss, EEDGaussMinus, EEDGaussPlus, C1, C3 PetscInt :: iDof, iGauss, numDofDisplacement, numDofDamage, numGauss character(len=MEF90MXSTRLEN) :: prefix class(MEF90HookesLaw), allocatable :: HookesLaw type(tVec), pointer :: damageLocal PetscCall(MEF90CtxGlobalOptionsSetFromOptions(MEF90DefMechCtx%MEF90Ctx%comm, trim(MEF90DefMechCtx%MEF90Ctx%prefix), MEF90CtxGlobalOptions, ierr)) PetscCall(MEF90DefMechGlobalOptionsSetFromOptions(MEF90DefMechCtx%comm, trim(MEF90DefMechCtx%prefix), MEF90DefMechGlobalOptions, ierr)) PetscCall(VecGetDM(MEF90DefMechCtx%displacementLocal, dmDisplacement, ierr)) PetscCall(VecGetDM(MEF90DefMechCtx%cohesiveDisplacement, dmCohesiveDisplacement, ierr)) PetscCall(VecGetDM(MEF90DefMechCtx%damageLocal, dmDamage, ierr)) PetscCall(VecGetDM(MEF90DefMechCtx%TemperatureLocal, dmTemperature, ierr)) PetscCall(VecGetDM(MEF90DefMechCtx%plasticStrain, dmPlasticStrain, ierr)) PetscCall(DMGetLocalSection(dmPlasticStrain, sectionPlasticStrain, ierr)) PetscCall(DMGetDimension(dmDisplacement, dim, ierr)) PetscCall(DMGetLocalVector(dmDamage, residualLoc, ierr)) !! Something subtle is going on here: !! I _have_ to use MEF90DefMechCtx%damageLocal because I need the constrained values, which would not be initialized if I were !! to create a new local Vec, or duplicate MEF90DefMechCtx%damageLocal !! This DOES create a side effect. if (MEF90DefMechGlobalOptions%multiPhaseField) then damageLocal => MEF90DefMechCtx%partialDamageLocal(MEF90DefMechCtx%currentSet) else damageLocal => MEF90DefMechCtx%damageLocal end if PetscCall(DMGlobalToLocal(dmDamage, damage, INSERT_VALUES, damageLocal, ierr)) PetscCall(VecSet(residual, 0.0_kr, ierr)) PetscCall(VecSet(residualLoc, 0.0_kr, ierr)) !! get IS for cell sets PetscCall(DMGetLabelIdIS(dmDamage, MEF90CellSetLabelName, setIS, ierr)) PetscCall(MEF90ISAllGatherMerge(MEF90DefMechCtx%comm, setIS, ierr)) if (.not. PetscObjectIsNull(setIS)) then PetscCall(ISGetIndices(setIS, setID, ierr)) PetscCall(VecGetArray(MEF90DefMechCtx%plasticStrain, plasticStrainArray, ierr)) do set = 1, size(setID) PetscCall(DMGetStratumIS(dmDamage, MEF90CellSetLabelName, setID(set), setPointIS, ierr)) if (.not. PetscObjectIsNull(setPointIS)) then write (prefix, '(A,"cs",I4.4,"_")') trim(MEF90DefMechCtx%prefix), setID(set) PetscCall(MEF90DefMechCellSetOptionsSetFromOptions(MEF90DefMechCtx%comm, prefix, MEF90_DIM, cellSetOptions, ierr)) select type (linearThermalExpansionMat => cellSetOptions%linearThermalExpansion) type is (MEF90_MATS) linearThermalExpansion = linearThermalExpansionMat end select PetscCall(ISGetIndices(setPointIS, setPointID, ierr)) PetscCall(DMPlexGetCellType(dmDisplacement, setPointID(1), cellGeometryVect, ierr)) PetscCall(DMPlexGetCellType(dmDamage, setPointID(1), cellGeometryScal, ierr)) PetscCall(MEF90ElementGetType(MEF90CtxGlobalOptions%elementFamily, MEF90CtxGlobalOptions%elementOrder, cellGeometryVect, elemVectType, ierr)) PetscCall(MEF90ElementGetType(MEF90CtxGlobalOptions%elementFamily, MEF90CtxGlobalOptions%elementOrder, cellGeometryScal, elemScalType, ierr)) !! get the ATModel and split objects PetscCall(MEF90DefMechGetATModel(MEF90DefMechCtx%comm, prefix, dim, ATModel, ierr)) PetscCall(ATModel%setFromOptions(ierr)) PetscCall(MEF90DefMechGetSplit(MEF90DefMechCtx%comm, prefix, Split, ierr)) PetscCall(Split%setFromOptions(ierr)) PetscCall(MEF90GetHookesLaw(MEF90DefMechCtx%Comm, prefix, MEF90_DIM, HookesLaw, ierr)) PetscCall(HookesLaw%setFromOptions(ierr)) !! Allocate elements QuadratureOrder = max(2 * elemVectType%order, Split%quadratureOrder) PetscCall(MEF90ElementCreate(dmDisplacement, setPointIS, elemVect, QuadratureOrder, elemVectType, ierr)) PetscCall(MEF90ElementCreate(dmDamage, setPointIS, elemScal, QuadratureOrder, elemScalType, ierr)) numDofDisplacement = size(elemVect(1)%BF(:, 1)) numDofDamage = size(elemScal(1)%BF(:, 1)) numGauss = size(elemScal(1)%Gauss_C) allocate (residualDof(numDofDamage)) do cell = 1, size(setPointID) residualDof = 0.0_kr PetscCall(DMPlexVecGetClosure(dmDamage, PETSC_NULL_SECTION, damageLocal, setPointID(cell), PETSC_NULL_INTEGER, damageDof, ierr)) PetscCall(DMPlexVecGetClosure(dmDisplacement, PETSC_NULL_SECTION, MEF90DefMechCtx%displacementLocal, setPointID(cell), PETSC_NULL_INTEGER, displacementDof, ierr)) PetscCall(DMPlexVecGetClosure(dmTemperature, PETSC_NULL_SECTION, MEF90DefMechCtx%TemperatureLocal, setPointID(cell), PETSC_NULL_INTEGER, temperatureDof, ierr)) do iGauss = 1, numGauss damageGauss = 0.0_kr gradDamageGauss = 0.0_kr do iDof = 1, numDofDamage damageGauss = damageGauss + damageDof(iDof) * elemScal(cell)%BF(iDof, iGauss) gradDamageGauss = gradDamageGauss + damageDof(iDof) * elemScal(cell)%Grad_BF(iDof, iGauss) end do ! iDof numDofDamage totalStrainGauss = 0.0_kr do iDof = 1, numDofDisplacement totalStrainGauss = totalStrainGauss + displacementDof(iDof) * elemVect(cell)%GradS_BF(iDof, iGauss) end do ! iDof numDofDisplacement temperatureGauss = 0.0_kr do iDof = 1, numDofDamage temperatureGauss = temperatureGauss + temperatureDof(iDof) * elemScal(cell)%BF(iDof, iGauss) end do ! iDof numDofDamage totalStrainGauss = totalStrainGauss - (temperatureGauss * linearThermalExpansion) ! #if MEF90_DIM == 2 !! We need something along these lines !! Adding terms in planestrain for plasticity with tr(p) = 0 ! If (.NOT. matProp%HookesLaw%isPlaneStress) Then ! stressGauss = stressGauss + stiffness * ( matProp%HookesLaw%lambda*trace(plasticStrainCell)*MEF90MatS2DIdentity ) ! End If ! #endif PetscCall(PetscSectionGetOffset(sectionPlasticStrain, setPointID(cell), vecOffset, ierr)) plasticStrainCell = plasticStrainArray(vecOffset + 1:vecOffset + 1 + SIZEOFMEF90_MATS) !! begin ugly hack !! The following won't work until I overload the LinAlg operations to work with the parent classes ! C2 = ATModel%fractureToughness / ATModel%cw * 0.5_kr * ATModel%internalLength * ATModel%toughnessAnisotropyMatrix select type (k => ATModel%toughnessAnisotropyMatrix) type is (MEF90_MATS) C2 = ATModel%fractureToughness / ATModel%cw * 0.5_kr * ATModel%internalLength * k end select !! end ugly hack EEDGaussPlus = 0.0_kr if ((.not. ATModel%isElastic) .and. ((.not. (MEF90DefMechGlobalOptions%multiPhaseField)) .or. (MEF90DefMechCtx%currentSet == set))) then call Split%setup(totalStrainGauss - plasticStrainCell, ierr) call Split%EED(HookesLaw, totalStrainGauss - plasticStrainCell, EEDGaussPlus, EEDGaussMinus, ierr) end if C1 = ATModel%fractureToughness / ATModel%cw * 0.25_kr / ATModel%internalLength C3 = ATModel%Da(damageGauss) * EEDGaussPlus + C1 * ATModel%Dw(damageGauss) do iDof = 1, numDofDamage residualDof(iDof) = residualDof(iDof) + elemScal(cell)%Gauss_C(iGauss) * ( & C3 * elemScal(cell)%BF(iDof, iGauss) & + (C2 * gradDamageGauss .DotP. elemScal(cell)%Grad_BF(iDof, iGauss)) & ) end do ! iDof numDofDamage end do ! iGauss PetscCall(DMPlexVecRestoreClosure(dmTemperature, PETSC_NULL_SECTION, MEF90DefMechCtx%TemperatureLocal, setPointID(cell), PETSC_NULL_INTEGER, temperatureDof, ierr)) PetscCall(DMPlexVecRestoreClosure(dmDisplacement, PETSC_NULL_SECTION, MEF90DefMechCtx%displacementLocal, setPointID(cell), PETSC_NULL_INTEGER, displacementDof, ierr)) PetscCall(DMPlexVecRestoreClosure(dmDamage, PETSC_NULL_SECTION, damageLocal, setPointID(cell), PETSC_NULL_INTEGER, damageDof, ierr)) PetscCall(DMPlexVecSetClosure(dmDamage, PETSC_NULL_SECTION, residualLoc, setPointID(cell), residualDof, ADD_VALUES, ierr)) end do ! cellresidualLoc ! !!! crack pressure ! If (cellSetOptions%crackPressure /= 0.0_Kr)) Then ! Do cell = 1,size(setPointID) ! residualDof = 0.0_Kr ! Do iGauss = 1,numGauss ! CrackPressureCell = 0.0_Kr ! PetscCall(DMPlexVecGetClosure(dmCrackPressure,PETSC_NULL_SECTION,MEF90DefMechCtx%crackPressure,setPointID(cell),PETSC_NULL_INTEGER,crackPressureDof,ierr)) ! Do iDof = 1,numDofDamage ! CrackPressureCell = CrackPressureCell + crackPressureDof(iDof) * elemScal(cell)%BF(iDof,iGauss) ! End Do ! iDof numDofDamage ! PetscCall(DMPlexVecRestoreClosure(dmCrackPressure,PETSC_NULL_SECTION,MEF90DefMechCtx%crackPressure,setPointID(cell),PETSC_NULL_INTEGER,crackPressureDof,ierr)) ! displacementCell = 0.0_Kr ! PetscCall(DMPlexVecGetClosure(dmDisplacement,PETSC_NULL_SECTION,MEF90DefMechCtx%displacementLocal,setPointID(cell),PETSC_NULL_INTEGER,displacementDof,ierr)) ! Do iDof = 1,numDofDisplacement ! displacementCell = displacementCell + displacementDof(iDof) * elemVect(cell)%BF(iDof,iGauss) ! End Do ! iDof numDofDisplacement ! PetscCall(DMPlexVecRestoreClosure(dmDisplacement,PETSC_NULL_SECTION,MEF90DefMechCtx%displacementLocal,setPointID(cell),PETSC_NULL_INTEGER,displacementDof,ierr)) ! Do iDof = 1,numDofDamage ! residualDof(iDof) = residualDof(iDof) + elemScal(cell)%Gauss_C(iGauss) * & ! CrackPressureCell * (displacementCell .DotP. elemScal(cell)%Grad_BF(iDof,iGauss)) ! End Do ! iDof numDofDamage ! End Do ! iGauss ! PetscCall(DMPlexVecSetClosure(dmDamage,PETSC_NULL_SECTION,residualLoc,setPointID(cell),residualDof,ADD_VALUES,ierr)) ! End Do ! cell ! End If ! crack Pressure deallocate (residualDof) PetscCall(MEF90ElementDestroy(elemVect, ierr)) PetscCall(MEF90ElementDestroy(elemScal, ierr)) PetscCall(ISRestoreIndices(setPointIS, setPointID, ierr)) PetscCall(ISDestroy(setPointIS, ierr)) end if ! setPointIS end do ! set PetscCall(ISRestoreIndices(setIS, setID, ierr)) PetscCall(VecRestoreArray(MEF90DefMechCtx%plasticStrain, plasticStrainArray, ierr)) PetscCall(ISDestroy(setIS, ierr)) end if ! setIS PetscCall(DMLocalToGlobalBegin(dmDamage, residualLoc, ADD_VALUES, residual, ierr)) PetscCall(DMLocalToGlobalEnd(dmDamage, residualLoc, ADD_VALUES, residual, ierr)) PetscCall(DMRestoreLocalVector(dmDamage, residualLoc, ierr)) end subroutine MEF90DefMechOperatorDamage #undef __FUNCT__ #define __FUNCT__ "MEF90DefMechTAOGradientDamage" !!! author: Alexis Marboeuf (2012-2022, marboeua@mcmaster.ca) !!! author: Blaise Bourdin (2023, bourdin@mcmaster.ca) !!! !!! MEF90DefMechTAOGradientDamage: !!! subroutine MEF90DefMechTAOGradientDamage(taoDamage, damage, residual, MEF90DefMechCtx, ierr) type(tTao), intent(IN) :: taoDamage type(tVec), intent(IN) :: damage type(tVec), intent(INOUT) :: residual type(MEF90DefMech_Type), intent(IN) :: MEF90DefMechCtx PetscErrorCode, intent(INOUT) :: ierr type(tSNES) :: dummySNES ! type(tDM) :: dmDamage ! PetscCall(VecGetDM(MEF90DefMechCtx%damageLocal(1), dmDamage, ierr)) ! PetscCall(DMGlobalToLocal(dmDamage, damage, INSERT_VALUES, MEF90DefMechCtx%damageLocal(1), ierr)) PetscCall(MEF90DefMechOperatorDamage(dummySNES, damage, residual, MEF90DefMechCtx, ierr)) end subroutine MEF90DefMechTAOGradientDamage #undef __FUNCT__ #define __FUNCT__ "MEF90DefMechBilinearFormDamage" !!! author: Blaise Bourdin (2012-19, bourdin@lsu.edu) !!! author: Erwan Tanne (2012-19, erwan.tanne@gmail.com) !!! author: Blaise Bourdin (2022, bourdin@mcmaster.ca) !!! author: Alexis Marboeuf (2022, marboeua@mcmaster.ca) !!! !!! MEF90DefMechBilinearFormDamage: !!! subroutine MEF90DefMechBilinearFormDamage(snesDamage, damage, A, M, MEF90DefMechCtx, ierr) type(tSNES), intent(IN) :: snesDamage type(tVec), intent(IN) :: damage type(tMat), intent(INOUT) :: A, M type(MEF90DefMech_Type), intent(IN) :: MEF90DefMechCtx PetscErrorCode, intent(INOUT) :: ierr type(tDM) :: dmDisplacement, dmDamage, dmTemperature, dmCohesiveDisplacement, dmPlasticStrain type(tPetscSection) :: sectionPlasticStrain PetscReal, dimension(:), pointer :: displacementDof, damageDof, temperatureDof, plasticStrainArray type(MEF90_MATS) :: plasticStrainCell type(tIS) :: setIS, setPointIS PetscInt, dimension(:), pointer :: setID, setPointID PetscInt :: set, QuadratureOrder, cell, dim, vecOffset type(MEF90DefMechCellSetOptions_Type) :: cellSetOptions type(MEF90_MATS) :: linearThermalExpansion type(MEF90_ELEMENT_ELAST), dimension(:), pointer :: elemVect type(MEF90_ELEMENT_SCAL), dimension(:), pointer :: elemScal type(eDMPolytopeType) :: cellGeometryVect, cellGeometryScal type(MEF90ElementType) :: elemVectType, elemScalType PetscReal, dimension(:,:), pointer :: matDof type(MEF90CtxGlobalOptions_Type) :: MEF90CtxGlobalOptions type(MEF90DefMechGlobalOptions_Type) :: MEF90DefMechGlobalOptions class(MEF90DefMechAT_Type), allocatable :: ATModel class(MEF90DefMechSplit), allocatable :: Split type(MEF90_MATS) :: totalStrainGauss, C2 PetscReal :: damageGauss, temperatureGauss, EEDGaussMinus, EEDGaussPlus, C1, C3 PetscInt :: iDof, jDof, iGauss, numDofDisplacement, numDofDamage, numGauss character(len=MEF90MXSTRLEN) :: prefix class(MEF90HookesLaw), allocatable :: HookesLaw type(tVec), pointer :: damageLocal PetscCall(MEF90CtxGlobalOptionsSetFromOptions(MEF90DefMechCtx%MEF90Ctx%comm, trim(MEF90DefMechCtx%MEF90Ctx%prefix), MEF90CtxGlobalOptions, ierr)) PetscCall(MEF90DefMechGlobalOptionsSetFromOptions(MEF90DefMechCtx%comm, trim(MEF90DefMechCtx%prefix), MEF90DefMechGlobalOptions, ierr)) PetscCall(VecGetDM(MEF90DefMechCtx%displacementLocal, dmDisplacement, ierr)) PetscCall(VecGetDM(MEF90DefMechCtx%cohesiveDisplacement, dmCohesiveDisplacement, ierr)) PetscCall(VecGetDM(MEF90DefMechCtx%damageLocal, dmDamage, ierr)) PetscCall(VecGetDM(MEF90DefMechCtx%TemperatureLocal, dmTemperature, ierr)) PetscCall(VecGetDM(MEF90DefMechCtx%plasticStrain, dmPlasticStrain, ierr)) PetscCall(DMGetLocalSection(dmPlasticStrain, sectionPlasticStrain, ierr)) PetscCall(DMGetDimension(dmDisplacement, dim, ierr)) !! Something subtle is going on here: !! I _have_ to use MEF90DefMechCtx%damageLocal because I need the constrained values, which would not be initialized if I were !! to create a new local Vec, or duplicate MEF90DefMechCtx%damageLocal !! This DOES create a side effect. if (MEF90DefMechGlobalOptions%multiPhaseField) then damageLocal => MEF90DefMechCtx%partialDamageLocal(MEF90DefMechCtx%currentSet) else damageLocal => MEF90DefMechCtx%damageLocal end if PetscCall(DMGlobalToLocal(dmDamage, damage, INSERT_VALUES, damageLocal, ierr)) PetscCall(MatZeroEntries(A, ierr)) !! get IS for cell sets PetscCall(DMGetLabelIdIS(dmDamage, MEF90CellSetLabelName, setIS, ierr)) PetscCall(MEF90ISAllGatherMerge(MEF90DefMechCtx%comm, setIS, ierr)) if (.not. PetscObjectIsNull(setIS)) then PetscCall(ISGetIndices(setIS, setID, ierr)) PetscCall(VecGetArray(MEF90DefMechCtx%plasticStrain, plasticStrainArray, ierr)) do set = 1, size(setID) PetscCall(DMGetStratumIS(dmDamage, MEF90CellSetLabelName, setID(set), setPointIS, ierr)) if (.not. PetscObjectIsNull(setPointIS)) then write (prefix, '(A,"cs",I4.4,"_")') trim(MEF90DefMechCtx%prefix), setID(set) PetscCall(MEF90DefMechCellSetOptionsSetFromOptions(MEF90DefMechCtx%comm, prefix, MEF90_DIM, cellSetOptions, ierr)) select type (linearThermalExpansionMat => cellSetOptions%linearThermalExpansion) type is (MEF90_MATS) linearThermalExpansion = linearThermalExpansionMat end select PetscCall(ISGetIndices(setPointIS, setPointID, ierr)) PetscCall(DMPlexGetCellType(dmDisplacement, setPointID(1), cellGeometryVect, ierr)) PetscCall(DMPlexGetCellType(dmDamage, setPointID(1), cellGeometryScal, ierr)) PetscCall(MEF90ElementGetType(MEF90CtxGlobalOptions%elementFamily, MEF90CtxGlobalOptions%elementOrder, cellGeometryVect, elemVectType, ierr)) PetscCall(MEF90ElementGetType(MEF90CtxGlobalOptions%elementFamily, MEF90CtxGlobalOptions%elementOrder, cellGeometryScal, elemScalType, ierr)) !! get the ATModel and split objects PetscCall(MEF90DefMechGetATModel(MEF90DefMechCtx%comm, prefix, dim, ATModel, ierr)) PetscCall(ATModel%setFromOptions(ierr)) PetscCall(MEF90DefMechGetSplit(MEF90DefMechCtx%comm, prefix, Split, ierr)) PetscCall(Split%setFromOptions(ierr)) PetscCall(MEF90GetHookesLaw(MEF90DefMechCtx%Comm, prefix, MEF90_DIM, HookesLaw, ierr)) PetscCall(HookesLaw%setFromOptions(ierr)) !! Allocate elements QuadratureOrder = max(2 * elemVectType%order, Split%quadratureOrder) PetscCall(MEF90ElementCreate(dmDisplacement, setPointIS, elemVect, QuadratureOrder, elemVectType, ierr)) PetscCall(MEF90ElementCreate(dmDamage, setPointIS, elemScal, QuadratureOrder, elemScalType, ierr)) numDofDisplacement = size(elemVect(1)%BF(:, 1)) numDofDamage = size(elemScal(1)%BF(:, 1)) numGauss = size(elemScal(1)%Gauss_C) allocate (matDof(numDofDamage, numDofDamage)) do cell = 1, size(setPointID) matDof = 0.0_kr PetscCall(DMPlexVecGetClosure(dmDamage, PETSC_NULL_SECTION, damageLocal, setPointID(cell), PETSC_NULL_INTEGER, damageDof, ierr)) PetscCall(DMPlexVecGetClosure(dmDisplacement, PETSC_NULL_SECTION, MEF90DefMechCtx%displacementLocal, setPointID(cell), PETSC_NULL_INTEGER, displacementDof, ierr)) PetscCall(DMPlexVecGetClosure(dmTemperature, PETSC_NULL_SECTION, MEF90DefMechCtx%TemperatureLocal, setPointID(cell), PETSC_NULL_INTEGER, temperatureDof, ierr)) do iGauss = 1, numGauss damageGauss = 0.0_kr do iDof = 1, numDofDamage damageGauss = damageGauss + damageDof(iDof) * elemScal(cell)%BF(iDof, iGauss) end do ! iDof numDofDamage totalStrainGauss = 0.0_kr do iDof = 1, numDofDisplacement totalStrainGauss = totalStrainGauss + displacementDof(iDof) * elemVect(cell)%GradS_BF(iDof, iGauss) end do ! iDof numDofDisplacement temperatureGauss = 0.0_kr do iDof = 1, numDofDamage temperatureGauss = temperatureGauss + temperatureDof(iDof) * elemScal(cell)%BF(iDof, iGauss) end do ! iDof numDofDamage totalStrainGauss = totalStrainGauss - (temperatureGauss * linearThermalExpansion) ! #if MEF90_DIM == 2 !! We need something along these lines !! Adding terms in planestrain for plasticity with tr(p) = 0 ! If (.NOT. matProp%HookesLaw%isPlaneStress) Then ! stressGauss = stressGauss + stiffness * ( matProp%HookesLaw%lambda*trace(plasticStrainCell)*MEF90MatS2DIdentity ) ! End If ! #endif PetscCall(PetscSectionGetOffset(sectionPlasticStrain, setPointID(cell), vecOffset, ierr)) plasticStrainCell = plasticStrainArray(vecOffset + 1:vecOffset + 1 + SIZEOFMEF90_MATS) !! begin ugly hack !! The following won't work until I overload the LinAlg operations to work with the parent classes ! C2 = ATModel%fractureToughness / ATModel%cw * 0.5_kr * ATModel%internalLength * ATModel%toughnessAnisotropyMatrix select type (k => ATModel%toughnessAnisotropyMatrix) type is (MEF90_MATS) C2 = ATModel%fractureToughness / ATModel%cw * 0.5_kr * ATModel%internalLength * k end select !! end ugly hack EEDGaussPlus = 0.0_kr if ((.not. ATModel%isElastic) .and. ((.not. (MEF90DefMechGlobalOptions%multiPhaseField)) .or. (MEF90DefMechCtx%currentSet == set))) then call Split%setup(totalStrainGauss - plasticStrainCell, ierr) call Split%EED(HookesLaw, totalStrainGauss - plasticStrainCell, EEDGaussPlus, EEDGaussMinus, ierr) end if C1 = ATModel%fractureToughness / ATModel%cw * 0.25_kr / ATModel%internalLength C3 = ATModel%D2a(damageGauss) * EEDGaussPlus + C1 * ATModel%D2w(damageGauss) do jDof = 1, numDofDamage do iDof = 1, numDofDamage matDof(jDof, iDof) = matDof(jDof, iDof) + elemScal(cell)%Gauss_C(iGauss) * ( & C3 * elemScal(cell)%BF(iDof, iGauss) * elemScal(cell)%BF(jDof, iGauss) & + (C2 * elemScal(cell)%Grad_BF(iDof, iGauss) .DotP. elemScal(cell)%Grad_BF(jDof, iGauss)) & ) end do ! iDof numDofDamage end do ! jDof numDofDamage end do ! iGauss PetscCall(DMPlexVecRestoreClosure(dmDisplacement, PETSC_NULL_SECTION, MEF90DefMechCtx%displacementLocal, setPointID(cell), PETSC_NULL_INTEGER, displacementDof, ierr)) PetscCall(DMPlexVecRestoreClosure(dmDamage, PETSC_NULL_SECTION, damageLocal, setPointID(cell), PETSC_NULL_INTEGER, damageDof, ierr)) PetscCall(DMPlexVecRestoreClosure(dmTemperature, PETSC_NULL_SECTION, MEF90DefMechCtx%TemperatureLocal, setPointID(cell), PETSC_NULL_INTEGER, temperatureDof, ierr)) PetscCall(DMPlexMatSetClosure(dmDamage, PETSC_NULL_SECTION, PETSC_NULL_SECTION, A, setPointID(cell), reshape(matDof, [numDofDamage*numDofDamage]), ADD_VALUES, ierr)) end do ! cell deallocate (matDof) PetscCall(MEF90ElementDestroy(elemVect, ierr)) PetscCall(MEF90ElementDestroy(elemScal, ierr)) PetscCall(ISRestoreIndices(setPointIS, setPointID, ierr)) PetscCall(ISDestroy(setPointIS, ierr)) end if ! setPointIS end do ! set PetscCall(ISRestoreIndices(setIS, setID, ierr)) PetscCall(VecRestoreArray(MEF90DefMechCtx%plasticStrain, plasticStrainArray, ierr)) end if ! setIS PetscCall(ISDestroy(setIS, ierr)) PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr)) PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr)) PetscCall(MatCopy(A, M, SAME_NONZERO_PATTERN, ierr)) end subroutine MEF90DefMechBilinearFormDamage #undef __FUNCT__ #define __FUNCT__ "MEF90DefMechTAOHessianDamage" !!! author: Alexis Marboeuf (2022, marboeua@mcmaster.ca) !!! author: Blaise Bourdin (2023, bourdin@mcmaster.ca) !!! !!! MEF90DefMechTAOHessianDamage: !!! subroutine MEF90DefMechTAOHessianDamage(taoDamage, damage, A, M, MEF90DefMechCtx, ierr) type(tTao), intent(IN) :: taoDamage type(tVec), intent(IN) :: damage type(tMat), intent(INOUT) :: A, M type(MEF90DefMech_Type), intent(IN) :: MEF90DefMechCtx PetscErrorCode, intent(INOUT) :: ierr type(tSNES) :: dummySNES ! type(tDM) :: dmDamage ! PetscCall(VecGetDM(MEF90DefMechCtx%damageLocal(1), dmDamage, ierr)) ! PetscCall(DMGlobalToLocal(dmDamage, damage, INSERT_VALUES, MEF90DefMechCtx%damageLocal(1), ierr)) PetscCall(MEF90DefMechBilinearFormDamage(dummySNES, damage, A, M, MEF90DefMechCtx, ierr)) end subroutine MEF90DefMechTAOHessianDamage #undef __FUNCT__ #define __FUNCT__ "MEF90DefMechSurfaceEnergy" !!! author: Blaise Bourdin (2014-2020, bourdin@lsu.edu) !!! author: Blaise Bourdin (2022, bourdin@mcmaster.ca) !!! author: Alexis Marboeuf (2022, marboeua@mcmaster.ca) !!! !!! MEF90DefMechSurfaceEnergy: !!! subroutine MEF90DefMechSurfaceEnergy(MEF90DefMechCtx, energy, ierr) type(MEF90DefMech_Type), intent(IN) :: MEF90DefMechCtx PetscReal, dimension(:), pointer :: energy PetscErrorCode, intent(INOUT) :: ierr type(tDM) :: dmDamage PetscReal, dimension(:), pointer :: damageDof type(tIS) :: setIS, setPointIS PetscInt, dimension(:), pointer :: setID, setPointID PetscInt :: set, QuadratureOrder, cell, dim type(MEF90DefMechCellSetOptions_Type) :: cellSetOptions type(MEF90_ELEMENT_SCAL), dimension(:), pointer :: elemScal type(eDMPolytopeType) :: cellGeometryScal type(MEF90ElementType) :: elemScalType type(MEF90CtxGlobalOptions_Type) :: MEF90CtxGlobalOptions type(MEF90DefMechGlobalOptions_Type) :: MEF90DefMechGlobalOptions class(MEF90DefMechAT_Type), allocatable :: ATModel type(MEF90_VECT) :: gradDamageGauss PetscReal :: damageGauss, myEnergy, C1, C2 PetscInt :: iDof, iGauss, numDofDamage, numGauss character(len=MEF90MXSTRLEN) :: prefix PetscCall(MEF90CtxGlobalOptionsSetFromOptions(MEF90DefMechCtx%MEF90Ctx%comm, trim(MEF90DefMechCtx%MEF90Ctx%prefix), MEF90CtxGlobalOptions, ierr)) PetscCall(MEF90DefMechGlobalOptionsSetFromOptions(MEF90DefMechCtx%comm, trim(MEF90DefMechCtx%prefix), MEF90DefMechGlobalOptions, ierr)) PetscCall(VecGetDM(MEF90DefMechCtx%damageLocal, dmDamage, ierr)) PetscCall(DMGetDimension(dmDamage, dim, ierr)) !! get IS for cell sets PetscCall(DMGetLabelIdIS(dmDamage, MEF90CellSetLabelName, setIS, ierr)) PetscCall(MEF90ISAllGatherMerge(MEF90DefMechCtx%comm, setIS, ierr)) energy = 0.0_kr if (.not. PetscObjectIsNull(setIS)) then PetscCall(ISGetIndices(setIS, setID, ierr)) do set = 1, size(setID) myEnergy = 0.0_kr PetscCall(DMGetStratumIS(dmDamage, MEF90CellSetLabelName, setID(set), setPointIS, ierr)) if (.not. PetscObjectIsNull(setPointIS)) then write (prefix, '(A,"cs",I4.4,"_")') trim(MEF90DefMechCtx%prefix), setID(set) PetscCall(MEF90DefMechCellSetOptionsSetFromOptions(MEF90DefMechCtx%comm, prefix, MEF90_DIM, cellSetOptions, ierr)) PetscCall(ISGetIndices(setPointIS, setPointID, ierr)) PetscCall(DMPlexGetCellType(dmDamage, setPointID(1), cellGeometryScal, ierr)) PetscCall(MEF90ElementGetType(MEF90CtxGlobalOptions%elementFamily, MEF90CtxGlobalOptions%elementOrder, cellGeometryScal, elemScalType, ierr)) !! get the ATModel object PetscCall(MEF90DefMechGetATModel(MEF90DefMechCtx%comm, prefix, dim, ATModel, ierr)) PetscCall(ATModel%setFromOptions(ierr)) !! Allocate elements QuadratureOrder = max(ATmodel%wOrder, 2 * (elemScalType%order - 1)) PetscCall(MEF90ElementCreate(dmDamage, setPointIS, elemScal, QuadratureOrder, elemScalType, ierr)) numDofDamage = size(elemScal(1)%BF(:, 1)) numGauss = size(elemScal(1)%Gauss_C) C1 = ATModel%fractureToughness / ATModel%cw * 0.25_kr / ATModel%internalLength C2 = ATModel%fractureToughness / ATModel%cw * 0.25_kr * ATModel%internalLength do cell = 1, size(setPointID) PetscCall(DMPlexVecGetClosure(dmDamage, PETSC_NULL_SECTION, MEF90DefMechCtx%damageLocal, setPointID(cell), PETSC_NULL_INTEGER, damageDof, ierr)) do iGauss = 1, numGauss damageGauss = 0.0_kr gradDamageGauss = 0.0_kr if (.not. ATModel%isElastic) then do iDof = 1, numDofDamage damageGauss = damageGauss + damageDof(iDof) * elemScal(cell)%BF(iDof, iGauss) gradDamageGauss = gradDamageGauss + damageDof(iDof) * elemScal(cell)%Grad_BF(iDof, iGauss) end do ! iDof numDofDamage end if !! UGLY but until I can figure out how to handle polymorphism, this has to be this way: select type (k => ATModel%toughnessAnisotropyMatrix) type is (MEF90_MATS) myEnergy = myEnergy + elemScal(cell)%Gauss_C(iGauss) * & (C1 * ATModel%w(damageGauss) + C2 * ((k * gradDamageGauss) .DotP. gradDamageGauss)) end select end do ! iGauss PetscCall(DMPlexVecRestoreClosure(dmDamage, PETSC_NULL_SECTION, MEF90DefMechCtx%damageLocal, setPointID(cell), PETSC_NULL_INTEGER, damageDof, ierr)) end do ! cell PetscCall(MEF90ElementDestroy(elemScal, ierr)) end if ! setPointIS PetscCall(ISDestroy(setPointIS, ierr)) PetscCallMPI(MPI_AllReduce(myEnergy, energy(set), 1, MPIU_SCALAR, MPI_SUM, MEF90DefMechCtx%comm, ierr)) end do ! set PetscCall(ISRestoreIndices(setIS, setID, ierr)) end if ! setIS PetscCall(ISDestroy(setIS, ierr)) end subroutine MEF90DefMechSurfaceEnergy #undef __FUNCT__ #define __FUNCT__ "MEF90DefMechTAObjectiveDamage" !!! author: Alexis Marboeuf (2022, marboeua@mcmaster.ca) !!! author: Blaise Bourdin (2023, bourdin@mcmaster.ca) !!! !!! MEF90DefMechTAOObjectiveDamage: !!! subroutine MEF90DefMechTAOObjectiveDamage(taoDamage, damage, energy, MEF90DefMechCtx, ierr) type(MEF90DefMech_Type), intent(IN) :: MEF90DefMechCtx type(tTao), intent(IN) :: taoDamage type(tVec), intent(IN) :: damage PetscReal, intent(INOUT) :: energy PetscErrorCode, intent(INOUT) :: ierr PetscReal, dimension(:), pointer :: surfaceEnergy, elasticEnergy type(tDM) :: dmDamage type(tIS) :: setIS PetscInt :: numSet PetscCall(VecGetDM(MEF90DefMechCtx%damageLocal, dmDamage, ierr)) PetscCall(DMGlobalToLocal(dmDamage, damage, INSERT_VALUES, MEF90DefMechCtx%damageLocal, ierr)) !! get IS for cell sets PetscCall(DMGetLabelIdIS(dmDamage, MEF90CellSetLabelName, setIS, ierr)) PetscCall(MEF90ISAllGatherMerge(MEF90DefMechCtx%comm, setIS, ierr)) numSet = 0 if (.not. PetscObjectIsNull(setIS)) then PetscCall(ISGetSize(setIS, numSet, ierr)) PetscCall(ISDestroy(setIS, ierr)) end if allocate (surfaceEnergy(numSet)) allocate (elasticEnergy(numSet)) surfaceEnergy = 0.0_kr elasticEnergy = 0.0_kr PetscCall(MEF90DefMechSurfaceEnergy(MEF90DefMechCtx, surfaceEnergy, ierr)) PetscCall(MEF90DefMechElasticEnergy(MEF90DefMechCtx, elasticEnergy, ierr)) energy = sum(surfaceEnergy) + sum(elasticEnergy) end subroutine MEF90DefMechTAOObjectiveDamage #undef __FUNCT__ #define __FUNCT__ "MEF90DefMechCrackVolume" !!! author: Erwan Tanne (2016-2021, erwan.tanne@gmail.com) !!! author: Blaise Bourdin (2016-2021, bourdin@lsu.edu) !!! author: Blaise Bourdin (2022, bourdin@mcmaster.ca) !!! author: Alexis Marboeuf (2022, marboeua@mcmaster.ca) !!! !!! MEF90DefMechCrackVolume: !!! subroutine MEF90DefMechCrackVolume(MEF90DefMechCtx, CrackVolume, ierr) type(MEF90DefMech_Type), intent(IN) :: MEF90DefMechCtx PetscReal, dimension(:), pointer :: CrackVolume PetscErrorCode, intent(INOUT) :: ierr type(tDM) :: dmDamage, dmDisplacement PetscReal, dimension(:), pointer :: damageDof, displacementDof type(tIS) :: setIS, setPointIS PetscInt, dimension(:), pointer :: setID, setPointID PetscInt :: set, QuadratureOrder, cell, dim type(MEF90DefMechCellSetOptions_Type) :: cellSetOptions type(MEF90_ELEMENT_ELAST), dimension(:), pointer :: elemVect type(MEF90_ELEMENT_SCAL), dimension(:), pointer :: elemScal type(eDMPolytopeType) :: cellGeometryVect type(eDMPolytopeType) :: cellGeometryScal type(MEF90ElementType) :: elemVectType, elemScalType type(MEF90CtxGlobalOptions_Type) :: MEF90CtxGlobalOptions class(MEF90DefMechAT_Type), allocatable :: ATModel type(MEF90_VECT) :: gradDamageGauss, displacementCell PetscReal :: myCrackVolume PetscInt :: iDof, iGauss, numDofDamage, numDofDisplacement, numGauss character(len=MEF90MXSTRLEN) :: prefix PetscCall(MEF90CtxGlobalOptionsSetFromOptions(MEF90DefMechCtx%MEF90Ctx%comm, trim(MEF90DefMechCtx%MEF90Ctx%prefix), MEF90CtxGlobalOptions, ierr)) PetscCall(VecGetDM(MEF90DefMechCtx%damageLocal, dmDamage, ierr)) PetscCall(VecGetDM(MEF90DefMechCtx%displacementLocal, dmDisplacement, ierr)) PetscCall(DMGetDimension(dmDamage, dim, ierr)) !! get IS for cell sets PetscCall(DMGetLabelIdIS(dmDamage, MEF90CellSetLabelName, setIS, ierr)) PetscCall(MEF90ISAllGatherMerge(MEF90DefMechCtx%comm, setIS, ierr)) CrackVolume = 0.0_kr if (.not. PetscObjectIsNull(setIS)) then PetscCall(ISGetIndices(setIS, setID, ierr)) do set = 1, size(setID) myCrackVolume = 0.0_kr PetscCall(DMGetStratumIS(dmDamage, MEF90CellSetLabelName, setID(set), setPointIS, ierr)) if (.not. PetscObjectIsNull(setPointIS)) then write (prefix, '(A,"cs",I4.4,"_")') trim(MEF90DefMechCtx%prefix), setID(set) PetscCall(MEF90DefMechCellSetOptionsSetFromOptions(MEF90DefMechCtx%comm, prefix, MEF90_DIM, cellSetOptions, ierr)) PetscCall(ISGetIndices(setPointIS, setPointID, ierr)) PetscCall(DMPlexGetCellType(dmDisplacement, setPointID(1), cellGeometryVect, ierr)) PetscCall(DMPlexGetCellType(dmDamage, setPointID(1), cellGeometryScal, ierr)) PetscCall(MEF90ElementGetType(MEF90CtxGlobalOptions%elementFamily, MEF90CtxGlobalOptions%elementOrder, cellGeometryVect, elemVectType, ierr)) PetscCall(MEF90ElementGetType(MEF90CtxGlobalOptions%elementFamily, MEF90CtxGlobalOptions%elementOrder, cellGeometryScal, elemScalType, ierr)) !! get the ATModel and split objects PetscCall(MEF90DefMechGetATModel(MEF90DefMechCtx%comm, prefix, dim, ATModel, ierr)) PetscCall(ATModel%setFromOptions(ierr)) !! Allocate elements QuadratureOrder = 2 * (elemScalType%order - 1) PetscCall(MEF90ElementCreate(dmDisplacement, setPointIS, elemVect, QuadratureOrder, elemVectType, ierr)) PetscCall(MEF90ElementCreate(dmDamage, setPointIS, elemScal, QuadratureOrder, elemScalType, ierr)) numDofDisplacement = size(elemVect(1)%BF(:, 1)) numDofDamage = size(elemScal(1)%BF(:, 1)) numGauss = size(elemScal(1)%Gauss_C) do cell = 1, size(setPointID) PetscCall(DMPlexVecGetClosure(dmDamage, PETSC_NULL_SECTION, MEF90DefMechCtx%damageLocal, setPointID(cell), PETSC_NULL_INTEGER, damageDof, ierr)) PetscCall(DMPlexVecGetClosure(dmDisplacement, PETSC_NULL_SECTION, MEF90DefMechCtx%displacementLocal, setPointID(cell), PETSC_NULL_INTEGER, displacementDof, ierr)) do iGauss = 1, numGauss gradDamageGauss = 0.0_kr if (.not. ATModel%isElastic) then do iDof = 1, numDofDamage gradDamageGauss = gradDamageGauss + damageDof(iDof) * elemScal(cell)%Grad_BF(iDof, iGauss) end do ! iDof numDofDamage end if displacementCell = 0.0_kr do iDof = 1, numDofDisplacement displacementCell = displacementCell + displacementDof(iDof) * elemVect(cell)%BF(iDof, iGauss) end do ! iDof numDofDisplacement myCrackVolume = myCrackVolume + elemScal(cell)%Gauss_C(iGauss) * ( & gradDamageGauss .DotP. displacementCell) end do ! iGauss PetscCall(DMPlexVecRestoreClosure(dmDamage, PETSC_NULL_SECTION, MEF90DefMechCtx%damageLocal, setPointID(cell), PETSC_NULL_INTEGER, damageDof, ierr)) PetscCall(DMPlexVecRestoreClosure(dmDisplacement, PETSC_NULL_SECTION, MEF90DefMechCtx%displacementLocal, setPointID(cell), PETSC_NULL_INTEGER, displacementDof, ierr)) end do ! cell PetscCall(MEF90ElementDestroy(elemScal, ierr)) PetscCall(MEF90ElementDestroy(elemVect, ierr)) PetscCall(ISDestroy(setPointIS, ierr)) end if ! setPointIS PetscCallMPI(MPI_AllReduce(myCrackVolume, CrackVolume(set), 1, MPIU_SCALAR, MPI_SUM, MEF90DefMechCtx%comm, ierr)) end do ! set PetscCall(ISRestoreIndices(setIS, setID, ierr)) PetscCall(ISDestroy(setIS, ierr)) end if ! setIS end subroutine MEF90DefMechCrackVolume end module MEF90_APPEND(m_MEF90_DefMechAssembly,MEF90_DIM)D