1 static char help[] = "Tests for point location\n\n";
2
3 #include <petscsf.h>
4 #include <petscdmplex.h>
5
6 /* To inspect the location process, use
7
8 -dm_plex_print_locate 5 -info :dm
9 */
10
11 typedef struct {
12 PetscBool centroids;
13 PetscBool custom;
14 } AppCtx;
15
ProcessOptions(MPI_Comm comm,AppCtx * options)16 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
17 {
18 PetscFunctionBeginUser;
19 options->centroids = PETSC_TRUE;
20 options->custom = PETSC_FALSE;
21
22 PetscOptionsBegin(comm, "", "Point Location Options", "DMPLEX");
23 PetscCall(PetscOptionsBool("-centroids", "Locate cell centroids", "ex17.c", options->centroids, &options->centroids, NULL));
24 PetscCall(PetscOptionsBool("-custom", "Locate user-defined points", "ex17.c", options->custom, &options->custom, NULL));
25 PetscOptionsEnd();
26 PetscFunctionReturn(PETSC_SUCCESS);
27 }
28
CreateMesh(MPI_Comm comm,DM * dm)29 static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm)
30 {
31 PetscFunctionBeginUser;
32 PetscCall(DMCreate(comm, dm));
33 PetscCall(DMSetType(*dm, DMPLEX));
34 PetscCall(DMSetFromOptions(*dm));
35 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
36 PetscFunctionReturn(PETSC_SUCCESS);
37 }
38
TestCentroidLocation(DM dm,AppCtx * user)39 static PetscErrorCode TestCentroidLocation(DM dm, AppCtx *user)
40 {
41 Vec points;
42 PetscSF cellSF = NULL;
43 const PetscSFNode *cells;
44 PetscScalar *a;
45 PetscInt cdim, n;
46 PetscInt cStart, cEnd, c;
47
48 PetscFunctionBeginUser;
49 if (!user->centroids) PetscFunctionReturn(PETSC_SUCCESS);
50 PetscCall(DMGetCoordinateDim(dm, &cdim));
51 PetscCall(DMGetCoordinatesLocalSetUp(dm));
52 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
53 /* Locate all centroids */
54 PetscCall(VecCreateSeq(PETSC_COMM_SELF, (cEnd - cStart) * cdim, &points));
55 PetscCall(VecSetBlockSize(points, cdim));
56 PetscCall(VecGetArray(points, &a));
57 for (c = cStart; c < cEnd; ++c) {
58 PetscReal centroid[3];
59 PetscInt off = (c - cStart) * cdim, d;
60
61 PetscCall(DMPlexComputeCellGeometryFVM(dm, c, NULL, centroid, NULL));
62 for (d = 0; d < cdim; ++d) a[off + d] = centroid[d];
63 }
64 PetscCall(VecRestoreArray(points, &a));
65 PetscCall(DMLocatePoints(dm, points, DM_POINTLOCATION_NONE, &cellSF));
66 PetscCall(VecDestroy(&points));
67 PetscCall(PetscSFGetGraph(cellSF, NULL, &n, NULL, &cells));
68 if (n != (cEnd - cStart)) {
69 for (c = 0; c < n; ++c) {
70 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));
71 }
72 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Located %" PetscInt_FMT " points instead of %" PetscInt_FMT, n, cEnd - cStart);
73 }
74 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);
75 PetscCall(PetscSFDestroy(&cellSF));
76 PetscFunctionReturn(PETSC_SUCCESS);
77 }
78
TestCustomLocation(DM dm,AppCtx * user)79 static PetscErrorCode TestCustomLocation(DM dm, AppCtx *user)
80 {
81 PetscSF cellSF = NULL;
82 const PetscSFNode *cells;
83 const PetscInt *found;
84 Vec points;
85 PetscScalar coords[2] = {0.5, 0.5};
86 PetscInt cdim, Np = 1, Nfd;
87 PetscMPIInt rank;
88 MPI_Comm comm;
89
90 PetscFunctionBeginUser;
91 if (!user->custom) PetscFunctionReturn(PETSC_SUCCESS);
92 PetscCall(DMGetCoordinateDim(dm, &cdim));
93
94 // Locate serially on each process
95 PetscCall(VecCreate(PETSC_COMM_SELF, &points));
96 PetscCall(VecSetBlockSize(points, cdim));
97 PetscCall(VecSetSizes(points, Np * cdim, PETSC_DETERMINE));
98 PetscCall(VecSetFromOptions(points));
99 for (PetscInt p = 0; p < Np; ++p) {
100 const PetscInt idx[2] = {p * cdim, p * cdim + 1};
101 PetscCall(VecSetValues(points, cdim, idx, coords, INSERT_VALUES));
102 }
103 PetscCall(VecAssemblyBegin(points));
104 PetscCall(VecAssemblyEnd(points));
105
106 PetscCall(DMLocatePoints(dm, points, DM_POINTLOCATION_NONE, &cellSF));
107
108 PetscCall(PetscSFGetGraph(cellSF, NULL, &Nfd, &found, &cells));
109 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
110 PetscCallMPI(MPI_Comm_rank(comm, &rank));
111 PetscCall(PetscSynchronizedPrintf(comm, "[%d] Found %" PetscInt_FMT " particles\n", rank, Nfd));
112 for (PetscInt p = 0; p < Nfd; ++p) {
113 const PetscInt point = found ? found[p] : p;
114 const PetscScalar *array;
115 PetscScalar *ccoords = NULL;
116 PetscInt numCoords;
117 PetscBool isDG;
118
119 // Since the v comm is SELF, rank is always 0
120 PetscCall(PetscSynchronizedPrintf(comm, " point %" PetscInt_FMT " cell %" PetscInt_FMT "\n", point, cells[p].index));
121 PetscCall(DMPlexGetCellCoordinates(dm, cells[p].index, &isDG, &numCoords, &array, &ccoords));
122 for (PetscInt c = 0; c < numCoords / cdim; ++c) {
123 PetscCall(PetscSynchronizedPrintf(comm, " "));
124 for (PetscInt d = 0; d < cdim; ++d) PetscCall(PetscSynchronizedPrintf(comm, " %g", (double)PetscRealPart(ccoords[c * cdim + d])));
125 PetscCall(PetscSynchronizedPrintf(comm, "\n"));
126 }
127 PetscCall(DMPlexRestoreCellCoordinates(dm, cells[p].index, &isDG, &numCoords, &array, &ccoords));
128 }
129 PetscCall(PetscSynchronizedFlush(comm, PETSC_STDOUT));
130
131 PetscCall(PetscSFDestroy(&cellSF));
132 PetscCall(VecDestroy(&points));
133 PetscFunctionReturn(PETSC_SUCCESS);
134 }
135
main(int argc,char ** argv)136 int main(int argc, char **argv)
137 {
138 DM dm;
139 AppCtx user;
140
141 PetscFunctionBeginUser;
142 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
143 PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
144 PetscCall(CreateMesh(PETSC_COMM_WORLD, &dm));
145 PetscCall(TestCentroidLocation(dm, &user));
146 PetscCall(TestCustomLocation(dm, &user));
147 PetscCall(DMDestroy(&dm));
148 PetscCall(PetscFinalize());
149 return 0;
150 }
151
152 /*TEST
153
154 testset:
155 args: -dm_plex_dim 1 -dm_plex_box_faces 10
156 output_file: output/empty.out
157
158 test:
159 suffix: seg
160
161 test:
162 suffix: seg_hash
163 args: -dm_refine 2 -dm_plex_hash_location
164
165 testset:
166 args: -dm_plex_box_faces 5,5
167 output_file: output/empty.out
168
169 test:
170 suffix: tri
171 requires: triangle
172
173 test:
174 suffix: tri_hash
175 requires: triangle
176 args: -dm_refine 2 -dm_plex_hash_location
177
178 test:
179 suffix: quad
180 args: -dm_plex_simplex 0
181
182 test:
183 suffix: quad_order_2
184 args: -dm_plex_simplex 0 -dm_coord_petscspace_degree 2
185
186 test:
187 suffix: quad_hash
188 args: -dm_plex_simplex 0 -dm_refine 2 -dm_plex_hash_location
189
190 testset:
191 args: -dm_plex_dim 3 -dm_plex_box_faces 3,3,3
192 output_file: output/empty.out
193
194 test:
195 suffix: tet
196 requires: ctetgen
197
198 test:
199 suffix: tet_hash
200 requires: ctetgen
201 args: -dm_refine 1 -dm_plex_hash_location
202
203 test:
204 suffix: hex
205 args: -dm_plex_simplex 0
206
207 test:
208 suffix: hex_hash
209 args: -dm_plex_simplex 0 -dm_refine 1 -dm_plex_hash_location
210
211 test:
212 suffix: hex_order_2
213 args: -dm_plex_simplex 0 -dm_refine 1 -dm_coord_petscspace_degree 2
214 nsize: 2
215
216 test:
217 suffix: hex_order_3
218 args: -dm_plex_simplex 0 -dm_coord_petscspace_degree 3
219
220 testset:
221 args: -centroids 0 -custom \
222 -dm_plex_simplex 0 -dm_plex_box_faces 21,21 -dm_distribute_overlap 4 -petscpartitioner_type simple
223 nsize: 2
224
225 test:
226 suffix: quad_overlap
227 args: -dm_plex_hash_location {{0 1}}
228
229 # Test location on a Monge Manifold
230 testset:
231 args: -dm_refine 3 -dm_coord_space 0 \
232 -dm_plex_option_phases proj_ -cdm_proj_dm_plex_coordinate_dim 3 -proj_dm_coord_space \
233 -proj_dm_coord_remap -proj_dm_coord_map sinusoid -proj_dm_coord_map_params 0.1,1.,1.
234 output_file: output/empty.out
235
236 test:
237 requires: triangle
238 suffix: tri_monge
239
240 test:
241 suffix: quad_monge
242 args: -dm_plex_simplex 0
243
244 TEST*/
245