Actual source code: ex17.c
1: static char help[] = "Tests for point location\n\n";
3: #include <petscsf.h>
4: #include <petscdmplex.h>
6: typedef struct {
7: PetscBool centroids;
8: PetscBool custom;
9: } AppCtx;
11: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
12: {
13: PetscFunctionBeginUser;
14: options->centroids = PETSC_TRUE;
15: options->custom = PETSC_FALSE;
17: PetscOptionsBegin(comm, "", "Point Location Options", "DMPLEX");
18: PetscCall(PetscOptionsBool("-centroids", "Locate cell centroids", "ex17.c", options->centroids, &options->centroids, NULL));
19: PetscCall(PetscOptionsBool("-custom", "Locate user-defined points", "ex17.c", options->custom, &options->custom, NULL));
20: PetscOptionsEnd();
21: PetscFunctionReturn(PETSC_SUCCESS);
22: }
24: static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm)
25: {
26: PetscFunctionBeginUser;
27: PetscCall(DMCreate(comm, dm));
28: PetscCall(DMSetType(*dm, DMPLEX));
29: PetscCall(DMSetFromOptions(*dm));
30: PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
31: PetscFunctionReturn(PETSC_SUCCESS);
32: }
34: static PetscErrorCode TestCentroidLocation(DM dm, AppCtx *user)
35: {
36: Vec points;
37: PetscSF cellSF = NULL;
38: const PetscSFNode *cells;
39: PetscScalar *a;
40: PetscInt cdim, n;
41: PetscInt cStart, cEnd, c;
43: PetscFunctionBeginUser;
44: if (!user->centroids) PetscFunctionReturn(PETSC_SUCCESS);
45: PetscCall(DMGetCoordinateDim(dm, &cdim));
46: PetscCall(DMGetCoordinatesLocalSetUp(dm));
47: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
48: /* Locate all centroids */
49: PetscCall(VecCreateSeq(PETSC_COMM_SELF, (cEnd - cStart) * cdim, &points));
50: PetscCall(VecSetBlockSize(points, cdim));
51: PetscCall(VecGetArray(points, &a));
52: for (c = cStart; c < cEnd; ++c) {
53: PetscReal centroid[3];
54: PetscInt off = (c - cStart) * cdim, d;
56: PetscCall(DMPlexComputeCellGeometryFVM(dm, c, NULL, centroid, NULL));
57: for (d = 0; d < cdim; ++d) a[off + d] = centroid[d];
58: }
59: PetscCall(VecRestoreArray(points, &a));
60: PetscCall(DMLocatePoints(dm, points, DM_POINTLOCATION_NONE, &cellSF));
61: PetscCall(VecDestroy(&points));
62: PetscCall(PetscSFGetGraph(cellSF, NULL, &n, NULL, &cells));
63: if (n != (cEnd - cStart)) {
64: for (c = 0; c < n; ++c) {
65: if (cells[c].index != c + cStart) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Could not locate centroid of cell %" PetscInt_FMT ", error %" PetscInt_FMT "\n", c + cStart, cells[c].index));
66: }
67: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Located %" PetscInt_FMT " points instead of %" PetscInt_FMT, n, cEnd - cStart);
68: }
69: for (c = cStart; c < cEnd; ++c) PetscCheck(cells[c - cStart].index == c, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not locate centroid of cell %" PetscInt_FMT ", instead found %" PetscInt_FMT, c, cells[c - cStart].index);
70: PetscCall(PetscSFDestroy(&cellSF));
71: PetscFunctionReturn(PETSC_SUCCESS);
72: }
74: static PetscErrorCode TestCustomLocation(DM dm, AppCtx *user)
75: {
76: PetscSF cellSF = NULL;
77: const PetscSFNode *cells;
78: const PetscInt *found;
79: Vec points;
80: PetscScalar coords[2] = {0.5, 0.5};
81: PetscInt cdim, Np = 1, Nfd;
82: PetscMPIInt rank;
83: MPI_Comm comm;
85: PetscFunctionBeginUser;
86: if (!user->custom) PetscFunctionReturn(PETSC_SUCCESS);
87: PetscCall(DMGetCoordinateDim(dm, &cdim));
89: // Locate serially on each process
90: PetscCall(VecCreate(PETSC_COMM_SELF, &points));
91: PetscCall(VecSetBlockSize(points, cdim));
92: PetscCall(VecSetSizes(points, Np * cdim, PETSC_DETERMINE));
93: PetscCall(VecSetFromOptions(points));
94: for (PetscInt p = 0; p < Np; ++p) {
95: const PetscInt idx[2] = {p * cdim, p * cdim + 1};
96: PetscCall(VecSetValues(points, cdim, idx, coords, INSERT_VALUES));
97: }
98: PetscCall(VecAssemblyBegin(points));
99: PetscCall(VecAssemblyEnd(points));
101: PetscCall(DMLocatePoints(dm, points, DM_POINTLOCATION_NONE, &cellSF));
103: PetscCall(PetscSFGetGraph(cellSF, NULL, &Nfd, &found, &cells));
104: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
105: PetscCallMPI(MPI_Comm_rank(comm, &rank));
106: PetscCall(PetscSynchronizedPrintf(comm, "[%d] Found %" PetscInt_FMT " particles\n", rank, Nfd));
107: for (PetscInt p = 0; p < Nfd; ++p) {
108: const PetscInt point = found ? found[p] : p;
109: const PetscScalar *array;
110: PetscScalar *ccoords = NULL;
111: PetscInt numCoords;
112: PetscBool isDG;
114: // Since the v comm is SELF, rank is always 0
115: PetscCall(PetscSynchronizedPrintf(comm, " point %" PetscInt_FMT " cell %" PetscInt_FMT "\n", point, cells[p].index));
116: PetscCall(DMPlexGetCellCoordinates(dm, cells[p].index, &isDG, &numCoords, &array, &ccoords));
117: for (PetscInt c = 0; c < numCoords / cdim; ++c) {
118: PetscCall(PetscSynchronizedPrintf(comm, " "));
119: for (PetscInt d = 0; d < cdim; ++d) PetscCall(PetscSynchronizedPrintf(comm, " %g", (double)PetscRealPart(ccoords[c * cdim + d])));
120: PetscCall(PetscSynchronizedPrintf(comm, "\n"));
121: }
122: PetscCall(DMPlexRestoreCellCoordinates(dm, cells[p].index, &isDG, &numCoords, &array, &ccoords));
123: }
124: PetscCall(PetscSynchronizedFlush(comm, PETSC_STDOUT));
126: PetscCall(PetscSFDestroy(&cellSF));
127: PetscCall(VecDestroy(&points));
128: PetscFunctionReturn(PETSC_SUCCESS);
129: }
131: int main(int argc, char **argv)
132: {
133: DM dm;
134: AppCtx user;
136: PetscFunctionBeginUser;
137: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
138: PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
139: PetscCall(CreateMesh(PETSC_COMM_WORLD, &dm));
140: PetscCall(TestCentroidLocation(dm, &user));
141: PetscCall(TestCustomLocation(dm, &user));
142: PetscCall(DMDestroy(&dm));
143: PetscCall(PetscFinalize());
144: return 0;
145: }
147: /*TEST
149: testset:
150: args: -dm_plex_dim 1 -dm_plex_box_faces 10
152: test:
153: suffix: seg
155: test:
156: suffix: seg_hash
157: args: -dm_refine 2 -dm_plex_hash_location
159: testset:
160: args: -dm_plex_box_faces 5,5
162: test:
163: suffix: tri
164: requires: triangle
166: test:
167: suffix: tri_hash
168: requires: triangle
169: args: -dm_refine 2 -dm_plex_hash_location
171: test:
172: suffix: quad
173: args: -dm_plex_simplex 0
175: test:
176: suffix: quad_hash
177: args: -dm_plex_simplex 0 -dm_refine 2 -dm_plex_hash_location
179: testset:
180: args: -dm_plex_dim 3 -dm_plex_box_faces 3,3,3
182: test:
183: suffix: tet
184: requires: ctetgen
186: test:
187: suffix: tet_hash
188: requires: ctetgen
189: args: -dm_refine 1 -dm_plex_hash_location
191: test:
192: suffix: hex
193: args: -dm_plex_simplex 0
195: test:
196: suffix: hex_hash
197: args: -dm_plex_simplex 0 -dm_refine 1 -dm_plex_hash_location
199: testset:
200: args: -centroids 0 -custom \
201: -dm_plex_simplex 0 -dm_plex_box_faces 21,21 -dm_distribute_overlap 4 -petscpartitioner_type simple
202: nsize: 2
204: test:
205: suffix: quad_overlap
206: args: -dm_plex_hash_location {{0 1}}
208: TEST*/