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
 10:   PetscInt, parameter              :: numFields = 2, numBC = 0 ! No boundary conditions
 11:   PetscMPIInt                      :: rank
 12:   PetscInt, dimension(2)           :: numComp
 13:   PetscInt, dimension(12)          :: numDof
 14:   PetscInt, dimension(:), pointer  :: remoteOffsets
 15:   type(tPetscSF)                   :: pointSF
 16:   type(tPetscSF)                   :: sectionSF
 17:   PetscScalar, dimension(:), pointer :: array
 18:   PetscReal                        :: val
 19:   PetscErrorCode                   :: ierr

 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:   numComp = [1_PETSC_INT_KIND, dim]
 31:   numDof = 0
 32:   !Let u be defined on cells
 33:   numDof(0*(dim + 1) + dim + 1) = 1
 34:   !Let v be defined on vertices
 35:   numDof(1*(dim + 1) + 1) = dim

 37:   !!! Create a PetscSection to handle the layout of the discretized variables
 38:   PetscCallA(DMSetNumFields(dm, numFields, ierr))
 39:   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))
 40:   ! Name the Field variables
 41:   PetscCallA(PetscSectionSetFieldName(section, 0_PETSC_INT_KIND, "u", ierr))
 42:   PetscCallA(PetscSectionSetFieldName(section, 1_PETSC_INT_KIND, "v", ierr))
 43:   ! Tell the DM to use this data layout
 44:   PetscCallA(DMSetLocalSection(dm, section, ierr))

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

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

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