1: #include <petsc/finclude/petsc.h>
  2: program ex14f90
  3:   use petsc
  4:   implicit none

  6:   type(tDM)                        :: dm
  7:   type(tVec)                       :: u
  8:   type(tPetscSection)              :: section
  9:   PetscInt                         :: dim, numFields, numBC
 10:   PetscMPIInt                      :: rank
 11:   PetscInt, dimension(2)            :: numComp
 12:   PetscInt, dimension(12)           :: numDof
 13:   PetscInt, dimension(:), pointer    :: remoteOffsets
 14:   type(tPetscSF)                   :: pointSF
 15:   type(tPetscSF)                   :: sectionSF
 16:   PetscScalar, dimension(:), pointer :: array
 17:   PetscReal                        :: val
 18:   PetscErrorCode                   :: ierr
 19:   PetscInt                         :: zero = 0, one = 1

 21:   PetscCallA(PetscInitialize(PETSC_NULL_CHARACTER, ierr))
 22:   PetscCallA(DMCreate(PETSC_COMM_WORLD, dm, ierr))
 23:   PetscCallA(DMSetType(dm, DMPLEX, ierr))
 24:   PetscCallA(DMSetFromOptions(dm, ierr))
 25:   PetscCallA(DMViewFromOptions(dm, PETSC_NuLL_OBJECT, "-dm_view", ierr))
 26:   PetscCallA(DMGetDimension(dm, dim, ierr))

 28:   !!! Describe the solution variables that are discretized on the mesh
 29:   ! Create scalar field u and a vector field v
 30:   numFields = 2
 31:   numComp = [one, dim]
 32:   numDof = 0
 33:   !Let u be defined on cells
 34:   numDof(0*(dim + 1) + dim + 1) = 1
 35:   !Let v be defined on vertices
 36:   numDof(1*(dim + 1) + 1) = dim
 37:   !No boundary conditions */
 38:   numBC = 0

 40:   !!! Create a PetscSection to handle the layout of the discretized variables
 41:   PetscCallA(DMSetNumFields(dm, numFields, ierr))
 42:   PetscCallA(DMPlexCreateSection(dm, PETSC_NULL_DMLABEL_ARRAY, numComp, numDof, numBC, PETSC_NULL_INTEGER_ARRAY, PETSC_NULL_IS_ARRAY, PETSC_NULL_IS_ARRAY, PETSC_NULL_IS, section, ierr))
 43:   ! Name the Field variables
 44:   PetscCallA(PetscSectionSetFieldName(section, zero, "u", ierr))
 45:   PetscCallA(PetscSectionSetFieldName(section, one, "v", ierr))
 46:   ! Tell the DM to use this data layout
 47:   PetscCallA(DMSetLocalSection(dm, section, ierr))

 49:   !!! Construct the communication pattern for halo exchange between local vectors */
 50:   ! Get the point SF: an object that says which copies of mesh points (cells,
 51:   ! vertices, faces, edges) are copies of points on other processes
 52:   PetscCallA(DMGetPointSF(dm, pointSF, ierr))
 53:   ! Relate the locations of ghost degrees of freedom on this process
 54:   ! to their locations of the non-ghost copies on a different process
 55:   PetscCallA(PetscSFCreateRemoteOffsets(pointSF, section, section, remoteOffsets, ierr))
 56:   ! Use that information to construct a star forest for halo exchange
 57:   ! for data described by the local section
 58:   PetscCallA(PetscSFCreateSectionSF(pointSF, section, remoteOffsets, section, sectionSF, ierr))
 59:   if (associated(remoteOffsets)) then
 60:     PetscCallA(PetscSFDestroyRemoteOffsets(remoteOffsets, ierr))
 61:   end if

 63:   !!! Demo of halo exchange
 64:   ! Create a Vec with this layout
 65:   PetscCallA(DMCreateLocalVector(dm, u, ierr))
 66:   PetscCallA(PetscObjectSetName(u, "Local vector", ierr))
 67:   PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))
 68:   ! Set all mesh values to the MPI rank
 69:   val = rank
 70:   PetscCallA(VecSet(u, val, ierr))
 71:   ! Get the raw array of values
 72:   PetscCallA(VecGetArrayWrite(u, array, ierr))
 73:   !!! HALO EXCHANGE
 74:   PetscCallA(PetscSFBcastBegin(sectionSF, MPIU_SCALAR, array, array, MPI_REPLACE, ierr))
 75:   ! local work can be done between Begin() and End()
 76:   PetscCallA(PetscSFBcastEnd(sectionSF, MPIU_SCALAR, array, array, MPI_REPLACE, ierr))
 77:   ! Restore the raw array of values
 78:   PetscCallA(VecRestoreArrayWrite(u, array, ierr))
 79:   ! View the results: should show which process has the non-ghost copy of each degree of freedom
 80:   PetscCallA(PetscSectionVecView(section, u, PETSC_VIEWER_STDOUT_WORLD, ierr))
 81:   PetscCallA(VecDestroy(u, ierr))

 83:   PetscCallA(PetscSFDestroy(sectionSF, ierr))
 84:   PetscCallA(PetscSectionDestroy(section, ierr))
 85:   PetscCallA(DMDestroy(dm, ierr))
 86:   PetscCallA(PetscFinalize(ierr))
 87: end program ex14f90
 88: !/*TEST
 89: !  # Test on a 1D mesh with overlap
 90: !  test:
 91: !    nsize: 3
 92: !    requires: !complex
 93: !    args: -dm_plex_dim 1 -dm_plex_box_faces 3 -dm_refine_pre 1 -petscpartitioner_type simple -dm_distribute_overlap 1
 94: !
 95: !TEST*/