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, §ion);
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, §ionSF);
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(§ionSF);
84: PetscSectionDestroy(§ion);
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*/