1: const char help[] = "Set up a PetscSF for halo exchange between local vectors";

  3: #include <petscdmplex.h>
  4: #include <petscsf.h>

  6: int main(int argc, char **argv)
  7: {
  8:   DM           dm;
  9:   Vec          u;
 10:   PetscSection section;
 11:   PetscInt     dim, numFields, numBC, i;
 12:   PetscMPIInt  rank;
 13:   PetscInt     numComp[2];
 14:   PetscInt     numDof[12];
 15:   PetscInt    *remoteOffsets;
 16:   PetscSF      pointSF;
 17:   PetscSF      sectionSF;
 18:   PetscScalar *array;

 20:   PetscFunctionBeginUser;
 21:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
 22:   /* Create a mesh */
 23:   PetscCall(DMCreate(PETSC_COMM_WORLD, &dm));
 24:   PetscCall(DMSetType(dm, DMPLEX));
 25:   PetscCall(DMSetFromOptions(dm));
 26:   PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
 27:   PetscCall(DMGetDimension(dm, &dim));

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

 42:   /** Create a PetscSection to handle the layout of the discretized variables */
 43:   PetscCall(DMSetNumFields(dm, numFields));
 44:   PetscCall(DMPlexCreateSection(dm, NULL, numComp, numDof, numBC, NULL, NULL, NULL, NULL, &section));
 45:   /* Name the Field variables */
 46:   PetscCall(PetscSectionSetFieldName(section, 0, "u"));
 47:   PetscCall(PetscSectionSetFieldName(section, 1, "v"));
 48:   /* Tell the DM to use this data layout */
 49:   PetscCall(DMSetLocalSection(dm, section));

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

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

 82:   /** Cleanup */
 83:   PetscCall(PetscSFDestroy(&sectionSF));
 84:   PetscCall(PetscSectionDestroy(&section));
 85:   PetscCall(DMDestroy(&dm));
 86:   PetscCall(PetscFinalize());
 87:   return 0;
 88: }

 90: /*TEST

 92:   # Test on a 1D mesh with overlap
 93:   test:
 94:     nsize: 3
 95:     requires: !complex
 96:     args: -dm_plex_dim 1 -dm_plex_box_faces 3 -dm_refine_pre 1 -petscpartitioner_type simple -dm_distribute_overlap 1

 98: TEST*/