module m_MEF90_DMPlex #include "petsc/finclude/petsc.h" #include "petsc/finclude/petscsf.h" use m_MEF90_Parameters use m_MEF90_Utils use m_MEF90_LinAlg use m_MEF90_Elements use m_MEF90_Ctx use, intrinsic :: iso_c_binding #ifdef MEF90_HAVE_SYMENGINEF90 use symengine #endif implicit none(type) enum, bind(c) enumerator :: MEF90CellSetType = 1, & MEF90FaceSetType, & !MEF90EdgeSetType, & MEF90VertexSetType end enum !!! Not sure why this has to be made explicitly public PetscEnum, dimension(3), parameter, public :: MEF90SetType = [ & MEF90CellSetType, & MEF90FaceSetType, & !MEF90EdgeSetType, & MEF90VertexSetType] character(len=MEF90MXSTRLEN), parameter :: MEF90CellSetLabelName = 'Cell Sets ' character(len=MEF90MXSTRLEN), parameter :: MEF90FaceSetLabelName = 'Face Sets ' !Character(len=MEF90MXSTRLEN),Parameter :: MEF90EdgeSetLabelName = 'Edge Sets ' character(len=MEF90MXSTRLEN), parameter :: MEF90VertexSetLabelName = 'Vertex Sets' character(len=MEF90MXSTRLEN), dimension(3), parameter :: MEF90SetLabelName = & [MEF90CellSetLabelName, & MEF90FaceSetLabelName, & !MEF90EdgeSetLabelName, & MEF90VertexSetLabelName] character(len=MEF90MXSTRLEN), parameter :: MEF90CellSetprefix = 'cs' character(len=MEF90MXSTRLEN), parameter :: MEF90FaceSetprefix = 'fs' !Character(len=MEF90MXSTRLEN),Parameter :: MEF90EdgeSetprefix = 'es' character(len=MEF90MXSTRLEN), parameter :: MEF90VertexSetprefix = 'vs' character(len=MEF90MXSTRLEN), dimension(3), parameter :: MEF90SetPrefix = & [MEF90CellSetPrefix, & MEF90FaceSetPrefix, & !MEF90EdgeSetPrefix, & MEF90VertexSetPrefix] private public :: MEF90CellSetLabelName, & MEF90FaceSetLabelName, & !MEF90EdgeSetLabelName, & MEF90VertexSetLabelName, & MEF90SetLabelName, & MEF90CellSetPrefix, & MEF90FaceSetPrefix, & !MEF90EdgeSetPrefix, & MEF90VertexSetPrefix, & MEF90SetPrefix, & MEF90CellSetType, & MEF90FaceSetType, & !MEF90EdgeSetType, & MEF90VertexSetType, & MEF90SectionAllocateDof, MEF90SectionAllocateDofSet, & MEF90SetupConstraintTableSet, MEF90SectionAllocateConstraint, & MEF90CellSectionCreate, & MEF90VecCopySF, MEF90IOSFCreate, MEF90FaceSetIOSFCreate, & MEF90ConstraintSFCreate, MEF90VecGlobalToLocalConstraint, & MEF90VecCreateIO, & MEF90CreateLocalVector, & MEF90CreateBoundaryLocalVector, & MEF90CreateCellVector, & MEF90CreateBoundaryCellVector, & MEF90VecSetBCValuesFromOptions, & MEF90VecSetValuesFromOptions, & MEF90VecSetBCValuesFromOptionsExpr, & MEF90VecSetValuesFromOptionsExpr, & MEF90VecGetClosureSize, & MEF90DMGetNumSets contains #undef __FUNCT__ #define __FUNCT__ "MEF90DMGetNumSets" !!! author: Blaise Bourdin (2026, bourdin@mcmaster.ca) !!! !!! MEF90DMGetNumSets: the total number of sets of a given label, over the whole communicator of dm !!! and not only the ones with a point on the local rank !!! subroutine MEF90DMGetNumSets(dm, labelName, numSet, ierr) type(tDM), intent(IN) :: dm character(len=*), intent(IN) :: labelName PetscInt, intent(OUT) :: numSet PetscErrorCode, intent(INOUT) :: ierr type(tIS) :: setIS MPIU_Comm :: comm PetscCall(PetscObjectGetComm(dm, comm, ierr)) PetscCall(DMGetLabelIdIS(dm, labelName, setIS, ierr)) PetscCall(MEF90ISAllGatherMerge(comm, setIS, ierr)) PetscCall(ISGetLocalSize(setIS, numSet, ierr)) PetscCall(ISDestroy(setIS, ierr)) end subroutine MEF90DMGetNumSets #undef __FUNCT__ #define __FUNCT__ "MEF90CreateLocalVector" !!! author: Blaise Bourdin (2022, bourdin@mcmaster.ca) !!! !!! MEF90CreateLocalVector: create a Vec associated with a FE space and constraints !!! cell set BC are obtained from the command line option -cs_BC [bool], [bool], ... !!! face set BC are obtained from the command line option -fs_BC [bool], [bool], ... !!! vertex set BC are obtained from the command line option -vs_BC [bool], [bool], ... !!! subroutine MEF90CreateLocalVector(dm, elemFamily, elemOrder, sDim, name, V, ierr) type(tDM), intent(IN) :: dm PetscEnum, intent(IN) :: elemFamily PetscInt, intent(IN) :: elemOrder, sDim character(len=MEF90MXSTRLEN), intent(IN) :: name type(tVec), intent(OUT) :: V PetscErrorCode, intent(INOUT) :: ierr type(tPetscSection) :: sectionV type(tDM) :: dmV PetscInt :: fieldV = 0 PetscInt :: pStart, pEnd MPIU_Comm :: comm PetscInt :: set PetscEnum :: setType PetscInt, dimension(:), pointer :: setID, pointID type(tIS) :: setIS, pointIS type(MEF90ElementType) :: elemType DMPolytopeType :: cellType PetscBool, dimension(:, :), pointer :: constraintTruthTable PetscBool, dimension(:), pointer :: setConstraints PetscInt :: numBC PetscBool :: flg character(len=MEF90MXSTRLEN) :: BCOptionName type(tPetscSF) :: naturalPointSF, naturalSF PetscCall(DMClone(dm, dmV, ierr)) PetscCall(PetscObjectSetName(dmv, name, ierr)) PetscCall(DMGetUseNatural(dm, flg, ierr)) PetscCall(DMSetUseNatural(dmV, flg, ierr)) PetscCall(PetscObjectGetComm(dmV, comm, ierr)) PetscCall(PetscSectionCreate(comm, sectionV, ierr)) PetscCall(PetscObjectSetName(sectionV, name, ierr)) PetscCall(PetscSectionSetNumFields(sectionV, 1_ki, ierr)) PetscCall(PetscSectionSetFieldName(sectionV, fieldV, trim(name), ierr)) PetscCall(PetscSectionSetFieldComponents(sectionV, fieldV, sdim, ierr)) PetscCall(DMPlexGetChart(dmV, pStart, pEnd, ierr)) PetscCall(PetscSectionSetChart(sectionV, pStart, pEnd, ierr)) PetscCall(DMGetLabelIdIS(dmV, MEF90CellSetLabelName, setIS, ierr)) !! Get a GLOBAL cell set IS ! PetscCall(MEF90ISAllGatherMerge(comm,setIS,ierr)) if (.not. PetscObjectIsNull(setIS)) then PetscCall(ISGetIndices(setIS, setID, ierr)) do set = 1, size(setID) !! Get cell type in order to pick the proper element type. !! We assume that all cells in a set have the same type, so all we need it to query the first cell in the set PetscCall(DMGetStratumIS(dmV, MEF90CellSetLabelName, setID(set), pointIS, ierr)) if (.not. PetscObjectIsNull(pointIS)) then PetscCall(ISGetIndices(pointIS, pointID, ierr)) PetscCall(DMPlexGetCellType(dmV, pointID(1), cellType, ierr)) PetscCall(MEF90ElementGetType(elemFamily, elemOrder, cellType, elemType, ierr)) PetscCall(MEF90SectionAllocateDofSet(dmV, MEF90CellSetType, setID(set), elemType, sdim, sectionV, ierr)) PetscCall(ISRestoreIndices(pointIS, pointID, ierr)) end if ! pointIS PetscCall(ISDestroy(pointIS, ierr)) end do ! set PetscCall(ISRestoreIndices(setIS, setID, ierr)) end if ! setIS PetscCall(ISDestroy(setIS, ierr)) allocate (ConstraintTruthTable(pEnd, sDim)) ConstraintTruthTable = .false. allocate (setConstraints(sDim)) do setType = 1, size(MEF90SetType) PetscCall(DMGetLabelIdIS(dm, MEF90SetLabelName(setType), setIS, ierr)) ! PetscCall(MEF90ISAllGatherMerge(comm,setIS,ierr)) if (.not. PetscObjectIsNull(setIS)) then PetscCall(ISGetIndices(setIS, setID, ierr)) do set = 1, size(setID) setConstraints = .false. write (BCOptionName, '("-",a2,I4.4,"_",a,"BC")') MEF90SetPrefix(setType), setID(set), trim(name) numBC = sDim PetscCall(PetscOptionsGetBoolArray(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, trim(BCOptionName), setConstraints, numBC, flg, ierr)) PetscCall(MEF90SetupConstraintTableSet(dmV, sectionV, MEF90SetType(setType), setID(set), setConstraints, ConstraintTruthTable, ierr)) end do PetscCall(ISRestoreIndices(setIS, setID, ierr)) end if ! setIS PetscCall(ISDestroy(setIS, ierr)) end do ! setType PetscCall(MEF90SectionAllocateConstraint(dmV, ConstraintTruthTable, sectionV, ierr)) deallocate (ConstraintTruthTable) deallocate (setConstraints) PetscCall(DMSetLocalSection(dmV, sectionV, ierr)) PetscCall(DMPlexGetMigrationSF(dm, naturalPointSF, ierr)) if (.not. PetscObjectIsNull(naturalPointSF)) then PetscCall(DMPlexSetMigrationSF(dmV, naturalPointSF, ierr)) PetscCall(DMPlexCreateGlobalToNaturalSF(dmV, PETSC_NULL_SECTION, naturalPointSF, naturalSF, ierr)) PetscCall(DMSetNaturalSF(dmV, naturalSF, ierr)) end if #ifdef PETSC_USE_DEBUG write (BCOptionName, '("-",a,"_section_view")') trim(name) PetscCall(PetscSectionViewFromOptions(sectionV, PETSC_NULL_OBJECT, BCOptionName, ierr)) if (.not. PetscObjectIsNull(naturalSF)) then write (BCOptionName, '("-",a,"_naturalSF_view")') trim(name) PetscCall(PetscSFViewFromOptions(naturalSF, PETSC_NULL_OBJECT, BCOptionName, ierr)) end if #endif PetscCall(DMCreateLocalVector(dmV, V, ierr)) PetscCall(PetscObjectSetName(V, name, ierr)) PetscCall(PetscSectionDestroy(sectionV, ierr)) if (.not. PetscObjectIsNull(naturalSF)) then PetscCall(PetscSFDestroy(naturalSF, ierr)) end if PetscCall(DMDestroy(dmV, ierr)) end subroutine MEF90CreateLocalVector #undef __FUNCT__ #define __FUNCT__ "MEF90CreateCellVector" !!! author: Blaise Bourdin (2022, bourdin@mcmaster.ca) !!! !!! MEF90CreateCellVector: create a Vec for a cell-based vector !!! subroutine MEF90CreateCellVector(dm, sDim, name, V, ierr) type(tDM), intent(IN) :: dm PetscInt, intent(IN) :: sDim character(len=MEF90MXSTRLEN), intent(IN) :: name type(tVec), intent(OUT) :: V PetscErrorCode, intent(INOUT) :: ierr PetscInt :: dim type(tPetscSection) :: sectionV type(tDM) :: dmV PetscInt :: fieldV = 0 PetscInt :: pStart, pEnd MPIU_Comm :: comm PetscInt :: set PetscInt, dimension(:), pointer :: setID type(tIS) :: setIS type(MEF90ElementType) :: elemType PetscBool :: flg type(tPetscSF) :: naturalPointSF, naturalSF PetscCall(DMClone(dm, dmV, ierr)) PetscCall(PetscObjectSetName(dmv, name, ierr)) PetscCall(DMGetUseNatural(dm, flg, ierr)) PetscCall(DMSetUseNatural(dmV, flg, ierr)) PetscCall(PetscObjectGetComm(dmV, comm, ierr)) PetscCall(PetscSectionCreate(comm, sectionV, ierr)) PetscCall(PetscObjectSetName(sectionV, name, ierr)) PetscCall(PetscSectionSetNumFields(sectionV, 1_ki, ierr)) PetscCall(PetscSectionSetFieldName(sectionV, fieldV, trim(name), ierr)) PetscCall(PetscSectionSetFieldComponents(sectionV, fieldV, sdim, ierr)) PetscCall(DMPlexGetChart(dmV, pStart, pEnd, ierr)) PetscCall(PetscSectionSetChart(sectionV, pStart, pEnd, ierr)) PetscCall(DMGetDimension(dm, dim, ierr)) if (dim == 2) then elemType = MEF90P0Lagrange2D else elemType = MEF90P0Lagrange3D end if PetscCall(DMGetLabelIdIS(dmV, MEF90CellSetLabelName, setIS, ierr)) !! Get a GLOBAL cell set IS ! PetscCall(MEF90ISAllGatherMerge(comm,setIS,ierr)) if (.not. PetscObjectIsNull(setIS)) then PetscCall(ISGetIndices(setIS, setID, ierr)) do set = 1, size(setID) PetscCall(MEF90SectionAllocateDofSet(dmV, MEF90CellSetType, setID(set), elemType, sdim, sectionV, ierr)) end do ! set PetscCall(ISRestoreIndices(setIS, setID, ierr)) end if ! setIS PetscCall(ISDestroy(setIS, ierr)) PetscCall(PetscSectionSetup(sectionV, ierr)) PetscCall(DMSetLocalSection(dmV, sectionV, ierr)) PetscCall(DMPlexGetMigrationSF(dm, naturalPointSF, ierr)) if (.not. PetscObjectIsNull(naturalPointSF)) then PetscCall(DMPlexSetMigrationSF(dmV, naturalPointSF, ierr)) PetscCall(DMPlexCreateGlobalToNaturalSF(dmV, PETSC_NULL_SECTION, naturalPointSF, naturalSF, ierr)) PetscCall(DMSetNaturalSF(dmV, naturalSF, ierr)) end if #ifdef PETSC_USE_DEBUG debugBlock: block character(len=MEF90MXSTRLEN) :: BCoptionName write (BCOptionName, '("-",a,"_section_view")') trim(name) PetscCall(PetscSectionViewFromOptions(sectionV, PETSC_NULL_OBJECT, BCOptionName, ierr)) if (.not. PetscObjectIsNull(naturalSF)) then write (BCOptionName, '("-",a,"_naturalSF_view")') trim(name) PetscCall(PetscSFViewFromOptions(naturalSF, PETSC_NULL_OBJECT, BCOptionName, ierr)) end if end block debugBlock #endif PetscCall(DMCreateLocalVector(dmV, V, ierr)) PetscCall(PetscObjectSetName(V, name, ierr)) PetscCall(DMDestroy(dmV, ierr)) PetscCall(PetscSectionDestroy(sectionV, ierr)) if (.not. PetscObjectIsNull(naturalSF)) then PetscCall(PetscSFDestroy(naturalSF, ierr)) end if end subroutine MEF90CreateCellVector #undef __FUNCT__ #define __FUNCT__ "MEF90CreateBoundaryLocalVector" !!! author: Blaise Bourdin (2022, bourdin@mcmaster.ca) !!! !!! MEF90CreateBoundaryLocalVector: create a Vec associated with a FE space over the boundary of a domain (face sets) !!! and constraints !!! face set BC are obtained from the command line option -fs_BC [bool], [bool], ... !!! vertex set BC are obtained from the command line option -vs_BC [bool], [bool], ... !!! subroutine MEF90CreateBoundaryLocalVector(dm, elemFamily, elemOrder, sDim, name, V, ierr) type(tDM), intent(IN) :: dm PetscEnum, intent(IN) :: elemFamily PetscInt, intent(IN) :: elemOrder, sDim character(len=MEF90MXSTRLEN), intent(IN) :: name type(tVec), intent(OUT) :: V PetscErrorCode, intent(INOUT) :: ierr type(tPetscSection) :: sectionV type(tDM) :: dmV PetscInt :: fieldV = 0 PetscInt :: pStart, pEnd MPIU_Comm :: comm PetscInt :: set PetscEnum :: setType PetscInt, dimension(:), pointer :: setID, pointID type(tIS) :: setIS, pointIS type(MEF90ElementType) :: elemType DMPolytopeType :: cellType PetscBool, dimension(:, :), pointer :: constraintTruthTable PetscBool, dimension(:), pointer :: setConstraints PetscInt :: numBC PetscBool :: flg character(len=MEF90MXSTRLEN) :: BCOptionName PetscCall(DMClone(dm, dmV, ierr)) PetscCall(PetscObjectSetName(dmv, name, ierr)) PetscCall(DMGetUseNatural(dm, flg, ierr)) PetscCall(DMSetUseNatural(dmV, flg, ierr)) PetscCall(PetscObjectGetComm(dmV, comm, ierr)) PetscCall(PetscSectionCreate(comm, sectionV, ierr)) PetscCall(PetscObjectSetName(sectionV, name, ierr)) PetscCall(PetscSectionSetNumFields(sectionV, 1_ki, ierr)) PetscCall(PetscSectionSetFieldName(sectionV, fieldV, trim(name), ierr)) PetscCall(PetscSectionSetFieldComponents(sectionV, fieldV, sdim, ierr)) PetscCall(DMPlexGetChart(dmV, pStart, pEnd, ierr)) PetscCall(PetscSectionSetChart(sectionV, pStart, pEnd, ierr)) PetscCall(DMGetLabelIdIS(dmV, MEF90FaceSetLabelName, setIS, ierr)) !! Get a GLOBAL cell set IS ! PetscCall(MEF90ISAllGatherMerge(comm,setIS,ierr)) if (.not. PetscObjectIsNull(setIS)) then PetscCall(ISGetIndices(setIS, setID, ierr)) do set = 1, size(setID) !! Get cell type in order to pick the proper element type. !! We assume that all cells in a set have the same type, so all we need it to query the first cell in the set PetscCall(DMGetStratumIS(dmV, MEF90FaceSetLabelName, setID(set), pointIS, ierr)) if (.not. PetscObjectIsNull(pointIS)) then PetscCall(ISGetIndices(pointIS, pointID, ierr)) PetscCall(DMPlexGetCellType(dmV, pointID(1), cellType, ierr)) PetscCall(MEF90ElementGetTypeBoundary(elemFamily, elemOrder, cellType, elemType, ierr)) PetscCall(MEF90SectionAllocateDofSet(dmV, MEF90FaceSetType, setID(set), elemType, sdim, sectionV, ierr)) PetscCall(ISRestoreIndices(pointIS, pointID, ierr)) end if ! pointIS PetscCall(ISDestroy(pointIS, ierr)) end do ! set PetscCall(ISRestoreIndices(setIS, setID, ierr)) end if ! setIS PetscCall(ISDestroy(setIS, ierr)) allocate (ConstraintTruthTable(pEnd, sDim)) ConstraintTruthTable = .false. allocate (setConstraints(sDim)) do setType = 1, size(MEF90SetType) PetscCall(DMGetLabelIdIS(dm, MEF90SetLabelName(setType), setIS, ierr)) ! PetscCall(MEF90ISAllGatherMerge(comm,setIS,ierr)) if (.not. PetscObjectIsNull(setIS)) then PetscCall(ISGetIndices(setIS, setID, ierr)) do set = 1, size(setID) setConstraints = .false. write (BCOptionName, '("-",a2,I4.4,"_",a,"BC")') MEF90SetPrefix(setType), setID(set), trim(name) numBC = sDim PetscCall(PetscOptionsGetBoolArray(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, trim(BCOptionName), setConstraints, numBC, flg, ierr)) PetscCall(MEF90SetupConstraintTableSet(dmV, sectionV, MEF90SetType(setType), setID(set), setConstraints, ConstraintTruthTable, ierr)) end do PetscCall(ISRestoreIndices(setIS, setID, ierr)) end if ! setIS PetscCall(ISDestroy(setIS, ierr)) end do ! setType PetscCall(MEF90SectionAllocateConstraint(dmV, ConstraintTruthTable, sectionV, ierr)) deallocate (ConstraintTruthTable) deallocate (setConstraints) PetscCall(DMSetLocalSection(dmV, sectionV, ierr)) #ifdef PETSC_USE_DEBUG debugBlock: block character(len=MEF90MXSTRLEN) :: BCoptionName write (BCOptionName, '("-",a,"_section_view")') trim(name) PetscCall(PetscSectionViewFromOptions(sectionV, PETSC_NULL_OBJECT, BCOptionName, ierr)) end block debugBlock #endif PetscCall(DMCreateLocalVector(dmV, V, ierr)) PetscCall(PetscObjectSetName(V, name, ierr)) PetscCall(PetscSectionDestroy(sectionV, ierr)) PetscCall(DMDestroy(dmV, ierr)) end subroutine MEF90CreateBoundaryLocalVector #undef __FUNCT__ #define __FUNCT__ "MEF90CreateBoundaryCellVector" !!! author: Blaise Bourdin (2022, bourdin@mcmaster.ca) !!! !!! MEF90CreateBoundaryCellVector: create a Vec for a cell-based vector !!! subroutine MEF90CreateBoundaryCellVector(dm, sDim, name, V, ierr) type(tDM), intent(IN) :: dm PetscInt, intent(IN) :: sDim character(len=MEF90MXSTRLEN), intent(IN) :: name type(tVec), intent(OUT) :: V PetscErrorCode, intent(INOUT) :: ierr PetscInt :: dim type(tPetscSection) :: sectionV type(tDM) :: dmV PetscInt :: fieldV = 0 PetscInt :: pStart, pEnd MPIU_Comm :: comm PetscInt :: set PetscInt, dimension(:), pointer :: setID type(tIS) :: setIS type(MEF90ElementType) :: elemType PetscBool :: flg PetscCall(DMClone(dm, dmV, ierr)) PetscCall(PetscObjectSetName(dmv, name, ierr)) PetscCall(DMGetUseNatural(dm, flg, ierr)) PetscCall(DMSetUseNatural(dmV, flg, ierr)) PetscCall(PetscObjectGetComm(dmV, comm, ierr)) PetscCall(PetscSectionCreate(comm, sectionV, ierr)) PetscCall(PetscObjectSetName(sectionV, name, ierr)) PetscCall(PetscSectionSetNumFields(sectionV, 1_ki, ierr)) PetscCall(PetscSectionSetFieldName(sectionV, fieldV, trim(name), ierr)) PetscCall(PetscSectionSetFieldComponents(sectionV, fieldV, sdim, ierr)) PetscCall(DMPlexGetChart(dmV, pStart, pEnd, ierr)) PetscCall(PetscSectionSetChart(sectionV, pStart, pEnd, ierr)) PetscCall(DMGetDimension(dm, dim, ierr)) if (dim == 2) then elemType = MEF90P0Lagrange2DBoundary else elemType = MEF90P0Lagrange3DBoundary end if PetscCall(DMGetLabelIdIS(dmV, MEF90FaceSetLabelName, setIS, ierr)) !! Get a GLOBAL cell set IS ! PetscCall(MEF90ISAllGatherMerge(comm,setIS,ierr)) if (.not. PetscObjectIsNull(setIS)) then PetscCall(ISGetIndices(setIS, setID, ierr)) do set = 1, size(setID) PetscCall(MEF90SectionAllocateDofSet(dmV, MEF90FaceSetType, setID(set), elemType, sdim, sectionV, ierr)) end do ! set PetscCall(ISRestoreIndices(setIS, setID, ierr)) end if ! setIS PetscCall(ISDestroy(setIS, ierr)) PetscCall(PetscSectionSetup(sectionV, ierr)) PetscCall(DMSetLocalSection(dmV, sectionV, ierr)) #ifdef PETSC_USE_DEBUG debugBlock: block character(len=MEF90MXSTRLEN) :: BCoptionName write (BCOptionName, '("-",a,"_section_view")') trim(name) PetscCall(PetscSectionViewFromOptions(sectionV, PETSC_NULL_OBJECT, BCOptionName, ierr)) end block debugBlock #endif PetscCall(DMCreateLocalVector(dmV, V, ierr)) PetscCall(PetscObjectSetName(V, name, ierr)) PetscCall(PetscSectionDestroy(sectionV, ierr)) PetscCall(DMDestroy(dmV, ierr)) end subroutine MEF90CreateBoundaryCellVector #undef __FUNCT__ #define __FUNCT__ "MEF90VecGetClosureSize" !!! author: Blaise Bourdin (2022, bourdin@mcmaster.ca) !!! !!! MEF90VecGetClosureSize: Associates dof to a section !!! subroutine MEF90VecGetClosureSize(v, p, clSize, ierr) type(tVec), intent(IN) :: v PetscInt, intent(IN) :: p PetscInt, intent(OUT) :: clSize PetscErrorCode, intent(INOUT) :: ierr type(tDM) :: dm type(tPetscSection) :: section PetscInt, dimension(:), pointer :: closure PetscInt :: point, numDof, numClosure clSize = 0 PetscCall(VecGetDM(v, dm, ierr)) PetscCall(DMGetLocalSection(dm, section, ierr)) PetscCall(DMPlexGetTransitiveClosure(dm, p, PETSC_TRUE, numClosure, closure, ierr)) if (size(closure) > 0) then do point = 1, size(closure), 2 PetscCall(PetscSectionGetDof(section, closure(point), numDof, ierr)) clSize = clSize + numDof end do end if PetscCall(DMPlexRestoreTransitiveClosure(dm, p, PETSC_TRUE, numClosure, closure, ierr)) end subroutine MEF90VecGetClosureSize #undef __FUNCT__ #define __FUNCT__ "MEF90SectionAllocateDof" !!! author: Blaise Bourdin (2022, bourdin@mcmaster.ca) !!! !!! MEF90SectionAllocateDof: Associates dof to a section !!! subroutine MEF90SectionAllocateDof(dm, setType, elemType, numComponents, section, ierr) type(tDM), intent(IN) :: dm PetscEnum, intent(IN) :: setType type(MEF90ElementType), intent(IN) :: elemType PetscInt, intent(IN) :: numComponents type(tPetscSection), intent(INOUT) :: section PetscErrorCode, intent(INOUT) :: ierr type(tIS) :: setIS PetscInt, dimension(:), pointer :: setID PetscInt :: set PetscCall(DMGetLabelIdIS(dm, MEF90SetLabelName(setType), setIS, ierr)) if (.not. PetscObjectIsNull(setIS)) then PetscCall(ISGetIndices(setIS, setID, ierr)) do set = 1, size(setID) PetscCall(MEF90SectionAllocateDofSet(dm, setType, setID(set), elemType, numComponents, Section, ierr)) end do PetscCall(ISRestoreIndices(setIS, setID, ierr)) end if ! setIS PetscCall(ISDestroy(setIS, ierr)) end subroutine MEF90SectionAllocateDof #undef __FUNCT__ #define __FUNCT__ "MEF90_SectionAllocateDofSet" !!! author: Blaise Bourdin (2022, bourdin@mcmaster.ca) !!! !!! MEF90SectionAllocateDofSet: Associates dof to a section in a set !!! subroutine MEF90SectionAllocateDofSet(dm, setType, setID, elemType, numComponents, section, ierr) type(tDM), intent(IN) :: dm PetscEnum, intent(IN) :: setType PetscInt, intent(IN) :: setID type(MEF90ElementType), intent(IN) :: elemType PetscInt, intent(IN) :: numComponents type(tPetscSection), intent(INOUT) :: section PetscErrorCode, intent(INOUT) :: ierr type(tIS) :: setPointIS PetscInt, dimension(:), pointer :: setPointID PetscInt :: point PetscInt, dimension(:), pointer :: closure PetscInt :: p, depth PetscInt :: field = 0 PetscCall(DMGetStratumIS(dm, MEF90SetLabelName(setType), setID, setPointIS, ierr)) if (.not. PetscObjectIsNull(setPointIS)) then PetscCall(ISGetIndices(setPointIS, setPointID, ierr)) !! This can probably be optimized by allocating closure outside of the loop !! But I can't figure out how it is done at the moment. nullify (closure) do point = 1, size(setPointID) PetscCall(DMPlexGetTransitiveClosure(dm, setPointID(point), PETSC_TRUE, PETSC_NULL_INTEGER, closure, ierr)) do p = 1, size(closure), 2 PetscCall(DMPlexGetPointDepth(dm, closure(p), depth, ierr)) if (elemType%numDofs(depth + 1) > 0) then PetscCall(PetscSectionSetDof(section, closure(p), elemType%numDofs(depth + 1) * numComponents, ierr)) PetscCall(PetscSectionSetFieldDof(section, closure(p), field, elemType%numDofs(depth + 1) * numComponents, ierr)) end if end do! p PetscCall(DMPlexRestoreTransitiveClosure(dm, setPointID(point), PETSC_TRUE, PETSC_NULL_INTEGER, closure, ierr)) end do! cell PetscCall(ISRestoreIndices(setPointIS, setPointID, ierr)) end if ! setPointIS PetscCall(ISDestroy(setPointIS, ierr)) end subroutine MEF90SectionAllocateDofSet #undef __FUNCT__ #define __FUNCT__ "MEF90CellSectionCreate" !!! author: Alexis Marboeuf (2022, marboeua@mcmaster.ca) !!! !!! MEF90CellSectionCreate: create a section with numComponent dof at cells !!! subroutine MEF90CellSectionCreate(dm, numComponents, section, ierr) type(tDM), intent(IN) :: dm PetscInt, intent(IN) :: numComponents type(tPetscSection), intent(INOUT) :: section PetscErrorCode, intent(INOUT) :: ierr type(tIS) :: setIS, setPointIS PetscInt, dimension(:), pointer :: setID, setPointID PetscInt :: set, point PetscInt, dimension(:), pointer :: closure PetscInt :: p, depth, dim PetscCall(DMGetDimension(dm, dim, ierr)) PetscCall(DMGetLabelIdIS(dm, 'Cell Sets', setIS, ierr)) if (.not. PetscObjectIsNull(setIS)) then PetscCall(ISGetIndices(setIS, setID, ierr)) do set = 1, size(setID) PetscCall(DMGetStratumIS(dm, 'Cell Sets', setID(set), setPointIS, ierr)) if (.not. PetscObjectIsNull(setPointIS)) then PetscCall(ISGetIndices(setPointIS, setPointID, ierr)) !! This can probably be optimized by allocating closure outside of the loop !! But I can't figure out how it is done at the moment. nullify (closure) do point = 1, size(setPointID) PetscCall(DMPlexGetTransitiveClosure(dm, setPointID(point), PETSC_TRUE, PETSC_NULL_INTEGER, closure, ierr)) do p = 1, size(closure), 2 PetscCall(DMPlexGetPointDepth(dm, closure(p), depth, ierr)) if (depth == dim) then PetscCall(PetscSectionSetDof(section, closure(p), numComponents, ierr)) end if end do! p PetscCall(DMPlexRestoreTransitiveClosure(dm, setPointID(point), PETSC_TRUE, PETSC_NULL_INTEGER, closure, ierr)) end do! cell PetscCall(ISRestoreIndices(setPointIS, setPointID, ierr)) end if ! setPointIS PetscCall(ISDestroy(setPointIS, ierr)) end do ! set PetscCall(ISRestoreIndices(setIS, setID, ierr)) end if ! setIS PetscCall(ISDestroy(setIS, ierr)) PetscCall(PetscSectionSetup(section, ierr)) end subroutine MEF90CellSectionCreate #undef __FUNCT__ #define __FUNCT__ "MEF90SetupConstraintTableSet" !!! author: Blaise Bourdin (2022, bourdin@mcmaster.ca) !!! !!! MEF90SetupConstraintTableSet: Build the contribution of a set to the constraint table !!! subroutine MEF90SetupConstraintTableSet(dm, section, setType, setID, constraints, table, ierr) type(tDM), intent(IN) :: dm type(tPetscSection), intent(INOUT) :: section PetscEnum, intent(IN) :: setType PetscInt, intent(IN) :: setID PetscBool, dimension(:), pointer :: constraints PetscBool, dimension(:, :), pointer :: table PetscErrorCode, intent(INOUT) :: ierr type(tIS) :: setPointIS PetscInt, dimension(:), pointer :: setPointID PetscInt :: point, numDof PetscInt, dimension(:), pointer :: closure PetscInt :: p PetscCall(DMGetStratumIS(dm, MEF90SetLabelName(setType), setID, setPointIS, ierr)) if (.not. PetscObjectIsNull(setPOintIS)) then PetscCall(ISGetIndices(setPointIS, setPointID, ierr)) !! This can probably be optimized by allocating closure outside of the loop !! But I can't figure out how it is done at the moment. nullify (closure) do point = 1, size(setPointID) PetscCall(DMPlexGetTransitiveClosure(dm, setPointID(point), PETSC_TRUE, PETSC_NULL_INTEGER, closure, ierr)) do p = 1, size(closure), 2 PetscCall(PetscSectionGetDoF(section, closure(p), numDof, ierr)) if (numDof > 0) then table(closure(p) + 1, :) = table(closure(p) + 1, :) .or. constraints end if ! numDof end do! p PetscCall(DMPlexRestoreTransitiveClosure(dm, setPointID(point), PETSC_TRUE, PETSC_NULL_INTEGER, closure, ierr)) end do! cell PetscCall(ISRestoreIndices(setPointIS, setPointID, ierr)) end if ! setPointIS PetscCall(ISDestroy(setPointIS, ierr)) end subroutine MEF90SetupConstraintTableSet #undef __FUNCT__ #define __FUNCT__ "MEF90SectionAllocateConstraint" !!! author: Blaise Bourdin (2022, bourdin@mcmaster.ca) !!! !!! MEF90SectionAllocateConstraint: Associates dof to a section !!! subroutine MEF90SectionAllocateConstraint(dm, table, section, ierr) type(tDM), intent(IN) :: dm PetscBool, dimension(:, :), pointer :: table type(tPetscSection), intent(INOUT) :: section PetscErrorCode, intent(INOUT) :: ierr PetscInt :: p, i, pStart, pEnd, numConstraints, numComponents PetscInt, dimension(:), pointer :: constraints PetscInt :: field = 0 PetscCall(DMPlexGetChart(dm, pStart, pEnd, ierr)) do p = 1, pEnd numConstraints = count(table(p, :)) if (numConstraints > 0) then PetscCall(PetscSectionSetConstraintDof(section, p - 1, numConstraints, ierr)) PetscCall(PetscSectionSetFieldConstraintDof(section, p - 1, field, numConstraints, ierr)) end if end do PetscCall(PetscSectionSetup(section, ierr)) numComponents = size(table, 2) do p = 1, pEnd numConstraints = count(table(p, :)) if (numConstraints > 0) then allocate (constraints(numConstraints)) constraints = pack([(i - 1, i=1, numComponents)], table(p, :)) PetscCall(PetscSectionSetConstraintIndices(section, p - 1, constraints, ierr)) PetscCall(PetscSectionSetFieldConstraintIndices(section, p - 1, field, constraints, ierr)) deallocate (constraints) end if end do end subroutine MEF90SectionAllocateConstraint #undef __FUNCT__ #define __FUNCT__ "MEF90VecCopySF" !!! author: Alexis Marboeuf (2022, marboeua@mcmaster.ca) !!! !!! MEF90VecCopySF: rearrange a Vec according to the given SF !!! subroutine MEF90VecCopySF(vin, vout, sf, ierr) type(tVec), intent(IN) :: vin type(tVec), intent(INOUT) :: vout type(tPetscSF), intent(IN) :: sf PetscErrorCode, intent(INOUT) :: ierr PetscScalar, dimension(:), pointer :: arrayin, arrayout PetscCall(VecGetArrayRead(vin, arrayin, ierr)) PetscCall(VecGetArray(vout, arrayout, ierr)) PetscCall(PetscSFBcastBegin(sf, MPIU_SCALAR, arrayin, arrayout, MPI_REPLACE, ierr)) PetscCall(PetscSFBcastEnd(sf, MPIU_SCALAR, arrayin, arrayout, MPI_REPLACE, ierr)) PetscCall(VecRestoreArrayRead(vin, arrayin, ierr)) PetscCall(VecRestoreArray(vout, arrayout, ierr)) end subroutine MEF90VecCopySF #undef __FUNCT__ #define __FUNCT__ "MEF90IOSFCreate" !!! author: Alexis Marboeuf (2022, marboeua@mcmaster.ca) !!! !!! MEF90IOSFCreate: sf mapping between local and IO ordering and distribution !!! subroutine MEF90IOSFCreate(MEF90Ctx, v, liosf, iolsf, ierr) type(tVec), intent(IN) :: v type(tPetscSF), intent(OUT) :: liosf, iolsf type(MEF90Ctx_type), intent(IN) :: MEF90Ctx PetscErrorCode, intent(INOUT) :: ierr type(tPetscSF) :: naturalSF, ioSF, lcgSF, cglSF, tempSF, invTempSF type(tDM) :: dm character(len=PETSC_MAX_PATH_LEN) :: vecname PetscCall(VecGetDM(v, dm, ierr)) PetscCall(CreateLocalToCGlobalSF_Private(MEF90Ctx, dm, lcgSF, ierr)) PetscCall(CreateCGlobalToLocalSF_Private(MEF90Ctx, dm, cglSF, ierr)) PetscCall(PetscObjectGetName(v, vecname, ierr)) if (MEF90Ctx%NumProcs > 1) then PetscCall(CreateNaturalToIOSF_Private(MEF90Ctx, dm, ioSF, ierr)) PetscCall(DMGetNaturalSF(dm, naturalSF, ierr)) PetscCall(PetscSFCompose(lcgSF, naturalSF, tempSF, ierr)) PetscCall(PetscSFCompose(tempSF, ioSF, liosf, ierr)) PetscCall(PetscSFSetUp(liosf, ierr)) PetscCall(PetscSFDestroy(tempSF, ierr)) PetscCall(PetscSFCompose(naturalSF, ioSF, tempSF, ierr)) PetscCall(PetscSFCreateInverseSF(tempSF, invTempSF, ierr)) PetscCall(PetscSFCompose(invTempSF, cglSF, iolsf, ierr)) PetscCall(PetscSFSetUp(iolsf, ierr)) PetscCall(PetscSFDestroy(ioSF, ierr)) PetscCall(PetscSFDestroy(cglSF, ierr)) PetscCall(PetscSFDestroy(lcgSF, ierr)) PetscCall(PetscSFDestroy(tempSF, ierr)) PetscCall(PetscSFDestroy(invTempSF, ierr)) else liosf = lcgSF iolsf = cglSF end if end subroutine MEF90IOSFCreate #undef __FUNCT__ #define __FUNCT__ "MEF90FaceSetIOSFCreate" !!! author: Alexis Marboeuf (2022, marboeua@mcmaster.ca) !!! !!! MEF90FaceSetIOSFCreate: sf mapping between local and IO ordering and distribution for !!! Vec defined on Face Sets !!! subroutine MEF90FaceSetIOSFCreate(MEF90Ctx, v, liosf, iolsf, ierr) type(tVec), intent(IN) :: v type(tPetscSF), intent(OUT) :: liosf, iolsf type(MEF90Ctx_type), intent(IN) :: MEF90Ctx PetscErrorCode, intent(INOUT) :: ierr type(tPetscSF) :: naturalSF, ioSF, lcgSF, cglSF, tempSF, temp2SF, invTempSF, iosideSF, sideioSF type(tDM) :: dm character(len=PETSC_MAX_PATH_LEN) :: vecname PetscCall(VecGetDM(v, dm, ierr)) PetscCall(CreateLocalToCGlobalSF_Private(MEF90Ctx, dm, lcgSF, ierr)) PetscCall(CreateCGlobalToLocalSF_Private(MEF90Ctx, dm, cglSF, ierr)) PetscCall(CreateSideSF_Private(MEF90Ctx, dm, iosideSF, sideioSF, ierr)) PetscCall(PetscObjectGetName(v, vecname, ierr)) if (MEF90Ctx%NumProcs > 1) then PetscCall(CreateNaturalToIOSF_Private(MEF90Ctx, dm, ioSF, ierr)) PetscCall(DMGetNaturalSF(dm, naturalSF, ierr)) PetscCall(PetscSFCompose(lcgSF, naturalSF, tempSF, ierr)) PetscCall(PetscSFCompose(tempSF, iosideSF, temp2SF, ierr)) PetscCall(PetscSFCompose(temp2SF, ioSF, liosf, ierr)) PetscCall(PetscSFDestroy(temp2SF, ierr)) PetscCall(PetscSFSetUp(liosf, ierr)) PetscCall(PetscSFDestroy(tempSF, ierr)) PetscCall(PetscSFCompose(naturalSF, iosideSF, temp2SF, ierr)) PetscCall(PetscSFCompose(temp2SF, ioSF, tempSF, ierr)) PetscCall(PetscSFDestroy(temp2SF, ierr)) PetscCall(PetscSFCreateInverseSF(tempSF, invTempSF, ierr)) PetscCall(PetscSFCompose(invTempSF, cglSF, iolsf, ierr)) PetscCall(PetscSFSetUp(iolsf, ierr)) PetscCall(PetscSFDestroy(ioSF, ierr)) PetscCall(PetscSFDestroy(cglSF, ierr)) PetscCall(PetscSFDestroy(lcgSF, ierr)) PetscCall(PetscSFDestroy(tempSF, ierr)) PetscCall(PetscSFDestroy(invTempSF, ierr)) else PetscCall(PetscSFCompose(lcgSF, iosideSF, liosf, ierr)) PetscCall(PetscSFCompose(sideioSF, cglSF, iolsf, ierr)) PetscCall(PetscSFDestroy(lcgSF, ierr)) PetscCall(PetscSFDestroy(cglSF, ierr)) end if PetscCall(PetscSFDestroy(iosideSF, ierr)) PetscCall(PetscSFDestroy(sideioSF, ierr)) end subroutine MEF90FaceSetIOSFCreate #undef __FUNCT__ #define __FUNCT__ "MEF90ConstraintSFCreate" !!! author: Alexis Marboeuf (2022, marboeua@mcmaster.ca) !!! author: Blaise Bourdin (2022, bourdin@mcmaster.ca) !!! !!! MEF90ConstraintSFCreate: sf mapping between local and constraint Vec ordering and distribution !!! subroutine MEF90ConstraintSFCreate(MEF90Ctx, v, vB, sf, invSF, ierr) type(tVec), intent(IN) :: v, vB type(tPetscSF), intent(OUT) :: sf, invSF type(MEF90Ctx_type), intent(IN) :: MEF90Ctx PetscErrorCode, intent(INOUT) :: ierr type(tDM) :: dm, dmB type(tPetscSection) :: locSection, locBSection type(sPetscSFNode), dimension(:), pointer :: remote PetscInt, dimension(:), pointer :: local, cindices PetscInt :: pStart, pEnd, p, d, nleaves, ldof, loff, cdof, coff, nsize, nroots nleaves = 0 nsize = 0 PetscCall(VecGetDM(v, dm, ierr)) PetscCall(VecGetDM(vB, dmB, ierr)) PetscCall(DMGetLocalSection(dm, locSection, ierr)) PetscCall(PetscSectionGetStorageSize(locSection, nroots, ierr)) PetscCall(DMGetLocalSection(dmB, locBSection, ierr)) PetscCall(PetscSectionGetChart(locBSection, pStart, pEnd, ierr)) do p = pStart, pEnd - 1 PetscCall(PetscSectionGetConstraintDof(locBSection, p, cdof, ierr)) nleaves = nleaves + cdof end do allocate (remote(nleaves)) allocate (local(nleaves)) do p = pStart, pEnd - 1 PetscCall(PetscSectionGetDoF(locSection, p, ldof, ierr)) PetscCall(PetscSectionGetOffset(locSection, p, loff, ierr)) PetscCall(PetscSectionGetConstraintDof(locBSection, p, cdof, ierr)) if (cdof > 0) then PetscCall(PetscSectionGetConstraintIndices(locBSection, p, cindices, ierr)) PetscCall(PetscSectionGetOffset(locBSection, p, coff, ierr)) if (coff >= 0) then do d = 1, cdof local(nsize + 1) = coff + cindices(d) remote(nsize + 1)%rank = MEF90Ctx%rank remote(nsize + 1)%index = loff + cindices(d) nsize = nsize + 1 end do ! d end if ! coff PetscCall(PetscSectionRestoreConstraintIndices(locBSection, p, cindices, ierr)) end if ! cdof end do PetscCall(PetscSFCreate(MEF90Ctx%Comm, sf, ierr)) PetscCall(PetscSFSetFromOptions(sf, ierr)) PetscCall(PetscSFSetGraph(sf, nroots, nleaves, local, PETSC_COPY_VALUES, remote, PETSC_COPY_VALUES, ierr)) deallocate (remote) deallocate (local) PetscCall(PetscSFSetUp(sf, ierr)) PetscCall(PetscSFCreateInverseSF(sf, invSF, ierr)) PetscCall(PetscSFSetUp(invSF, ierr)) end subroutine MEF90ConstraintSFCreate #undef __FUNCT__ #define __FUNCT__ "MEF90VecGlobalToLocalConstraint" !!! author: Blaise Bourdin (2022, bourdin@mcmaster.ca) !!! !!! MEF90VecGlobalToLocalConstraint: do a VecGlobaltoLocal then copy constrained values !!! subroutine MEF90VecGlobalToLocalConstraint(g, c, l, ierr) type(tVec), intent(IN) :: g, c type(tVec), intent(INOUT) :: l PetscErrorCode, intent(INOUT) :: ierr type(tDM) :: dm type(tPetscSection) :: s PetscInt :: numConstraint PetscInt :: p, pStart, pEnd PetscReal, dimension(:), pointer :: vArray PetscCall(VecGetDM(g, dm, ierr)) PetscCall(DMGlobalToLocal(dm, g, INSERT_VALUES, l, ierr)) if (.not. PetscObjectIsNull(c)) then PetscCall(DMGetLocalSection(dm, s, ierr)) PetscCall(PetscSectionGetChart(s, pStart, pEnd, ierr)) do p = pStart, pEnd - 1 PetscCall(PetscSectionGetConstraintDof(s, p, numConstraint, ierr)) if (numConstraint > 0) then PetscCall(VecGetValuesSection(c, s, p, vArray, ierr)) PetscCall(VecSetValuesSection(l, s, p, vArray, INSERT_ALL_VALUES, ierr)) PetscCall(VecRestoreValuesSection(c, s, p, vArray, ierr)) end if ! numConstraint end do ! p end if end subroutine MEF90VecGlobalToLocalConstraint #undef __FUNCT__ #define __FUNCT__ "MEF90VecCreateIO" !!! author: Alexis Marboeuf (2022, marboeua@mcmaster.ca) !!! !!! MEF90VecCreateIO: create IO Vec !!! subroutine MEF90VecCreateIO(v, bs, sf, ierr) type(tPetscSF), intent(IN) :: sf type(tVec), intent(INOUT) :: v PetscInt, intent(IN) :: bs PetscErrorCode, intent(INOUT) :: ierr MPIU_Comm :: comm PetscInt :: nroots, nleaves type(sPetscSFNode), dimension(:), pointer :: iremote PetscInt, dimension(:), pointer :: ilocal PetscCall(PetscSFGetGraph(sf, nroots, nleaves, ilocal, iremote, ierr)) PetscCall(PetscObjectGetComm(sf, comm, ierr)) PetscCall(VecCreateMPI(comm, nleaves, PETSC_DETERMINE, v, ierr)) PetscCall(VecSetBlockSize(v, bs, ierr)) PetscCall(PetscSFrestoreGraph(sf, nroots, nleaves, ilocal, iremote, ierr)) end subroutine MEF90VecCreateIO #undef __FUNCT__ #define __FUNCT__ "MEF90VecSetValuesFromOptions" !!! author: Blaise Bourdin (2022, bourdin@mcmaster.ca) !!! !!! MEF90VecSetValuesFromOptions: Fill values of a Vec using command line options !!! subroutine MEF90VecSetValuesFromOptions(v, scalingFactor, ierr) type(tVec), intent(INOUT) :: v PetscReal, intent(IN) :: scalingFactor PetscErrorCode, intent(INOUT) :: ierr type(tDM) :: dm PetscEnum :: setType PetscInt :: set, point type(tIS) :: setIS, pointIS PetscInt, dimension(:), pointer :: setID, pointID character(len=MEF90MXSTRLEN) :: ValueKey, name PetscBool :: flg PetscInt :: dim, numOpt, bs, numDofClosure, i PetscReal, dimension(:), pointer :: Val, vArray type(tPetscSection) :: section PetscCall(VecGetDM(v, dm, ierr)) PetscCall(PetscObjectGetName(v, name, ierr)) PetscCall(DMGetLocalSection(dm, section, ierr)) PetscCall(DMGetDimension(dm, dim, ierr)) PetscCall(VecGetBlockSize(v, bs, ierr)) allocate (Val(bs)) do setType = 1, size(MEF90SetType) PetscCall(DMGetLabelIdIS(dm, MEF90SetLabelName(setType), setIS, ierr)) ! PetscCall(MEF90ISAllGatherMerge(comm,setIS,ierr)) if (.not. PetscObjectIsNull(setIS)) then PetscCall(ISGetIndices(setIS, setID, ierr)) do set = 1, size(setID) write (ValueKey, '("-",a2,I4.4,"_",a)') MEF90SetPrefix(setType), setID(set), trim(name) numOpt = bs PetscCall(PetscOptionsGetRealArray(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, trim(ValueKey), Val, numOpt, flg, ierr)) if (numOpt > 0) then PetscCall(DMGetStratumIS(dm, MEF90SetLabelName(setType), setID(set), pointIS, ierr)) !! Set the values on the closure of the current point if (.not. PetscObjectIsNull(pointIS)) then PetscCall(ISGetIndices(pointIS, pointID, ierr)) do point = 1, size(pointID) PetscCall(MEF90VecGetClosureSize(v, pointID(point), numDofClosure, ierr)) if (numDofClosure > 0) then PetscCall(DMPlexVecGetClosure(dm, section, v, pointID(point), PETSC_NULL_INTEGER, vArray, ierr)) do i = 1, numDofClosure / bs vArray((i - 1) * bs + 1:i * bs) = scalingFactor * Val end do PetscCall(DMPlexVecSetClosure(dm, section, v, pointID(point), vArray, INSERT_ALL_VALUES, ierr)) PetscCall(DMPlexVecRestoreClosure(dm, section, v, pointID(point), PETSC_NULL_INTEGER, vArray, ierr)) end if ! numDofClosure end do ! point PetscCall(ISRestoreIndices(pointIS, pointID, ierr)) end if ! pointIS PetscCall(ISDestroy(pointIS, ierr)) end if ! numOpt end do ! set PetscCall(ISRestoreIndices(setIS, setID, ierr)) end if ! setIS PetscCall(ISDestroy(setIS, ierr)) end do ! setType deallocate (Val) end subroutine MEF90VecSetValuesFromOptions #undef __FUNCT__ #define __FUNCT__ "MEF90VecSetBCValuesFromOptions" !!! author: Blaise Bourdin (2022, bourdin@mcmaster.ca) !!! !!! MEF90VecSetBCValuesFromOptions: Fill boundary values of a Vec using command line options !!! subroutine MEF90VecSetBCValuesFromOptions(v, scalingFactor, ierr) type(tVec), intent(INOUT) :: v PetscReal, intent(IN) :: scalingFactor PetscErrorCode, intent(INOUT) :: ierr type(tDM) :: dm PetscEnum :: setType PetscInt :: set, point, c type(tIS) :: setIS, pointIS PetscInt, dimension(:), pointer :: setID, pointID character(len=MEF90MXSTRLEN) :: BCOptionKey, BCValueKey, name PetscBool, dimension(:), pointer :: setBC PetscBool :: flg PetscInt :: dim, numBC, bs, numDofClosure PetscReal, dimension(:), pointer :: BCVal, vArray type(tPetscSection) :: section PetscCall(VecGetDM(v, dm, ierr)) PetscCall(PetscObjectGetName(v, name, ierr)) PetscCall(DMGetLocalSection(dm, section, ierr)) PetscCall(DMGetDimension(dm, dim, ierr)) PetscCall(VecGetBlockSize(v, bs, ierr)) allocate (setBC(bs)) allocate (BCVal(bs)) do setType = 1, size(MEF90SetType) PetscCall(DMGetLabelIdIS(dm, MEF90SetLabelName(setType), setIS, ierr)) ! PetscCall(MEF90ISAllGatherMerge(comm,setIS,ierr)) if (.not. PetscObjectIsNull(setIS)) then PetscCall(ISGetIndices(setIS, setID, ierr)) do set = 1, size(setID) setBC = .false. write (BCOptionKey, '("-",a2,I4.4,"_",a,"BC")') MEF90SetPrefix(setType), setID(set), trim(name) numBC = bs PetscCall(PetscOptionsGetBoolArray(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, trim(BCOptionKey), setBC, numBC, flg, ierr)) if (any(setBC)) then !! At least 1 dof has a boundary condition !! Get the unit BC value on the set write (BCValueKey, '("-",a2,I4.4,"_Boundary",a)') MEF90SetPrefix(setType), setID(set), trim(name) numBC = bs PetscCall(PetscOptionsGetRealArray(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, trim(BCValueKey), BCVal, numBC, flg, ierr)) PetscCall(DMGetStratumIS(dm, MEF90SetLabelName(setType), setID(set), pointIS, ierr)) !! Set the boundary values on the closure of the current point if (.not. PetscObjectIsNull(pointIS)) then PetscCall(ISGetIndices(pointIS, pointID, ierr)) do point = 1, size(pointID) PetscCall(MEF90VecGetClosureSize(v, pointID(point), numDofClosure, ierr)) if (numDofClosure > 0) then PetscCall(DMPlexVecGetClosure(dm, section, v, pointID(point), PETSC_NULL_INTEGER, vArray, ierr)) do c = 1, bs if (setBC(c)) then vArray(c::bs) = scalingFactor * BCVal(c) end if ! setBC end do ! c PetscCall(DMPlexVecSetClosure(dm, section, v, pointID(point), vArray, INSERT_ALL_VALUES, ierr)) PetscCall(DMPlexVecRestoreClosure(dm, section, v, pointID(point), PETSC_NULL_INTEGER, vArray, ierr)) end if end do ! point PetscCall(ISRestoreIndices(pointIS, pointID, ierr)) end if ! pointIS PetscCall(ISDestroy(pointIS, ierr)) end if ! setBC end do ! set PetscCall(ISRestoreIndices(setIS, setID, ierr)) end if ! setIS PetscCall(ISDestroy(setIS, ierr)) end do ! setType deallocate (setBC) deallocate (BCVal) end subroutine MEF90VecSetBCValuesFromOptions #undef __FUNCT__ #define __FUNCT__ "MEF90VecSetValuesFromOptionsExpr" !!! author: Blaise Bourdin (2025, bourdin@mcmaster.ca) !!! !!! MEF90VecSetValuesFromOptionsExpr: Fill values of a Vec using expressions passed as petsc options !!! subroutine MEF90VecSetValuesFromOptionsExpr(v, t, ierr) type(tVec), intent(INOUT) :: v PetscReal, intent(IN) :: t PetscErrorCode, intent(INOUT) :: ierr #ifdef MEF90_HAVE_SYMENGINEF90 type(tDM) :: dm PetscEnum :: setType PetscInt :: set, point, p type(tIS) :: setIS, pointIS PetscInt, dimension(:), pointer :: setID, pointID character(len=MEF90MXSTRLEN) :: ValueKey, name, ExprStr, IOBuffer PetscBool :: flg PetscInt :: dim, numOpt, bs, numDofClosure, numDof, i, c, dof PetscInt, dimension(:), pointer :: closure PetscReal, dimension(:), pointer :: vArray type(tPetscSection) :: section type(tPetscSection) :: coordSection type(tVec) :: coordVec PetscReal, dimension(:), pointer :: coordArray PetscReal, dimension(3) :: xyz type(Basic), dimension(:), allocatable :: exprs type(Basic) :: tmpExpr type(symbol), dimension(3) :: vars character :: delim = ';' character(len=MEF90MXSTRLEN), dimension(:), allocatable :: ExprStrComp PetscCall(VecGetDM(v, dm, ierr)) PetscCall(PetscObjectGetName(v, name, ierr)) PetscCall(DMGetLocalSection(dm, section, ierr)) PetscCall(DMGetDimension(dm, dim, ierr)) PetscCall(VecGetBlockSize(v, bs, ierr)) vars = [Symbol("x"), Symbol("y"), Symbol("z")] allocate (exprs(bs)) PetscCall(DMGetCoordinateSection(dm, coordSection, ierr)) PetscCall(DMGetCoordinatesLocal(dm, coordVec, ierr)) do setType = 1, size(MEF90SetType) PetscCall(DMGetLabelIdIS(dm, MEF90SetLabelName(setType), setIS, ierr)) ! PetscCall(MEF90ISAllGatherMerge(comm,setIS,ierr)) if (.not. PetscObjectIsNull(setIS)) then PetscCall(ISGetIndices(setIS, setID, ierr)) do set = 1, size(setID) write (ValueKey, '("-",a2,I4.4,"_",a,"Expr")') MEF90SetPrefix(setType), setID(set), trim(name) PetscCall(PetscOptionsGetString(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, trim(ValueKey), ExprStr, flg, ierr)) if (flg) then call MEF90StrTokenize(ExprStr, delim, ExprStrComp) numOpt = size(ExprStrComp) if (numOpt < bs) then write (IOBuffer, "(A,' was expecting ',I2,' expressions but got ',I2,' for key ',A)") __FUNCT__, bs, size(ExprStrComp), trim(ValueKey) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, IOBuffer) end if !! create parsers do c = 1, bs exprs(c) = parse(ExprStrComp(c)) exprs(c) = exprs(c)%subs(Symbol("t"), RealDouble(t)) end do PetscCall(DMGetStratumIS(dm, MEF90SetLabelName(setType), setID(set), pointIS, ierr)) !! Set the values on the closure of the current point if (.not. PetscObjectIsNull(pointIS)) then PetscCall(ISGetIndices(pointIS, pointID, ierr)) do point = 1, size(pointID) PetscCall(DMPlexGetTransitiveClosure(dm, pointID(point), PETSC_TRUE, numDofClosure, closure, ierr)) if (numDofClosure > 0) then PetscCall(DMPlexVecGetClosure(dm, section, v, pointID(point), PETSC_NULL_INTEGER, vArray, ierr)) dof = 0 do p = 1, size(closure), 2 PetscCall(PetscSectionGetDof(section, closure(p), numDof, ierr)) if (numDof > 0) then dof = dof + 1 !! Get the coordinates of the dof associated with the point !! trick: the coordinate of a point is the average of the coordinates of the points in its closure PetscCall(DMPlexVecGetClosure(dm, coordSection, coordVec, closure(p), PETSC_NULL_INTEGER, coordArray, ierr)) do i = 1, dim xyz(i) = sum(coordArray(i:size(coordArray):dim)) * dim / size(coordArray) end do PetscCall(DMPlexVecRestoreClosure(dm, coordSection, coordVec, closure(p), PETSC_NULL_INTEGER, coordArray, ierr)) do c = 1, bs tmpExpr = Exprs(c) do i = 1, dim tmpExpr = tmpExpr%subs(vars(i), RealDouble(xyz(i))) end do ! i tmpExpr = tmpExpr%evalf() vArray((dof - 1) * bs + c) = tmpExpr%dbl() end do ! c end if !bnumDof end do !p PetscCall(DMPlexVecSetClosure(dm, section, v, pointID(point), vArray, INSERT_ALL_VALUES, ierr)) PetscCall(DMPlexVecRestoreClosure(dm, section, v, pointID(point), PETSC_NULL_INTEGER, vArray, ierr)) end if ! numDofClosure PetscCall(DMPlexRestoreTransitiveClosure(dm, pointID(point), PETSC_TRUE, numDofClosure, closure, ierr)) end do ! point PetscCall(ISRestoreIndices(pointIS, pointID, ierr)) end if ! pointIS PetscCall(ISDestroy(pointIS, ierr)) deallocate (ExprStrComp) end if ! flg end do ! set PetscCall(ISRestoreIndices(setIS, setID, ierr)) end if ! setIS PetscCall(ISDestroy(setIS, ierr)) end do ! setType deallocate (exprs) #else write (*, *) "ERROR: ", __FUNCT__, " requires symengine-f90 support" stop #endif end subroutine MEF90VecSetValuesFromOptionsExpr #undef __FUNCT__ #define __FUNCT__ "MEF90VecSetBCValuesFromOptionsExpr" !!! author: Blaise Bourdin (2025, bourdin@mcmaster.ca) !!! !!! MEF90VecSetBCValuesFromOptionsExpr: Fill boundary values of a Vec using expressions passed as petsc options !!! subroutine MEF90VecSetBCValuesFromOptionsExpr(v, t, ierr) type(tVec), intent(INOUT) :: v PetscReal, intent(IN) :: t PetscErrorCode, intent(INOUT) :: ierr #ifdef MEF90_HAVE_SYMENGINEF90 type(tDM) :: dm PetscEnum :: setType PetscInt :: set, point, c, p, i type(tIS) :: setIS, pointIS PetscInt, dimension(:), pointer :: setID, pointID character(len=MEF90MXSTRLEN) :: BCOptionKey, BCValueKey, name, ExprStr, IOBuffer PetscBool, dimension(:), pointer :: setBC PetscBool :: flg PetscInt :: dim, numOpt, numBC, bs, numDofClosure, numDof, dof PetscInt, dimension(:), pointer :: closure PetscReal, dimension(:), pointer :: vArray type(tPetscSection) :: section type(tPetscSection) :: coordSection type(tVec) :: coordVec PetscReal, dimension(:), pointer :: coordArray PetscReal, dimension(3) :: xyz type(Basic), dimension(:), allocatable :: exprs type(Basic) :: tmpExpr type(symbol), dimension(3) :: vars character :: delim = ';' character(len=MEF90MXSTRLEN), dimension(:), allocatable :: ExprStrComp PetscCall(VecGetDM(v, dm, ierr)) PetscCall(PetscObjectGetName(v, name, ierr)) PetscCall(DMGetLocalSection(dm, section, ierr)) PetscCall(DMGetDimension(dm, dim, ierr)) PetscCall(VecGetBlockSize(v, bs, ierr)) allocate (setBC(bs)) vars = [Symbol("x"), Symbol("y"), Symbol("z")] allocate (exprs(bs)) PetscCall(DMGetCoordinateSection(dm, coordSection, ierr)) PetscCall(DMGetCoordinatesLocal(dm, coordVec, ierr)) do setType = 1, size(MEF90SetType) PetscCall(DMGetLabelIdIS(dm, MEF90SetLabelName(setType), setIS, ierr)) ! PetscCall(MEF90ISAllGatherMerge(comm,setIS,ierr)) if (.not. PetscObjectIsNull(setIS)) then PetscCall(ISGetIndices(setIS, setID, ierr)) do set = 1, size(setID) setBC = .false. write (BCOptionKey, '("-",a2,I4.4,"_",a,"BC")') MEF90SetPrefix(setType), setID(set), trim(name) numBC = bs PetscCall(PetscOptionsGetBoolArray(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, trim(BCOptionKey), setBC, numBC, flg, ierr)) if (any(setBC)) then !! At least 1 dof has a boundary condition !! Get the unit BC value on the set write (BCValueKey, '("-",a2,I4.4,"_Boundary",a,"Expr")') MEF90SetPrefix(setType), setID(set), trim(name) PetscCall(PetscOptionsGetString(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, trim(BCValueKey), ExprStr, flg, ierr)) numOpt = MEF90StrCount(ExprStr, delim) call MEF90StrTokenize(ExprStr, delim, ExprStrComp) if (size(ExprStrComp) < bs) then write (IOBuffer, "(A,' was expecting ',I2,' expressions but got ',I2,' for key ',A)") __FUNCT__, bs, size(ExprStrComp), trim(BCValueKey) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, IOBuffer) end if !! create parsers do c = 1, bs exprs(c) = parse(ExprStrComp(c)) exprs(c) = exprs(c)%subs(Symbol("t"), RealDouble(t)) end do PetscCall(DMGetStratumIS(dm, MEF90SetLabelName(setType), setID(set), pointIS, ierr)) !! Set the boundary values on the closure of the current point if (.not. PetscObjectIsNull(pointIS)) then PetscCall(ISGetIndices(pointIS, pointID, ierr)) do point = 1, size(pointID) PetscCall(DMPlexGetTransitiveClosure(dm, pointID(point), PETSC_TRUE, numDofClosure, closure, ierr)) if (numDofClosure > 0) then PetscCall(DMPlexVecGetClosure(dm, section, v, pointID(point), PETSC_NULL_INTEGER, vArray, ierr)) dof = 0 do p = 1, size(closure), 2 PetscCall(PetscSectionGetDof(section, closure(p), numDof, ierr)) if (numDof > 0) then dof = dof + 1 !! Get the coordinates of the dof associated with the point !! trick: the coordinate of a point is the average of the coordinates of the points in its closure PetscCall(DMPlexVecGetClosure(dm, coordSection, coordVec, closure(p), PETSC_NULL_INTEGER, coordArray, ierr)) do i = 1, dim xyz(i) = sum(coordArray(i:size(coordArray):dim)) * dim / size(coordArray) end do PetscCall(DMPlexVecRestoreClosure(dm, coordSection, coordVec, closure(p), PETSC_NULL_INTEGER, coordArray, ierr)) do c = 1, bs if (setBC(c)) then tmpExpr = Exprs(c) do i = 1, dim tmpExpr = tmpExpr%subs(vars(i), RealDouble(xyz(i))) end do ! i tmpExpr = tmpExpr%evalf() vArray((dof - 1) * bs + c) = tmpExpr%dbl() end if end do ! c end if !bnumDof end do !p PetscCall(DMPlexVecSetClosure(dm, section, v, pointID(point), vArray, INSERT_ALL_VALUES, ierr)) PetscCall(DMPlexVecRestoreClosure(dm, section, v, pointID(point), PETSC_NULL_INTEGER, vArray, ierr)) end if ! numDofClosure PetscCall(DMPlexRestoreTransitiveClosure(dm, pointID(point), PETSC_TRUE, numDofClosure, closure, ierr)) end do ! point PetscCall(ISRestoreIndices(pointIS, pointID, ierr)) end if ! pointIS PetscCall(ISDestroy(pointIS, ierr)) deallocate (ExprStrComp) end if ! setBC end do ! set PetscCall(ISRestoreIndices(setIS, setID, ierr)) end if ! setIS PetscCall(ISDestroy(setIS, ierr)) end do ! setType deallocate (setBC) deallocate (exprs) #else write (*, *) "ERROR: ", __FUNCT__, " requires symengine-f90 support" stop #endif end subroutine MEF90VecSetBCValuesFromOptionsExpr !! Private functions below #undef __FUNCT__ #define __FUNCT__ "CreateNaturalToIOSF_Private" !!! author: Alexis Marboeuf (2022, marboeua@mcmaster.ca) !!! !!! CreateNaturalToIOSF_Private: !!! subroutine CreateNaturalToIOSF_Private(MEF90Ctx, dm, sf, ierr) type(tDM), intent(IN) :: dm type(MEF90Ctx_type), intent(IN) :: MEF90Ctx type(tPetscSF), intent(OUT) :: sf PetscErrorCode, intent(INOUT) :: ierr type(tVec) :: vnat, vio type(tPetscLayout) :: ioMap, natMap type(sPetscSFNode), dimension(:), pointer :: remote PetscInt, dimension(:), pointer :: ioRange PetscInt :: nroots, nleaves, globalIndex, i, globalSize, bs PetscCall(DMPlexCreateNaturalVector(dm, vnat, ierr)) PetscCall(VecGetSize(vnat, globalSize, ierr)) PetscCall(VecGetBlockSize(vnat, bs, ierr)) PetscCall(VecCreate(MEF90Ctx%Comm, vio, ierr)) PetscCall(VecSetBlockSize(vio, bs, ierr)) PetscCall(VecSetSizes(vio, PETSC_DETERMINE, globalSize, ierr)) PetscCall(VecSetFromOptions(vio, ierr)) PetscCall(VecGetLayout(vio, ioMap, ierr)) PetscCall(VecGetLayout(vnat, natMap, ierr)) PetscCall(PetscLayoutGetLocalSize(natMap, nroots, ierr)) PetscCall(PetscLayoutGetLocalSize(ioMap, nleaves, ierr)) PetscCall(PetscLayoutGetRanges(ioMap, ioRange, ierr)) allocate (remote(nleaves)) do i = 0, nleaves - 1 globalIndex = ioRange(MEF90Ctx%rank + 1) + i remote(i + 1)%rank = 0 remote(i + 1)%index = globalIndex end do PetscCall(PetscSFCreate(MEF90Ctx%Comm, sf, ierr)) ! PetscCall(PetscObjectSetName(sf,"Natural-To-IO SF",ierr)) PetscCall(PetscSFSetFromOptions(sf, ierr)) PetscCall(PetscSFSetGraph(sf, nroots, nleaves, PETSC_NULL_INTEGER_ARRAY, PETSC_COPY_VALUES, remote, PETSC_COPY_VALUES, ierr)) PetscCall(PetscSFSetUp(sf, ierr)) ! PetscCall(PetscSFViewFromOptions(sf,PETSC_NULL_OBJECT,"-naturaltoio_sf_view",ierr)) PetscCall(VecDestroy(vio, ierr)) PetscCall(VecDestroy(vnat, ierr)) deallocate (remote) end subroutine CreateNaturalToIOSF_Private #undef __FUNCT__ #define __FUNCT__ "CreateLocalToCGlobalSF_Private" !!! author: Alexis Marboeuf (2022, marboeua@mcmaster.ca) !!! !!! CreateLocalToCGlobalSF_Private: !!! subroutine CreateLocalToCGlobalSF_Private(MEF90Ctx, dm, sf, ierr) type(tDM), intent(IN) :: dm type(MEF90Ctx_type), intent(IN) :: MEF90Ctx type(tPetscSF), intent(OUT) :: sf PetscErrorCode, intent(INOUT) :: ierr type(tPetscSection) :: locSection, gSection type(tPetscSF) :: overlapSF, idSF type(sPetscSFNode), dimension(:), pointer :: remote PetscInt, dimension(:), pointer :: remoteOffsets PetscInt :: pStart, pEnd, p, n PetscCall(DMGetLocalSection(dm, locSection, ierr)) PetscCall(DMGetPointSF(dm, overlapSF, ierr)) PetscCall(PetscSectionCreateGlobalSection(locSection, overlapSF, PETSC_TRUE, PETSC_TRUE, PETSC_TRUE, gSection, ierr)) PetscCall(PetscSectionGetChart(locSection, pStart, pEnd, ierr)) n = pEnd - pStart allocate (remote(n)) do p = 1, n remote(p)%rank = MEF90Ctx%rank remote(p)%index = p - 1 end do PetscCall(PetscSFCreate(MEF90Ctx%Comm, idSF, ierr)) PetscCall(PetscSFSetFromOptions(idSF, ierr)) PetscCall(PetscSFSetGraph(idSF, n, n, PETSC_NULL_INTEGER_ARRAY, PETSC_COPY_VALUES, remote, PETSC_COPY_VALUES, ierr)) PetscCall(PetscSFSetUp(idSF, ierr)) PetscCall(PetscSFCreateRemoteOffsets(idSF, locSection, gSection, remoteOffsets, ierr)) PetscCall(PetscSFCreateSectionSF(idSF, locSection, remoteOffsets, gSection, sf, ierr)) if (associated(remoteOffsets)) then PetscCall(PetscSFDestroyRemoteOffsets(remoteOffsets, ierr)) end if PetscCall(PetscSFSetUp(sf, ierr)) ! PetscCall(PetscObjectSetName(sf,"Local-To-CGlobal SF",ierr)) PetscCall(PetscSFSetFromOptions(sf, ierr)) ! PetscCall(PetscSFViewFromOptions(sf,PETSC_NULL_OBJECT,"-localtocglobal_sf_view",ierr)) PetscCall(PetscSectionDestroy(gSection, ierr)) PetscCall(PetscSFDestroy(idSF, ierr)) deallocate (remote) end subroutine CreateLocalToCGlobalSF_Private #undef __FUNCT__ #define __FUNCT__ "CreateCGlobalToLocalSF_Private" !!! author: Alexis Marboeuf (2022, marboeua@mcmaster.ca) !!! !!! CreateCGlobalToLocalSF_Private: !!! subroutine CreateCGlobalToLocalSF_Private(MEF90Ctx, dm, sf, ierr) type(tDM), intent(IN) :: dm type(MEF90Ctx_type), intent(IN) :: MEF90Ctx type(tPetscSF), intent(OUT) :: sf PetscErrorCode, intent(INOUT) :: ierr type(tPetscSection) :: locSection, gSection type(tPetscSF) :: overlapSF, idSF, tempSF, ttempSF type(sPetscSFNode), dimension(:), pointer :: remote, tempRemote, lgRemote, glRemote PetscInt, dimension(:), pointer :: tempLocal, lgLocal, glLocal, remoteOffsets PetscInt :: pStart, pEnd, p, n, lgNRoots, lgNLeaves, tempNRoots, tempNLeaves, glNRoots, glNLeaves PetscCall(DMGetLocalSection(dm, locSection, ierr)) PetscCall(DMGetPointSF(dm, overlapSF, ierr)) PetscCall(PetscSectionCreateGlobalSection(locSection, overlapSF, PETSC_TRUE, PETSC_TRUE, PETSC_TRUE, gSection, ierr)) PetscCall(PetscSectionGetChart(locSection, pStart, pEnd, ierr)) n = pEnd - pStart allocate (remote(n)) do p = 1, n remote(p)%rank = MEF90Ctx%rank remote(p)%index = p - 1 end do PetscCall(PetscSFCreate(MEF90Ctx%Comm, idSF, ierr)) PetscCall(PetscSFSetFromOptions(idSF, ierr)) PetscCall(PetscSFSetGraph(idSF, n, n, PETSC_NULL_INTEGER_ARRAY, PETSC_COPY_VALUES, remote, PETSC_COPY_VALUES, ierr)) PetscCall(PetscSFSetUp(idSF, ierr)) PetscCall(PetscSFCreateRemoteOffsets(idSF, locSection, gSection, remoteOffsets, ierr)) PetscCall(PetscSFCreateSectionSF(idSF, locSection, remoteOffsets, gSection, tempSF, ierr)) PetscCall(PetscSFCreateInverseSF(tempSF, sf, ierr)) PetscCall(PetscSFDestroyRemoteOffsets(remoteOffsets, ierr)) PetscCall(PetscSFDestroy(tempSF, ierr)) if (MEF90Ctx%NumProcs > 1) then PetscCall(PetscSFGetGraph(sf, lgNRoots, lgNLeaves, lgLocal, lgRemote, ierr)) PetscCall(PetscSFCreateRemoteOffsets(overlapSF, locSection, gSection, remoteOffsets, ierr)) PetscCall(PetscSFCreateSectionSF(overlapSF, locSection, remoteOffsets, gSection, ttempSF, ierr)) PetscCall(PetscSFCreateInverseSF(ttempSF, tempSF, ierr)) PetscCall(PetscSFDestroyRemoteOffsets(remoteOffsets, ierr)) PetscCall(PetscSFDestroy(ttempSF, ierr)) PetscCall(PetscSFGetGraph(tempSF, tempNRoots, tempNLeaves, tempLocal, tempRemote, ierr)) glNRoots = lgNRoots glNLeaves = lgNLeaves + tempNLeaves allocate (glLocal(glNLeaves)) allocate (glRemote(glNLeaves)) if (loc(lgLocal) /= loc(PETSC_NULL_INTEGER)) then do p = 1, lgNLeaves glLocal(p) = lgLocal(p) glRemote(p)%rank = lgRemote(p)%rank glRemote(p)%index = lgRemote(p)%index end do else do p = 1, lgNLeaves glLocal(p) = p - 1 glRemote(p)%rank = lgRemote(p)%rank glRemote(p)%index = lgRemote(p)%index end do end if if (loc(tempLocal) /= loc(PETSC_NULL_INTEGER)) then do p = 1, tempNLeaves glLocal(p + lgNLeaves) = tempLocal(p) glRemote(p + lgNLeaves)%rank = tempRemote(p)%rank glRemote(p + lgNLeaves)%index = tempRemote(p)%index end do else do p = 1, tempNLeaves glLocal(p + lgNLeaves) = p + lgNLeaves - 1 glRemote(p + lgNLeaves)%rank = tempRemote(p)%rank glRemote(p + lgNLeaves)%index = tempRemote(p)%index end do end if PetscCall(PetscSFDestroy(sf, ierr)) PetscCall(PetscSFCreate(MEF90Ctx%Comm, sf, ierr)) PetscCall(PetscSFSetFromOptions(sf, ierr)) PetscCall(PetscSFSetGraph(sf, glNRoots, glNLeaves, glLocal, PETSC_COPY_VALUES, glRemote, PETSC_COPY_VALUES, ierr)) PetscCall(PetscSFSetUp(sf, ierr)) PetscCall(PetscSFRestoreGraph(tempSF, tempNRoots, tempNLeaves, tempLocal, tempRemote, ierr)) deallocate (glLocal) deallocate (glRemote) ! PetscCall(PetscSFRestoreGraph(tempSF, tempNRoots, tempNLeaves, tempLocal, tempRemote, ierr)) PetscCall(PetscSFRestoreGraph(sf, lgNRoots, lgNLeaves, lgLocal, lgRemote, ierr)) PetscCall(PetscSFDestroy(tempSF, ierr)) end if PetscCall(PetscSectionDestroy(gSection, ierr)) PetscCall(PetscSFDestroy(idSF, ierr)) ! PetscCall(PetscObjectSetName(sf,"CGlobal-To-Local SF",ierr)) ! PetscCall(PetscSFViewFromOptions(sf,PETSC_NULL_OBJECT,"-cglobaltolocal_sf_view",ierr)) deallocate (remote) end subroutine CreateCGlobalToLocalSF_Private #undef __FUNCT__ #define __FUNCT__ "CreateSideSF_Private" !!! author: Alexis Marboeuf (2022, marboeua@mcmaster.ca) !!! !!! CreateSideSF_Private: !!! subroutine CreateSideSF_Private(MEF90Ctx, dm, sf, invSF, ierr) type(tDM), intent(IN) :: dm type(MEF90Ctx_type), intent(IN) :: MEF90Ctx type(tPetscSF), intent(OUT) :: sf, invSF PetscErrorCode, intent(INOUT) :: ierr type(tIS) :: ssIS, gssIS, faceIS, facesIS PetscInt, dimension(:), pointer :: ssID, faceID type(tPetscSF) :: migrationSF, tempSF type(tVec) :: localVec PetscInt :: set, face, nroots, nleaves, i, totalleaves, numComponent, uNumComponent, j, numSS, key, numFaces type(sPetscSFNode), dimension(:), pointer :: iremote type(tIS), dimension(:), pointer :: locfacesIS PetscInt, dimension(:), pointer :: ilocal, permIndices, emptyInd, facesID, procSSID nleaves = 0_ki totalleaves = 0_ki PetscCall(DMGetLabelIdIS(dm, "Face Sets", gssIS, ierr)) PetscCall(MEF90ISAllGatherMerge(MEF90Ctx%comm, gssIS, ierr)) PetscCall(ISGetSize(gssIS, numSS, ierr)) PetscCall(DMGetLabelIdIS(dm, "Face Sets", ssIS, ierr)) PetscCall(DMPlexGetMigrationSF(dm, migrationSF, ierr)) PetscCall(ISGetIndices(ssIS, ssID, ierr)) allocate (locfacesIS(numSS)) allocate (procSSID(numSS)) do set = 1, numSS procSSID(set) = -1 end do do set = 1, size(ssID) PetscCall(ISLocate(gssIS, ssID(set), key, ierr)) if (key >= 0) then procSSID(key + 1) = ssID(set) end if end do do set = 1, numSS PetscCall(DMGetStratumIS(dm, "Face Sets", procSSID(set), faceIS, ierr)) if (.not. PetscObjectIsNull(faceIS)) then PetscCall(ISGetIndices(faceIS, faceID, ierr)) else allocate (faceID(0)) end if numFaces = size(faceID) allocate (ilocal(numFaces)) allocate (iremote(numFaces)) do face = 1, numFaces iremote(face)%rank = MEF90Ctx%Rank iremote(face)%index = face - 1 ilocal(face) = faceID(face) end do PetscCall(PetscSFCreate(MEF90Ctx%Comm, tempSF, ierr)) PetscCall(PetscSFSetFromOptions(tempSF, ierr)) PetscCall(PetscSFSetGraph(tempSF, numFaces, numFaces, ilocal, PETSC_COPY_VALUES, iremote, PETSC_COPY_VALUES, ierr)) PetscCall(PetscSFSetUp(tempSF, ierr)) if (MEF90Ctx%NumProcs > 1) then PetscCall(PetscSFComposeInverse(migrationSF, tempSF, sf, ierr)) PetscCall(PetscSFDestroy(tempSF, ierr)) else PetscCall(PetscSFCreateInverseSF(tempSF, sf, ierr)) PetscCall(PetscSFDestroy(tempSF, ierr)) end if PetscCall(PetscSFSetUp(sf, ierr)) deallocate (ilocal) deallocate (iremote) if (.not. PetscObjectIsNull(faceIS)) then PetscCall(ISRestoreIndices(faceIS, faceID, ierr)) else deallocate (faceID) end if PetscCall(PetscSFGetGraph(sf, nroots, nleaves, emptyInd, iremote, ierr)) allocate (ilocal(nleaves)) do i = 1, nleaves ilocal(i) = iremote(i)%index end do PetscCall(ISCreateGeneral(MEF90Ctx%comm, nleaves, ilocal, PETSC_COPY_VALUES, locfacesIS(set), ierr)) PetscCall(ISDestroy(faceIS, ierr)) PetscCall(PetscSFRestoreGraph(sf, nroots, nleaves, emptyInd, iremote, ierr)) PetscCall(PetscSFDestroy(sf, ierr)) deallocate (ilocal) end do do set = 1, numSS PetscCall(MEF90ISAllGatherMerge(MEF90Ctx%comm, locfacesIS(set), ierr)) end do PetscCall(ISConcatenate(MEF90Ctx%comm, numSS, locfacesIS, facesIS, ierr)) PetscCall(ISGetIndices(facesIS, faceID, ierr)) if (MEF90Ctx%rank == 0) then totalleaves = size(faceID) end if allocate (facesID(totalleaves)) allocate (permIndices(totalleaves)) do i = 1, totalleaves permIndices(i) = i - 1 facesID(i) = faceID(i) end do PetscCall(ISRestoreIndices(facesIS, faceID, ierr)) PetscCall(ISRestoreIndices(ssIS, ssID, ierr)) PetscCall(PetscSortIntWithPermutation(totalleaves, facesID, permIndices, ierr)) PetscCall(DMGetLocalVector(dm, localVec, ierr)) PetscCall(VecGetBlockSize(localVec, numComponent, ierr)) PetscCall(MPI_Allreduce(numComponent, uNumComponent, 1, MPIU_INTEGER, MPI_MAX, MEF90Ctx%comm, ierr)) if ((numComponent == 1) .and. (uNumComponent > 1)) then numComponent = uNumComponent end if PetscCall(DMRestoreLocalVector(dm, localVec, ierr)) allocate (iremote(numComponent * totalleaves)) allocate (ilocal(numComponent * totalleaves)) do i = 1, totalleaves do j = 1, numComponent ilocal((i - 1) * numComponent + j) = (i - 1) * numComponent + j - 1 iremote((i - 1) * numComponent + j)%rank = MEF90Ctx%rank iremote((i - 1) * numComponent + j)%index = numComponent * permIndices(i) + j - 1 end do end do PetscCall(PetscSFCreate(MEF90Ctx%Comm, invSF, ierr)) PetscCall(PetscSFSetFromOptions(invSF, ierr)) PetscCall(PetscSFSetGraph(invSF, numComponent * totalleaves, numComponent * totalleaves, ilocal, PETSC_COPY_VALUES, iremote, PETSC_COPY_VALUES, ierr)) PetscCall(PetscSFSetUp(invSF, ierr)) PetscCall(PetscSFCreateInverseSF(invSF, sf, ierr)) PetscCall(PetscSFSetUp(sf, ierr)) do set = 1, numSS PetscCall(ISDestroy(locfacesIS(set), ierr)) end do deallocate (locfacesIS) deallocate (permIndices) deallocate (procSSID) deallocate (facesID) deallocate (ilocal) deallocate (iremote) PetscCall(ISDestroy(gssIS, ierr)) PetscCall(ISDestroy(ssIS, ierr)) PetscCall(ISDestroy(facesIS, ierr)) end subroutine CreateSideSF_Private end module m_MEF90_DMPlex