Actual source code: ex14.c

  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;

 21:   PetscInitialize(&argc, &argv, NULL, help);
 22:   /* Create a mesh */
 23:   DMCreate(PETSC_COMM_WORLD, &dm);
 24:   DMSetType(dm, DMPLEX);
 25:   DMSetFromOptions(dm);
 26:   DMViewFromOptions(dm, NULL, "-dm_view");
 27:   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:   DMSetNumFields(dm, numFields);
 44:   DMPlexCreateSection(dm, NULL, numComp, numDof, numBC, NULL, NULL, NULL, NULL, &section);
 45:   /* Name the Field variables */
 46:   PetscSectionSetFieldName(section, 0, "u");
 47:   PetscSectionSetFieldName(section, 1, "v");
 48:   /* Tell the DM to use this data layout */
 49:   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:   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:   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:   PetscSFCreateSectionSF(pointSF, section, remoteOffsets, section, &sectionSF);
 61:   PetscFree(remoteOffsets);

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

 82:   /** Cleanup */
 83:   PetscSFDestroy(&sectionSF);
 84:   PetscSectionDestroy(&section);
 85:   DMDestroy(&dm);
 86:   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*/