Actual source code: swarmpic_sort.c
1: #include <petscdmda.h>
2: #include <petscdmplex.h>
3: #include <petsc/private/dmswarmimpl.h>
5: static int sort_CompareSwarmPoint(const void *dataA, const void *dataB)
6: {
7: SwarmPoint *pointA = (SwarmPoint *)dataA;
8: SwarmPoint *pointB = (SwarmPoint *)dataB;
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: }
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: }
26: static PetscErrorCode DMSwarmSortCreate(DMSwarmSort *_ctx)
27: {
28: DMSwarmSort ctx;
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: }
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;
49: PetscFunctionBegin;
50: if (!ctx) PetscFunctionReturn(PETSC_SUCCESS);
51: if (ctx->isvalid) PetscFunctionReturn(PETSC_SUCCESS);
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));
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));
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));
79: /* sum points per cell */
80: for (p = 0; p < ctx->npoints; p++) ctx->pcell_offsets[ctx->list[p].cell_index]++;
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;
91: ctx->isvalid = PETSC_TRUE;
92: PetscCall(PetscLogEventEnd(DMSWARM_Sort, 0, 0, 0, 0));
93: PetscFunctionReturn(PETSC_SUCCESS);
94: }
96: PetscErrorCode DMSwarmSortDestroy(DMSwarmSort *_ctx)
97: {
98: DMSwarmSort ctx;
100: PetscFunctionBegin;
101: if (!_ctx) PetscFunctionReturn(PETSC_SUCCESS);
102: if (!*_ctx) PetscFunctionReturn(PETSC_SUCCESS);
103: ctx = *_ctx;
104: if (ctx->list) PetscCall(PetscFree(ctx->list));
105: if (ctx->pcell_offsets) PetscCall(PetscFree(ctx->pcell_offsets));
106: PetscCall(PetscFree(ctx));
107: *_ctx = NULL;
108: PetscFunctionReturn(PETSC_SUCCESS);
109: }
111: /*@
112: DMSwarmSortGetNumberOfPointsPerCell - Returns the number of points in a cell
114: Not Collective
116: Input Parameters:
117: + sw - a `DMSWARM` objects
118: - cell - the cell number in the cell `DM`
120: Output Parameter:
121: . npoints - the number of points in the cell
123: Level: advanced
125: Notes:
126: You must call `DMSwarmSortGetAccess()` before you can call `DMSwarmSortGetNumberOfPointsPerCell()`
128: .seealso: `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSortGetAccess()`, `DMSwarmSortGetPointsPerCell()`
129: @*/
130: PetscErrorCode DMSwarmSortGetNumberOfPointsPerCell(DM sw, PetscInt cell, PetscInt *npoints)
131: {
132: DMSwarmCellDM celldm;
133: DMSwarmSort ctx;
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: }
146: /*@C
147: DMSwarmSortGetPointsPerCell - Creates an array of point indices for all points in a cell
149: Not Collective
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
157: Level: advanced
159: Note:
160: You must call `DMSwarmSortGetAccess()` before you can call `DMSwarmSortGetPointsPerCell()`, and call `DMSwarmRestorePointsPerCell()` afterwards
162: .seealso: `DMSWARM`, `DMSwarmSetType()`, `DMSwarmRestorePointsPerCell()`, `DMSwarmSortGetAccess()`, `DMSwarmSortGetNumberOfPointsPerCell()`
163: @*/
164: PetscErrorCode DMSwarmSortGetPointsPerCell(DM sw, PetscInt cell, PetscInt *npoints, PetscInt **pidlist)
165: {
166: DMSwarmCellDM celldm;
167: PetscInt pid, pid_unsorted;
168: DMSwarmSort ctx;
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: }
184: /*@C
185: DMSwarmSortRestorePointsPerCell - Restores an array of point indices for all points in a cell
187: Not Collective
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
195: Level: advanced
197: Note:
198: You must call `DMSwarmSortGetAccess()` and `DMSwarmSortGetPointsPerCell()` before you can call `DMSwarmSortRestorePointsPerCell()`
200: .seealso: `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSortGetPointsPerCell()`, `DMSwarmSortGetAccess()`, `DMSwarmSortGetNumberOfPointsPerCell()`
201: @*/
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: }
209: /*@
210: DMSwarmSortGetAccess - Setups up a `DMSWARM` point sort context for efficient traversal of points within a cell
212: Not Collective
214: Input Parameter:
215: . sw - a `DMSWARM` object
217: Level: advanced
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.
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.
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()`.
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`
239: You must call `DMSwarmSortGetAccess()` before you can call `DMSwarmSortGetPointsPerCell()` or `DMSwarmSortGetNumberOfPointsPerCell()`
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.
244: You must call `DMSwarmSortRestoreAccess()` when you no longer need access to the sort context
246: .seealso: `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSortRestoreAccess()`
247: @*/
248: PetscErrorCode DMSwarmSortGetAccess(DM sw)
249: {
250: DM dm;
251: DMSwarmCellDM celldm;
252: DMSwarmSort ctx;
253: PetscInt ncells = 0;
254: PetscBool isda, isplex, isshell;
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: }
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;
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;
279: PetscCall(DMPlexGetHeightStratum(dm, 0, &ps, &pe));
280: ncells = pe - ps;
281: } else if (isshell) {
282: PetscErrorCode (*method_DMShellGetNumberOfCells)(DM, PetscInt *);
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");
291: /* setup */
292: PetscCall(DMSwarmSortSetup(ctx, sw, ncells));
293: PetscFunctionReturn(PETSC_SUCCESS);
294: }
296: /*@
297: DMSwarmSortRestoreAccess - Invalidates the `DMSWARM` point sorting context previously computed with `DMSwarmSortGetAccess()`
299: Not Collective
301: Input Parameter:
302: . sw - a `DMSWARM` object
304: Level: advanced
306: Note:
307: You must call `DMSwarmSortGetAccess()` before calling `DMSwarmSortRestoreAccess()`
309: .seealso: `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSortGetAccess()`
310: @*/
311: PetscErrorCode DMSwarmSortRestoreAccess(DM sw)
312: {
313: DMSwarmCellDM celldm;
314: DMSwarmSort ctx;
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: }
325: /*@
326: DMSwarmSortGetIsValid - Gets the isvalid flag associated with a `DMSWARM` point sorting context
328: Not Collective
330: Input Parameter:
331: . sw - a `DMSWARM` object
333: Output Parameter:
334: . isvalid - flag indicating whether the sort context is up-to-date
336: Level: advanced
338: .seealso: `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSortGetAccess()`
339: @*/
340: PetscErrorCode DMSwarmSortGetIsValid(DM sw, PetscBool *isvalid)
341: {
342: DMSwarmCellDM celldm;
343: DMSwarmSort ctx;
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: }
356: /*@
357: DMSwarmSortGetSizes - Gets the sizes associated with a `DMSWARM` point sorting context
359: Not Collective
361: Input Parameter:
362: . sw - a `DMSWARM` object
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)
368: Level: advanced
370: .seealso: `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSortGetAccess()`
371: @*/
372: PetscErrorCode DMSwarmSortGetSizes(DM sw, PetscInt *ncells, PetscInt *npoints)
373: {
374: DMSwarmCellDM celldm;
375: DMSwarmSort ctx;
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: }