ExtractSubMesh Subroutine

  subroutine ExtractSubMesh(globalMesh,iset,setOfElem,setGlobalId,setOffsets,comm,subMesh)
    !! Build the sub-mesh of element set iset from this rank's slice of the global mesh.
    !!
    !! The sub-mesh holds this rank's elements of set iset, in global element order, with
    !! the global mesh's geometry, materials and boundary-condition table. Its
    !! decomposition is built from setOffsets(:,iset), so element ownership is unchanged;
    !! a rank may own zero elements. Faces are renumbered as described in the module
    !! documentation, and mortars inside the set are copied with renumbered element ids.
    !!
    !! The decomposition is initialised on comm, which must be a duplicate of
    !! globalMesh%decomp%mpiComm made for this sub-mesh alone (MPI_Comm_dup). The
    !! decomposition does not own it: the caller frees comm with MPI_Comm_free after
    !! subMesh%Free. A separate communicator per sub-mesh keeps the halo messages of
    !! different sets from matching each other.
    !!
    !! Collective over comm: every rank calls it for every set.
    !!
    !!  Input
    !!    - globalMesh : the decomposed global mesh
    !!    - iset : the set to extract, 1 <= iset <= size(setOffsets,2)
    !!    - setOfElem, setGlobalId, setOffsets : from GatherSetOfElem and
    !!      BuildSetNumbering, with the partition validated by ValidateSetPartition
    !!    - comm : this sub-mesh's duplicate communicator
    !!  Output
    !!    - subMesh : the sub-mesh, device copy up to date
    implicit none
    type(Mesh3D),intent(in) :: globalMesh
    integer,intent(in) :: iset
    integer,intent(in) :: setOfElem(1:)
    integer,intent(in) :: setGlobalId(1:)
    integer,intent(in) :: setOffsets(1:,1:)
    integer,intent(in) :: comm
    type(Mesh3D),intent(out) :: subMesh
    ! Local
    integer :: nGeo,nBCs,nLocal,rank,goff,el,eg,e2,s,n,k,m,q,a,nm,na
    integer,allocatable :: subLocal(:)
    integer,allocatable :: mortarNew(:),anisoNew(:)

    nGeo = globalMesh%nGeo
    nBCs = globalMesh%nBCs

    ! The decomposition pointers carry no null default, and GenerateDecompositionFromOffsets
    ! tests their association; subMesh may live in uninitialised (heap) memory.
    nullify(subMesh%decomp%elemToRank,subMesh%decomp%offSetElem)
    call subMesh%decomp%Init(comm=comm)
    ! Start the request scratch at one entry: the global mesh's size would be copied into
    ! every set on every rank. MultiSetModel3D%SetupInterfaces reserves each set's actual
    ! requirement once its variable count is known, and every exchange reserves before
    ! posting in any case.
    call subMesh%decomp%GenerateDecompositionFromOffsets(setOffsets(:,iset),1)

    rank = subMesh%decomp%rankId
    nLocal = setOffsets(rank+2,iset)-setOffsets(rank+1,iset)

    call subMesh%Init(nGeo,nLocal,6*nLocal,8*nLocal,nBCs)
    subMesh%nGlobalElem = setOffsets(subMesh%decomp%nRanks+1,iset)
    subMesh%nUniqueSides = globalMesh%nUniqueSides ! global side count: the CPU halo tag stride
    subMesh%nUniqueNodes = globalMesh%nUniqueNodes
    subMesh%nCornerNodes = globalMesh%nCornerNodes
    subMesh%quadrature = globalMesh%quadrature
    subMesh%nLayers = 0 ! a set need not contain whole columns of an extruded mesh

    ! Boundary-condition and material tables (replicated on every rank).
    if(nBCs > 0) then
      subMesh%BCType(1:4,1:nBCs) = globalMesh%BCType(1:4,1:nBCs)
      do k = 1,nBCs
        subMesh%BCNames(k) = globalMesh%BCNames(k)
      enddo
    endif
    subMesh%nMaterials = globalMesh%nMaterials
    if(allocated(subMesh%materialNames)) deallocate(subMesh%materialNames)
    allocate(subMesh%materialNames(1:globalMesh%nMaterials))
    subMesh%materialNames(1:globalMesh%nMaterials) = globalMesh%materialNames(1:globalMesh%nMaterials)

    ! Per-element data. subLocal(el) is the sub-mesh local id of global-mesh local
    ! element el, or 0 when el belongs to another set.
    goff = globalMesh%decomp%offsetElem(rank+1)
    allocate(subLocal(1:max(globalMesh%nElem,1)))
    subLocal = 0
    n = 0
    do el = 1,globalMesh%nElem
      eg = goff+el
      if(setOfElem(eg) /= iset) cycle
      n = n+1
      subLocal(el) = n
      subMesh%nodeCoords(1:3,1:nGeo+1,1:nGeo+1,1:nGeo+1,n) = &
        globalMesh%nodeCoords(1:3,1:nGeo+1,1:nGeo+1,1:nGeo+1,el)
      subMesh%globalNodeIDs(1:nGeo+1,1:nGeo+1,1:nGeo+1,n) = &
        globalMesh%globalNodeIDs(1:nGeo+1,1:nGeo+1,1:nGeo+1,el)
      subMesh%elemInfo(1:6,n) = globalMesh%elemInfo(1:6,el)
      subMesh%elemMaterial(n) = globalMesh%elemMaterial(el)

      do s = 1,6
        e2 = globalMesh%sideInfo(3,s,el)
        subMesh%sideInfo(1:5,s,n) = globalMesh%sideInfo(1:5,s,el)
        if(e2 > 0) then
          if(setOfElem(e2) == iset) then
            ! Same-set neighbour: renumber; side, flip and global side id are unchanged.
            subMesh%sideInfo(3,s,n) = setGlobalId(e2)
          else
            ! Neighbour in another set: an interface face of this sub-mesh.
            subMesh%sideInfo(1,s,n) = 0
            subMesh%sideInfo(3,s,n) = 0
            subMesh%sideInfo(4,s,n) = 0
            subMesh%sideInfo(5,s,n) = SELF_BC_SET_INTERFACE
          endif
        endif
      enddo
    enddo

    if(n /= nLocal) then
      error stop "ExtractSubMesh : local element count disagrees with setOffsets"
    endif

    ! Mortars wholly inside this set (ValidateSetPartition rejects any other), renumbered
    ! in their global order. The tables are replicated, so every rank copies all of them.
    nm = 0
    allocate(mortarNew(1:max(globalMesh%nMortars,1)))
    mortarNew = 0
    do m = 1,globalMesh%nMortars
      if(setOfElem(globalMesh%mortarInfo(1,m)) == iset) then
        nm = nm+1
        mortarNew(m) = nm
      endif
    enddo
    na = 0
    allocate(anisoNew(1:max(globalMesh%nAnisoMortars,1)))
    anisoNew = 0
    do a = 1,globalMesh%nAnisoMortars
      if(setOfElem(globalMesh%anisoMortarInfo(1,a)) == iset) then
        na = na+1
        anisoNew(a) = na
      endif
    enddo

    subMesh%nMortars = nm
    if(nm > 0) then
      allocate(subMesh%mortarInfo(1:14,1:nm))
      do m = 1,globalMesh%nMortars
        if(mortarNew(m) == 0) cycle
        subMesh%mortarInfo(1:14,mortarNew(m)) = globalMesh%mortarInfo(1:14,m)
        subMesh%mortarInfo(1,mortarNew(m)) = setGlobalId(globalMesh%mortarInfo(1,m))
        do q = 1,4
          subMesh%mortarInfo(1+2*q,mortarNew(m)) = setGlobalId(globalMesh%mortarInfo(1+2*q,m))
        enddo
        ! Mortar index on this rank's faces of the mortar (big face, then sub-faces).
        call MarkMortarFace(globalMesh%mortarInfo(1,m),globalMesh%mortarInfo(2,m),mortarNew(m))
        do q = 1,4
          call MarkMortarFace(globalMesh%mortarInfo(1+2*q,m), &
                              globalMesh%mortarInfo(2+2*q,m)/10,mortarNew(m))
        enddo
      enddo
    endif

    subMesh%nAnisoMortars = na
    if(na > 0) then
      allocate(subMesh%anisoMortarInfo(1:9,1:na))
      do a = 1,globalMesh%nAnisoMortars
        if(anisoNew(a) == 0) cycle
        subMesh%anisoMortarInfo(1:9,anisoNew(a)) = globalMesh%anisoMortarInfo(1:9,a)
        subMesh%anisoMortarInfo(1,anisoNew(a)) = setGlobalId(globalMesh%anisoMortarInfo(1,a))
        do q = 1,2
          subMesh%anisoMortarInfo(1+2*q,anisoNew(a)) = &
            setGlobalId(globalMesh%anisoMortarInfo(1+2*q,a))
        enddo
        ! Anisotropic mortar faces carry sideInfo(1) = nMortars + a.
        call MarkMortarFace(globalMesh%anisoMortarInfo(1,a),globalMesh%anisoMortarInfo(2,a), &
                            nm+anisoNew(a))
        do q = 1,2
          call MarkMortarFace(globalMesh%anisoMortarInfo(1+2*q,a), &
                              globalMesh%anisoMortarInfo(2+2*q,a)/10,nm+anisoNew(a))
        enddo
      enddo
    endif

    deallocate(subLocal,mortarNew,anisoNew)

    call subMesh%UpdateDevice()

  contains

    subroutine MarkMortarFace(eGlobal,face,mortarIndex)
      !! Set sideInfo(1) = mortarIndex on local side face of global element eGlobal when this
      !! rank owns it; elements on other ranks are marked by their owners.
      implicit none
      integer,intent(in) :: eGlobal
      integer,intent(in) :: face
      integer,intent(in) :: mortarIndex
      ! Local
      integer :: elLocal

      elLocal = eGlobal-goff
      if(elLocal < 1 .or. elLocal > globalMesh%nElem) return
      subMesh%sideInfo(1,face,subLocal(elLocal)) = mortarIndex

    endsubroutine MarkMortarFace

  endsubroutine ExtractSubMesh