GatherSetOfElem Subroutine

  subroutine GatherSetOfElem(globalMesh,setOfElemLocal,setOfElem)
    !! Gather the set id of every global element onto every rank.
    !!
    !! On return setOfElem(e) = s(e) for e = 1..nGlobalElem, in global element order,
    !! where rank r contributed setOfElemLocal(1:nElem) for its own global range. This
    !! is one MPI_Allgatherv over globalMesh%decomp%mpiComm, run at setup only.
    !!
    !!  Input
    !!    - globalMesh : the decomposed global mesh
    !!    - setOfElemLocal(1:globalMesh%nElem) : set id of each rank-local element
    !!  Output
    !!    - setOfElem(1:nGlobalElem) : set id of each global element, where
    !!      nGlobalElem = globalMesh%decomp%offsetElem(nRanks+1)
    implicit none
    type(Mesh3D),intent(in) :: globalMesh
    integer,intent(in) :: setOfElemLocal(1:globalMesh%nElem)
    integer,intent(out) :: setOfElem(1:)
    ! Local
    integer :: r,nRanks,ierror
    integer,allocatable :: counts(:),displs(:)

    nRanks = globalMesh%decomp%nRanks
    if(size(setOfElem) /= globalMesh%decomp%offsetElem(nRanks+1)) then
      error stop "GatherSetOfElem : setOfElem must have one entry per global element"
    endif

    if(nRanks == 1) then
      setOfElem(1:globalMesh%nElem) = setOfElemLocal(1:globalMesh%nElem)
      return
    endif

    allocate(counts(1:nRanks),displs(1:nRanks))
    do r = 1,nRanks
      counts(r) = globalMesh%decomp%offsetElem(r+1)-globalMesh%decomp%offsetElem(r)
      displs(r) = globalMesh%decomp%offsetElem(r)
    enddo
    call mpi_allgatherv(setOfElemLocal,counts(globalMesh%decomp%rankId+1),MPI_INTEGER, &
                        setOfElem,counts,displs,MPI_INTEGER, &
                        globalMesh%decomp%mpiComm,ierror)
    deallocate(counts,displs)

  endsubroutine GatherSetOfElem