Actual source code: swarmpic_sort.c
1: #include <petscdmda.h>
2: #include <petscdmplex.h>
3: #include <petsc/private/dmswarmimpl.h>
5: PetscClassId DMSWARMSORT_CLASSID;
7: static int sort_CompareSwarmPoint(const void *dataA, const void *dataB)
8: {
9: SwarmPoint *pointA = (SwarmPoint *)dataA;
10: SwarmPoint *pointB = (SwarmPoint *)dataB;
12: if (pointA->cell_index < pointB->cell_index) {
13: return -1;
14: } else if (pointA->cell_index > pointB->cell_index) {
15: return 1;
16: } else {
17: return 0;
18: }
19: }
21: static PetscErrorCode DMSwarmSortApplyCellIndexSort(DMSwarmSort ctx)
22: {
23: PetscFunctionBegin;
24: if (ctx->list) qsort(ctx->list, ctx->npoints, sizeof(SwarmPoint), sort_CompareSwarmPoint);
25: PetscFunctionReturn(PETSC_SUCCESS);
26: }
28: static PetscErrorCode DMSwarmSortCreate(DMSwarmSort *sort)
29: {
30: DMSwarmSort s;
32: PetscFunctionBegin;
33: PetscCall(DMInitializePackage());
34: PetscCall(PetscHeaderCreate(s, DMSWARMSORT_CLASSID, "DMSwarmSort", "Sort context for a DMSwarm", "DM", PETSC_COMM_SELF, DMSwarmSortDestroy, DMSwarmSortView));
35: PetscCall(PetscObjectSetName((PetscObject)s, "Sort"));
36: s->isvalid = PETSC_FALSE;
37: s->ncells = 0;
38: s->npoints = 0;
39: PetscCall(PetscMalloc1(1, &s->pcell_offsets));
40: PetscCall(PetscMalloc1(1, &s->list));
41: *sort = s;
42: PetscFunctionReturn(PETSC_SUCCESS);
43: }
45: static PetscErrorCode DMSwarmSortSetup(DMSwarmSort ctx, DM dm, PetscInt ncells)
46: {
47: DMSwarmCellDM celldm;
48: PetscInt *swarm_cellid;
49: PetscInt p, npoints;
50: PetscInt tmp, c, count;
51: const char *cellid;
53: PetscFunctionBegin;
54: if (!ctx) PetscFunctionReturn(PETSC_SUCCESS);
55: if (ctx->isvalid) PetscFunctionReturn(PETSC_SUCCESS);
57: PetscCall(PetscLogEventBegin(DMSWARM_Sort, 0, 0, 0, 0));
58: /* check the number of cells */
59: if (ncells != ctx->ncells) {
60: PetscCall(PetscRealloc(sizeof(PetscInt) * (ncells + 1), &ctx->pcell_offsets));
61: ctx->ncells = ncells;
62: }
63: PetscCall(PetscArrayzero(ctx->pcell_offsets, ctx->ncells + 1));
65: /* get the number of points */
66: PetscCall(DMSwarmGetLocalSize(dm, &npoints));
67: if (npoints != ctx->npoints) {
68: PetscCall(PetscRealloc(sizeof(SwarmPoint) * npoints, &ctx->list));
69: ctx->npoints = npoints;
70: }
71: PetscCall(PetscArrayzero(ctx->list, npoints));
73: PetscCall(DMSwarmGetCellDMActive(dm, &celldm));
74: PetscCall(DMSwarmCellDMGetCellID(celldm, &cellid));
75: PetscCall(DMSwarmGetField(dm, cellid, NULL, NULL, (void **)&swarm_cellid));
76: for (p = 0; p < ctx->npoints; p++) {
77: ctx->list[p].point_index = p;
78: ctx->list[p].cell_index = swarm_cellid[p];
79: PetscCheck(ctx->list[p].cell_index >= 0 && ctx->list[p].cell_index < ctx->ncells, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cell index %" PetscInt_FMT " is not in [0, %" PetscInt_FMT ")", swarm_cellid[p], ctx->ncells);
80: }
81: PetscCall(DMSwarmRestoreField(dm, cellid, NULL, NULL, (void **)&swarm_cellid));
82: PetscCall(DMSwarmSortApplyCellIndexSort(ctx));
84: /* sum points per cell */
85: for (p = 0; p < ctx->npoints; p++) ctx->pcell_offsets[ctx->list[p].cell_index]++;
87: /* create offset list */
88: count = 0;
89: for (c = 0; c < ctx->ncells; c++) {
90: tmp = ctx->pcell_offsets[c];
91: ctx->pcell_offsets[c] = count;
92: count = count + tmp;
93: }
94: ctx->pcell_offsets[c] = count;
96: ctx->isvalid = PETSC_TRUE;
97: PetscCall(PetscLogEventEnd(DMSWARM_Sort, 0, 0, 0, 0));
98: PetscFunctionReturn(PETSC_SUCCESS);
99: }
101: /*@
102: DMSwarmSortDestroy - destroy a `DMSwarmSort`
104: Collective
106: Input Parameter:
107: . sort - address of `DMSwarmSort`
109: Level: advanced
111: .seealso: `DMSwarmSort`, `DMSwarmSortCreate()`
112: @*/
113: PetscErrorCode DMSwarmSortDestroy(DMSwarmSort *sort)
114: {
115: PetscFunctionBegin;
116: if (!sort) PetscFunctionReturn(PETSC_SUCCESS);
117: if (!*sort) PetscFunctionReturn(PETSC_SUCCESS);
119: if (--((PetscObject)*sort)->refct > 0) {
120: *sort = NULL;
121: PetscFunctionReturn(PETSC_SUCCESS);
122: }
123: PetscTryTypeMethod(*sort, destroy);
124: PetscCall(PetscFree((*sort)->list));
125: PetscCall(PetscFree((*sort)->pcell_offsets));
126: PetscCall(PetscHeaderDestroy(sort));
127: PetscFunctionReturn(PETSC_SUCCESS);
128: }
130: /*@
131: DMSwarmSortView - view a `DMSwarmSort`
133: Collective
135: Input Parameters:
136: + sort - `DMSwarmSort`
137: - viewer - viewer to display sort context, for example `PETSC_VIEWER_STDOUT_WORLD`
139: Level: advanced
141: .seealso: `DMSwarmSort`, `DMSwarmSortCreate()`, `PetscViewer`
142: @*/
143: PetscErrorCode DMSwarmSortView(DMSwarmSort sort, PetscViewer viewer)
144: {
145: PetscBool iascii;
147: PetscFunctionBegin;
149: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sort), &viewer));
151: PetscCheckSameComm(sort, 1, viewer, 2);
152: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
153: if (iascii) {
154: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)sort, viewer));
155: PetscCall(PetscViewerASCIIPushTab(viewer));
156: }
157: PetscTryTypeMethod(sort, view, viewer);
158: if (iascii) PetscCall(PetscViewerASCIIPopTab(viewer));
159: PetscFunctionReturn(PETSC_SUCCESS);
160: }
162: /*@
163: DMSwarmSortGetNumberOfPointsPerCell - Returns the number of points in a cell
165: Not Collective
167: Input Parameters:
168: + sw - a `DMSWARM` objects
169: - cell - the cell number in the cell `DM`
171: Output Parameter:
172: . npoints - the number of points in the cell
174: Level: advanced
176: Notes:
177: You must call `DMSwarmSortGetAccess()` before you can call `DMSwarmSortGetNumberOfPointsPerCell()`
179: .seealso: `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSortGetAccess()`, `DMSwarmSortGetPointsPerCell()`
180: @*/
181: PetscErrorCode DMSwarmSortGetNumberOfPointsPerCell(DM sw, PetscInt cell, PetscInt *npoints)
182: {
183: DMSwarmCellDM celldm;
184: DMSwarmSort ctx;
186: PetscFunctionBegin;
187: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
188: PetscCall(DMSwarmCellDMGetSort(celldm, &ctx));
189: PetscCheck(ctx, PetscObjectComm((PetscObject)sw), PETSC_ERR_USER, "The DMSwarmSort context has not been created. Must call DMSwarmSortGetAccess() first");
190: PetscCheck(ctx->isvalid, PETSC_COMM_SELF, PETSC_ERR_USER, "SwarmPointSort container is not valid. Must call DMSwarmSortGetAccess() first");
191: 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);
192: PetscCheck(cell >= 0, PETSC_COMM_SELF, PETSC_ERR_USER, "Cell index (%" PetscInt_FMT ") cannot be negative", cell);
193: *npoints = ctx->pcell_offsets[cell + 1] - ctx->pcell_offsets[cell];
194: PetscFunctionReturn(PETSC_SUCCESS);
195: }
197: /*@C
198: DMSwarmSortGetPointsPerCell - Creates an array of point indices for all points in a cell
200: Not Collective
202: Input Parameters:
203: + sw - a `DMSWARM` object
204: . cell - the cell number in the cell `DM`
205: . npoints - the number of points in the cell
206: - pidlist - array of the indices identifying all points in cell e
208: Level: advanced
210: Note:
211: You must call `DMSwarmSortGetAccess()` before you can call `DMSwarmSortGetPointsPerCell()`, and call `DMSwarmRestorePointsPerCell()` afterwards
213: .seealso: `DMSWARM`, `DMSwarmSetType()`, `DMSwarmRestorePointsPerCell()`, `DMSwarmSortGetAccess()`, `DMSwarmSortGetNumberOfPointsPerCell()`
214: @*/
215: PetscErrorCode DMSwarmSortGetPointsPerCell(DM sw, PetscInt cell, PetscInt *npoints, PetscInt **pidlist)
216: {
217: DMSwarmCellDM celldm;
218: PetscInt pid, pid_unsorted;
219: DMSwarmSort ctx;
221: PetscFunctionBegin;
222: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
223: PetscCall(DMSwarmCellDMGetSort(celldm, &ctx));
224: PetscCheck(ctx, PetscObjectComm((PetscObject)sw), PETSC_ERR_USER, "The DMSwarmSort context has not been created. Must call DMSwarmSortGetAccess() first");
225: PetscCall(DMSwarmSortGetNumberOfPointsPerCell(sw, cell, npoints));
226: PetscCall(DMGetWorkArray(sw, *npoints, MPIU_SCALAR, pidlist));
227: for (PetscInt p = 0; p < *npoints; ++p) {
228: pid = ctx->pcell_offsets[cell] + p;
229: pid_unsorted = ctx->list[pid].point_index;
230: (*pidlist)[p] = pid_unsorted;
231: }
232: PetscFunctionReturn(PETSC_SUCCESS);
233: }
235: /*@C
236: DMSwarmSortRestorePointsPerCell - Restores an array of point indices for all points in a cell
238: Not Collective
240: Input Parameters:
241: + dm - a `DMSWARM` object
242: . e - the index of the cell
243: . npoints - the number of points in the cell
244: - pidlist - array of the indices identifying all points in cell e
246: Level: advanced
248: Note:
249: You must call `DMSwarmSortGetAccess()` and `DMSwarmSortGetPointsPerCell()` before you can call `DMSwarmSortRestorePointsPerCell()`
251: .seealso: `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSortGetPointsPerCell()`, `DMSwarmSortGetAccess()`, `DMSwarmSortGetNumberOfPointsPerCell()`
252: @*/
253: PetscErrorCode DMSwarmSortRestorePointsPerCell(DM dm, PetscInt e, PetscInt *npoints, PetscInt **pidlist)
254: {
255: PetscFunctionBegin;
256: PetscCall(DMRestoreWorkArray(dm, *npoints, MPIU_SCALAR, pidlist));
257: PetscFunctionReturn(PETSC_SUCCESS);
258: }
260: /*@
261: DMSwarmSortGetAccess - Setups up a `DMSWARM` point sort context for efficient traversal of points within a cell
263: Not Collective
265: Input Parameter:
266: . sw - a `DMSWARM` object
268: Level: advanced
270: Notes:
271: Calling `DMSwarmSortGetAccess()` creates a list which enables easy identification of all points contained in a
272: given cell. This method does not explicitly sort the data within the `DMSWARM` based on the cell index associated
273: with a `DMSWARM` point.
275: The sort context is valid only for the `DMSWARM` points defined at the time when `DMSwarmSortGetAccess()` was called.
276: For example, suppose the swarm contained NP points when `DMSwarmSortGetAccess()` was called. If the user subsequently
277: adds 10 additional points to the swarm, the sort context is still valid, but only for the first NP points.
278: The indices associated with the 10 new additional points will not be contained within the sort context.
279: This means that the user can still safely perform queries via `DMSwarmSortGetPointsPerCell()` and
280: `DMSwarmSortGetPointsPerCell()`, however the results return will be based on the first NP points.
282: If any` DMSWARM` re-sizing method is called after `DMSwarmSortGetAccess()` which modifies any of the first NP entries
283: in the `DMSWARM`, the sort context will become invalid. Currently there are no guards to prevent the user from
284: invalidating the sort context. For this reason, we highly recommend you do not use `DMSwarmRemovePointAtIndex()` in
285: between calls to `DMSwarmSortGetAccess()` and `DMSwarmSortRestoreAccess()`.
287: To facilitate safe removal of points using the sort context, we suggest a "two pass" strategy in which the
288: first pass "marks" points for removal, and the second pass actually removes the points from the `DMSWARM`
290: You must call `DMSwarmSortGetAccess()` before you can call `DMSwarmSortGetPointsPerCell()` or `DMSwarmSortGetNumberOfPointsPerCell()`
292: The sort context may become invalid if any re-sizing methods are applied which alter the first NP points
293: within swarm at the time `DMSwarmSortGetAccess()` was called.
295: You must call `DMSwarmSortRestoreAccess()` when you no longer need access to the sort context
297: .seealso: `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSortRestoreAccess()`
298: @*/
299: PetscErrorCode DMSwarmSortGetAccess(DM sw)
300: {
301: DM dm;
302: DMSwarmCellDM celldm;
303: DMSwarmSort ctx;
304: PetscInt ncells = 0;
305: PetscBool isda, isplex, isshell;
307: PetscFunctionBegin;
308: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
309: PetscCall(DMSwarmCellDMGetSort(celldm, &ctx));
310: if (!ctx) {
311: PetscCall(DMSwarmSortCreate(&ctx));
312: PetscCall(DMSwarmCellDMSetSort(celldm, ctx));
313: PetscCall(DMSwarmSortDestroy(&ctx));
314: PetscCall(DMSwarmCellDMGetSort(celldm, &ctx));
315: }
317: /* get the number of cells */
318: PetscCall(DMSwarmGetCellDM(sw, &dm));
319: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMDA, &isda));
320: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMPLEX, &isplex));
321: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMSHELL, &isshell));
322: if (isda) {
323: const PetscInt *element;
324: PetscInt nel, npe;
326: PetscCall(DMDAGetElements(dm, &nel, &npe, &element));
327: ncells = nel;
328: PetscCall(DMDARestoreElements(dm, &nel, &npe, &element));
329: } else if (isplex) {
330: PetscInt ps, pe;
332: PetscCall(DMPlexGetHeightStratum(dm, 0, &ps, &pe));
333: ncells = pe - ps;
334: } else if (isshell) {
335: PetscErrorCode (*method_DMShellGetNumberOfCells)(DM, PetscInt *);
337: PetscCall(PetscObjectQueryFunction((PetscObject)dm, "DMGetNumberOfCells_C", &method_DMShellGetNumberOfCells));
338: if (method_DMShellGetNumberOfCells) {
339: PetscCall(method_DMShellGetNumberOfCells(dm, &ncells));
340: } else
341: 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);");
342: } else SETERRQ(PetscObjectComm((PetscObject)sw), PETSC_ERR_SUP, "Cannot determine the number of cells for a DM not of type DA, PLEX or SHELL");
344: /* setup */
345: PetscCall(DMSwarmSortSetup(ctx, sw, ncells));
346: PetscFunctionReturn(PETSC_SUCCESS);
347: }
349: /*@
350: DMSwarmSortRestoreAccess - Invalidates the `DMSWARM` point sorting context previously computed with `DMSwarmSortGetAccess()`
352: Not Collective
354: Input Parameter:
355: . sw - a `DMSWARM` object
357: Level: advanced
359: Note:
360: You must call `DMSwarmSortGetAccess()` before calling `DMSwarmSortRestoreAccess()`
362: .seealso: `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSortGetAccess()`
363: @*/
364: PetscErrorCode DMSwarmSortRestoreAccess(DM sw)
365: {
366: DMSwarmCellDM celldm;
367: DMSwarmSort ctx;
369: PetscFunctionBegin;
370: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
371: PetscCall(DMSwarmCellDMGetSort(celldm, &ctx));
372: if (!ctx) PetscFunctionReturn(PETSC_SUCCESS);
373: PetscCheck(ctx->isvalid, PetscObjectComm((PetscObject)sw), PETSC_ERR_SUP, "You must call DMSwarmSortGetAccess() before calling DMSwarmSortRestoreAccess()");
374: ctx->isvalid = PETSC_FALSE;
375: PetscFunctionReturn(PETSC_SUCCESS);
376: }
378: /*@
379: DMSwarmSortGetIsValid - Gets the isvalid flag associated with a `DMSWARM` point sorting context
381: Not Collective
383: Input Parameter:
384: . sw - a `DMSWARM` object
386: Output Parameter:
387: . isvalid - flag indicating whether the sort context is up-to-date
389: Level: advanced
391: .seealso: `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSortGetAccess()`
392: @*/
393: PetscErrorCode DMSwarmSortGetIsValid(DM sw, PetscBool *isvalid)
394: {
395: DMSwarmCellDM celldm;
396: DMSwarmSort ctx;
398: PetscFunctionBegin;
399: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
400: PetscCall(DMSwarmCellDMGetSort(celldm, &ctx));
401: if (!ctx) {
402: *isvalid = PETSC_FALSE;
403: PetscFunctionReturn(PETSC_SUCCESS);
404: }
405: *isvalid = ctx->isvalid;
406: PetscFunctionReturn(PETSC_SUCCESS);
407: }
409: /*@
410: DMSwarmSortGetSizes - Gets the sizes associated with a `DMSWARM` point sorting context
412: Not Collective
414: Input Parameter:
415: . sw - a `DMSWARM` object
417: Output Parameters:
418: + ncells - number of cells within the sort context (pass `NULL` to ignore)
419: - npoints - number of points used to create the sort context (pass `NULL` to ignore)
421: Level: advanced
423: .seealso: `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSortGetAccess()`
424: @*/
425: PetscErrorCode DMSwarmSortGetSizes(DM sw, PetscInt *ncells, PetscInt *npoints)
426: {
427: DMSwarmCellDM celldm;
428: DMSwarmSort ctx;
430: PetscFunctionBegin;
431: PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
432: PetscCall(DMSwarmCellDMGetSort(celldm, &ctx));
433: if (!ctx) {
434: if (ncells) *ncells = 0;
435: if (npoints) *npoints = 0;
436: PetscFunctionReturn(PETSC_SUCCESS);
437: }
438: if (ncells) *ncells = ctx->ncells;
439: if (npoints) *npoints = ctx->npoints;
440: PetscFunctionReturn(PETSC_SUCCESS);
441: }