1 #include <petscdmda.h> /*I "petscdmda.h" I*/
2 #include <petscdmplex.h> /*I "petscdmplex.h" I*/
3 #include <petsc/private/dmswarmimpl.h> /*I "petscdmswarm.h" I*/
4
sort_CompareSwarmPoint(const void * dataA,const void * dataB)5 static int sort_CompareSwarmPoint(const void *dataA, const void *dataB)
6 {
7 SwarmPoint *pointA = (SwarmPoint *)dataA;
8 SwarmPoint *pointB = (SwarmPoint *)dataB;
9
10 if (pointA->cell_index < pointB->cell_index) {
11 return -1;
12 } else if (pointA->cell_index > pointB->cell_index) {
13 return 1;
14 } else {
15 return 0;
16 }
17 }
18
DMSwarmSortApplyCellIndexSort(DMSwarmSort ctx)19 static PetscErrorCode DMSwarmSortApplyCellIndexSort(DMSwarmSort ctx)
20 {
21 PetscFunctionBegin;
22 if (ctx->list) qsort(ctx->list, ctx->npoints, sizeof(SwarmPoint), sort_CompareSwarmPoint);
23 PetscFunctionReturn(PETSC_SUCCESS);
24 }
25
DMSwarmSortCreate(DMSwarmSort * _ctx)26 static PetscErrorCode DMSwarmSortCreate(DMSwarmSort *_ctx)
27 {
28 DMSwarmSort ctx;
29
30 PetscFunctionBegin;
31 PetscCall(PetscNew(&ctx));
32 ctx->isvalid = PETSC_FALSE;
33 ctx->ncells = 0;
34 ctx->npoints = 0;
35 PetscCall(PetscMalloc1(1, &ctx->pcell_offsets));
36 PetscCall(PetscMalloc1(1, &ctx->list));
37 *_ctx = ctx;
38 PetscFunctionReturn(PETSC_SUCCESS);
39 }
40
DMSwarmSortSetup(DMSwarmSort ctx,DM dm,PetscInt ncells)41 static PetscErrorCode DMSwarmSortSetup(DMSwarmSort ctx, DM dm, PetscInt ncells)
42 {
43 DMSwarmCellDM celldm;
44 PetscInt *swarm_cellid;
45 PetscInt p, npoints;
46 PetscInt tmp, c, count;
47 const char *cellid;
48
49 PetscFunctionBegin;
50 if (!ctx) PetscFunctionReturn(PETSC_SUCCESS);
51 if (ctx->isvalid) PetscFunctionReturn(PETSC_SUCCESS);
52
53 PetscCall(PetscLogEventBegin(DMSWARM_Sort, 0, 0, 0, 0));
54 /* check the number of cells */
55 if (ncells != ctx->ncells) {
56 PetscCall(PetscRealloc(sizeof(PetscInt) * (ncells + 1), &ctx->pcell_offsets));
57 ctx->ncells = ncells;
58 }
59 PetscCall(PetscArrayzero(ctx->pcell_offsets, ctx->ncells + 1));
60
61 /* get the number of points */
62 PetscCall(DMSwarmGetLocalSize(dm, &npoints));
63 if (npoints != ctx->npoints) {
64 PetscCall(PetscRealloc(sizeof(SwarmPoint) * npoints, &ctx->list));
65 ctx->npoints = npoints;
66 }
67 PetscCall(PetscArrayzero(ctx->list, npoints));
68
69 PetscCall(DMSwarmGetCellDMActive(dm, &celldm));
70 PetscCall(DMSwarmCellDMGetCellID(celldm, &cellid));
71 PetscCall(DMSwarmGetField(dm, cellid, NULL, NULL, (void **)&swarm_cellid));
72 for (p = 0; p < ctx->npoints; p++) {
73 ctx->list[p].point_index = p;
74 ctx->list[p].cell_index = swarm_cellid[p];
75 }
76 PetscCall(DMSwarmRestoreField(dm, cellid, NULL, NULL, (void **)&swarm_cellid));
77 PetscCall(DMSwarmSortApplyCellIndexSort(ctx));
78
79 /* sum points per cell */
80 for (p = 0; p < ctx->npoints; p++) ctx->pcell_offsets[ctx->list[p].cell_index]++;
81
82 /* create offset list */
83 count = 0;
84 for (c = 0; c < ctx->ncells; c++) {
85 tmp = ctx->pcell_offsets[c];
86 ctx->pcell_offsets[c] = count;
87 count = count + tmp;
88 }
89 ctx->pcell_offsets[c] = count;
90
91 ctx->isvalid = PETSC_TRUE;
92 PetscCall(PetscLogEventEnd(DMSWARM_Sort, 0, 0, 0, 0));
93 PetscFunctionReturn(PETSC_SUCCESS);
94 }
95
DMSwarmSortDestroy(DMSwarmSort * ctx)96 PetscErrorCode DMSwarmSortDestroy(DMSwarmSort *ctx)
97 {
98 DMSwarmSort ictx;
99
100 PetscFunctionBegin;
101 if (!ctx) PetscFunctionReturn(PETSC_SUCCESS);
102 if (!*ctx) PetscFunctionReturn(PETSC_SUCCESS);
103 ictx = *ctx;
104 if (ictx->list) PetscCall(PetscFree(ictx->list));
105 if (ictx->pcell_offsets) PetscCall(PetscFree(ictx->pcell_offsets));
106 PetscCall(PetscFree(ictx));
107 *ctx = NULL;
108 PetscFunctionReturn(PETSC_SUCCESS);
109 }
110
111 /*@
112 DMSwarmSortGetNumberOfPointsPerCell - Returns the number of points in a cell
113
114 Not Collective
115
116 Input Parameters:
117 + sw - a `DMSWARM` objects
118 - cell - the cell number in the cell `DM`
119
120 Output Parameter:
121 . npoints - the number of points in the cell
122
123 Level: advanced
124
125 Notes:
126 You must call `DMSwarmSortGetAccess()` before you can call `DMSwarmSortGetNumberOfPointsPerCell()`
127
128 .seealso: `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSortGetAccess()`, `DMSwarmSortGetPointsPerCell()`
129 @*/
DMSwarmSortGetNumberOfPointsPerCell(DM sw,PetscInt cell,PetscInt * npoints)130 PetscErrorCode DMSwarmSortGetNumberOfPointsPerCell(DM sw, PetscInt cell, PetscInt *npoints)
131 {
132 DMSwarmCellDM celldm;
133 DMSwarmSort ctx;
134
135 PetscFunctionBegin;
136 PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
137 PetscCall(DMSwarmCellDMGetSort(celldm, &ctx));
138 PetscCheck(ctx, PetscObjectComm((PetscObject)sw), PETSC_ERR_USER, "The DMSwarmSort context has not been created. Must call DMSwarmSortGetAccess() first");
139 PetscCheck(ctx->isvalid, PETSC_COMM_SELF, PETSC_ERR_USER, "SwarmPointSort container is not valid. Must call DMSwarmSortGetAccess() first");
140 PetscCheck(cell < ctx->ncells, PETSC_COMM_SELF, PETSC_ERR_USER, "Cell index (%" PetscInt_FMT ") is greater than max number of local cells (%" PetscInt_FMT ")", cell, ctx->ncells);
141 PetscCheck(cell >= 0, PETSC_COMM_SELF, PETSC_ERR_USER, "Cell index (%" PetscInt_FMT ") cannot be negative", cell);
142 *npoints = ctx->pcell_offsets[cell + 1] - ctx->pcell_offsets[cell];
143 PetscFunctionReturn(PETSC_SUCCESS);
144 }
145
146 /*@C
147 DMSwarmSortGetPointsPerCell - Creates an array of point indices for all points in a cell
148
149 Not Collective
150
151 Input Parameters:
152 + sw - a `DMSWARM` object
153 . cell - the cell number in the cell `DM`
154 . npoints - the number of points in the cell
155 - pidlist - array of the indices identifying all points in cell e
156
157 Level: advanced
158
159 Note:
160 You must call `DMSwarmSortGetAccess()` before you can call `DMSwarmSortGetPointsPerCell()`, and call `DMSwarmRestorePointsPerCell()` afterwards
161
162 .seealso: `DMSWARM`, `DMSwarmSetType()`, `DMSwarmRestorePointsPerCell()`, `DMSwarmSortGetAccess()`, `DMSwarmSortGetNumberOfPointsPerCell()`
163 @*/
DMSwarmSortGetPointsPerCell(DM sw,PetscInt cell,PetscInt * npoints,PetscInt ** pidlist)164 PetscErrorCode DMSwarmSortGetPointsPerCell(DM sw, PetscInt cell, PetscInt *npoints, PetscInt **pidlist)
165 {
166 DMSwarmCellDM celldm;
167 PetscInt pid, pid_unsorted;
168 DMSwarmSort ctx;
169
170 PetscFunctionBegin;
171 PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
172 PetscCall(DMSwarmCellDMGetSort(celldm, &ctx));
173 PetscCheck(ctx, PetscObjectComm((PetscObject)sw), PETSC_ERR_USER, "The DMSwarmSort context has not been created. Must call DMSwarmSortGetAccess() first");
174 PetscCall(DMSwarmSortGetNumberOfPointsPerCell(sw, cell, npoints));
175 PetscCall(DMGetWorkArray(sw, *npoints, MPIU_SCALAR, pidlist));
176 for (PetscInt p = 0; p < *npoints; ++p) {
177 pid = ctx->pcell_offsets[cell] + p;
178 pid_unsorted = ctx->list[pid].point_index;
179 (*pidlist)[p] = pid_unsorted;
180 }
181 PetscFunctionReturn(PETSC_SUCCESS);
182 }
183
184 /*@C
185 DMSwarmSortRestorePointsPerCell - Restores an array of point indices for all points in a cell
186
187 Not Collective
188
189 Input Parameters:
190 + dm - a `DMSWARM` object
191 . e - the index of the cell
192 . npoints - the number of points in the cell
193 - pidlist - array of the indices identifying all points in cell e
194
195 Level: advanced
196
197 Note:
198 You must call `DMSwarmSortGetAccess()` and `DMSwarmSortGetPointsPerCell()` before you can call `DMSwarmSortRestorePointsPerCell()`
199
200 .seealso: `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSortGetPointsPerCell()`, `DMSwarmSortGetAccess()`, `DMSwarmSortGetNumberOfPointsPerCell()`
201 @*/
DMSwarmSortRestorePointsPerCell(DM dm,PetscInt e,PetscInt * npoints,PetscInt ** pidlist)202 PetscErrorCode DMSwarmSortRestorePointsPerCell(DM dm, PetscInt e, PetscInt *npoints, PetscInt **pidlist)
203 {
204 PetscFunctionBegin;
205 PetscCall(DMRestoreWorkArray(dm, *npoints, MPIU_SCALAR, pidlist));
206 PetscFunctionReturn(PETSC_SUCCESS);
207 }
208
209 /*@
210 DMSwarmSortGetAccess - Setups up a `DMSWARM` point sort context for efficient traversal of points within a cell
211
212 Not Collective
213
214 Input Parameter:
215 . sw - a `DMSWARM` object
216
217 Level: advanced
218
219 Notes:
220 Calling `DMSwarmSortGetAccess()` creates a list which enables easy identification of all points contained in a
221 given cell. This method does not explicitly sort the data within the `DMSWARM` based on the cell index associated
222 with a `DMSWARM` point.
223
224 The sort context is valid only for the `DMSWARM` points defined at the time when `DMSwarmSortGetAccess()` was called.
225 For example, suppose the swarm contained NP points when `DMSwarmSortGetAccess()` was called. If the user subsequently
226 adds 10 additional points to the swarm, the sort context is still valid, but only for the first NP points.
227 The indices associated with the 10 new additional points will not be contained within the sort context.
228 This means that the user can still safely perform queries via `DMSwarmSortGetPointsPerCell()` and
229 `DMSwarmSortGetPointsPerCell()`, however the results return will be based on the first NP points.
230
231 If any` DMSWARM` re-sizing method is called after `DMSwarmSortGetAccess()` which modifies any of the first NP entries
232 in the `DMSWARM`, the sort context will become invalid. Currently there are no guards to prevent the user from
233 invalidating the sort context. For this reason, we highly recommend you do not use `DMSwarmRemovePointAtIndex()` in
234 between calls to `DMSwarmSortGetAccess()` and `DMSwarmSortRestoreAccess()`.
235
236 To facilitate safe removal of points using the sort context, we suggest a "two pass" strategy in which the
237 first pass "marks" points for removal, and the second pass actually removes the points from the `DMSWARM`
238
239 You must call `DMSwarmSortGetAccess()` before you can call `DMSwarmSortGetPointsPerCell()` or `DMSwarmSortGetNumberOfPointsPerCell()`
240
241 The sort context may become invalid if any re-sizing methods are applied which alter the first NP points
242 within swarm at the time `DMSwarmSortGetAccess()` was called.
243
244 You must call `DMSwarmSortRestoreAccess()` when you no longer need access to the sort context
245
246 .seealso: `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSortRestoreAccess()`
247 @*/
DMSwarmSortGetAccess(DM sw)248 PetscErrorCode DMSwarmSortGetAccess(DM sw)
249 {
250 DM dm;
251 DMSwarmCellDM celldm;
252 DMSwarmSort ctx;
253 PetscInt ncells = 0;
254 PetscBool isda, isplex, isshell;
255
256 PetscFunctionBegin;
257 PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
258 PetscCall(DMSwarmCellDMGetSort(celldm, &ctx));
259 if (!ctx) {
260 PetscCall(DMSwarmSortCreate(&ctx));
261 PetscCall(DMSwarmCellDMSetSort(celldm, ctx));
262 }
263
264 /* get the number of cells */
265 PetscCall(DMSwarmGetCellDM(sw, &dm));
266 PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMDA, &isda));
267 PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMPLEX, &isplex));
268 PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMSHELL, &isshell));
269 if (isda) {
270 const PetscInt *element;
271 PetscInt nel, npe;
272
273 PetscCall(DMDAGetElements(dm, &nel, &npe, &element));
274 ncells = nel;
275 PetscCall(DMDARestoreElements(dm, &nel, &npe, &element));
276 } else if (isplex) {
277 PetscInt ps, pe;
278
279 PetscCall(DMPlexGetHeightStratum(dm, 0, &ps, &pe));
280 ncells = pe - ps;
281 } else if (isshell) {
282 PetscErrorCode (*method_DMShellGetNumberOfCells)(DM, PetscInt *);
283
284 PetscCall(PetscObjectQueryFunction((PetscObject)dm, "DMGetNumberOfCells_C", &method_DMShellGetNumberOfCells));
285 if (method_DMShellGetNumberOfCells) {
286 PetscCall(method_DMShellGetNumberOfCells(dm, &ncells));
287 } else
288 SETERRQ(PetscObjectComm((PetscObject)sw), PETSC_ERR_SUP, "Cannot determine the number of cells for the DMSHELL object. User must provide a method via PetscObjectComposeFunction( (PetscObject)shelldm, \"DMGetNumberOfCells_C\", your_function_to_compute_number_of_cells);");
289 } else SETERRQ(PetscObjectComm((PetscObject)sw), PETSC_ERR_SUP, "Cannot determine the number of cells for a DM not of type DA, PLEX or SHELL");
290
291 /* setup */
292 PetscCall(DMSwarmSortSetup(ctx, sw, ncells));
293 PetscFunctionReturn(PETSC_SUCCESS);
294 }
295
296 /*@
297 DMSwarmSortRestoreAccess - Invalidates the `DMSWARM` point sorting context previously computed with `DMSwarmSortGetAccess()`
298
299 Not Collective
300
301 Input Parameter:
302 . sw - a `DMSWARM` object
303
304 Level: advanced
305
306 Note:
307 You must call `DMSwarmSortGetAccess()` before calling `DMSwarmSortRestoreAccess()`
308
309 .seealso: `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSortGetAccess()`
310 @*/
DMSwarmSortRestoreAccess(DM sw)311 PetscErrorCode DMSwarmSortRestoreAccess(DM sw)
312 {
313 DMSwarmCellDM celldm;
314 DMSwarmSort ctx;
315
316 PetscFunctionBegin;
317 PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
318 PetscCall(DMSwarmCellDMGetSort(celldm, &ctx));
319 if (!ctx) PetscFunctionReturn(PETSC_SUCCESS);
320 PetscCheck(ctx->isvalid, PetscObjectComm((PetscObject)sw), PETSC_ERR_SUP, "You must call DMSwarmSortGetAccess() before calling DMSwarmSortRestoreAccess()");
321 ctx->isvalid = PETSC_FALSE;
322 PetscFunctionReturn(PETSC_SUCCESS);
323 }
324
325 /*@
326 DMSwarmSortGetIsValid - Gets the isvalid flag associated with a `DMSWARM` point sorting context
327
328 Not Collective
329
330 Input Parameter:
331 . sw - a `DMSWARM` object
332
333 Output Parameter:
334 . isvalid - flag indicating whether the sort context is up-to-date
335
336 Level: advanced
337
338 .seealso: `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSortGetAccess()`
339 @*/
DMSwarmSortGetIsValid(DM sw,PetscBool * isvalid)340 PetscErrorCode DMSwarmSortGetIsValid(DM sw, PetscBool *isvalid)
341 {
342 DMSwarmCellDM celldm;
343 DMSwarmSort ctx;
344
345 PetscFunctionBegin;
346 PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
347 PetscCall(DMSwarmCellDMGetSort(celldm, &ctx));
348 if (!ctx) {
349 *isvalid = PETSC_FALSE;
350 PetscFunctionReturn(PETSC_SUCCESS);
351 }
352 *isvalid = ctx->isvalid;
353 PetscFunctionReturn(PETSC_SUCCESS);
354 }
355
356 /*@
357 DMSwarmSortGetSizes - Gets the sizes associated with a `DMSWARM` point sorting context
358
359 Not Collective
360
361 Input Parameter:
362 . sw - a `DMSWARM` object
363
364 Output Parameters:
365 + ncells - number of cells within the sort context (pass `NULL` to ignore)
366 - npoints - number of points used to create the sort context (pass `NULL` to ignore)
367
368 Level: advanced
369
370 .seealso: `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSortGetAccess()`
371 @*/
DMSwarmSortGetSizes(DM sw,PetscInt * ncells,PetscInt * npoints)372 PetscErrorCode DMSwarmSortGetSizes(DM sw, PetscInt *ncells, PetscInt *npoints)
373 {
374 DMSwarmCellDM celldm;
375 DMSwarmSort ctx;
376
377 PetscFunctionBegin;
378 PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
379 PetscCall(DMSwarmCellDMGetSort(celldm, &ctx));
380 if (!ctx) {
381 if (ncells) *ncells = 0;
382 if (npoints) *npoints = 0;
383 PetscFunctionReturn(PETSC_SUCCESS);
384 }
385 if (ncells) *ncells = ctx->ncells;
386 if (npoints) *npoints = ctx->npoints;
387 PetscFunctionReturn(PETSC_SUCCESS);
388 }
389