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: PetscInt *swarm_cellid;
44: PetscInt p, npoints;
45: PetscInt tmp, c, count;
47: PetscFunctionBegin;
48: if (!ctx) PetscFunctionReturn(PETSC_SUCCESS);
49: if (ctx->isvalid) PetscFunctionReturn(PETSC_SUCCESS);
51: PetscCall(PetscLogEventBegin(DMSWARM_Sort, 0, 0, 0, 0));
52: /* check the number of cells */
53: if (ncells != ctx->ncells) {
54: PetscCall(PetscRealloc(sizeof(PetscInt) * (ncells + 1), &ctx->pcell_offsets));
55: ctx->ncells = ncells;
56: }
57: PetscCall(PetscArrayzero(ctx->pcell_offsets, ctx->ncells + 1));
59: /* get the number of points */
60: PetscCall(DMSwarmGetLocalSize(dm, &npoints));
61: if (npoints != ctx->npoints) {
62: PetscCall(PetscRealloc(sizeof(SwarmPoint) * npoints, &ctx->list));
63: ctx->npoints = npoints;
64: }
65: PetscCall(PetscArrayzero(ctx->list, npoints));
67: PetscCall(DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid));
68: for (p = 0; p < ctx->npoints; p++) {
69: ctx->list[p].point_index = p;
70: ctx->list[p].cell_index = swarm_cellid[p];
71: }
72: PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid));
73: PetscCall(DMSwarmSortApplyCellIndexSort(ctx));
75: /* sum points per cell */
76: for (p = 0; p < ctx->npoints; p++) ctx->pcell_offsets[ctx->list[p].cell_index]++;
78: /* create offset list */
79: count = 0;
80: for (c = 0; c < ctx->ncells; c++) {
81: tmp = ctx->pcell_offsets[c];
82: ctx->pcell_offsets[c] = count;
83: count = count + tmp;
84: }
85: ctx->pcell_offsets[c] = count;
87: ctx->isvalid = PETSC_TRUE;
88: PetscCall(PetscLogEventEnd(DMSWARM_Sort, 0, 0, 0, 0));
89: PetscFunctionReturn(PETSC_SUCCESS);
90: }
92: PetscErrorCode DMSwarmSortDestroy(DMSwarmSort *_ctx)
93: {
94: DMSwarmSort ctx;
96: PetscFunctionBegin;
97: if (!_ctx) PetscFunctionReturn(PETSC_SUCCESS);
98: if (!*_ctx) PetscFunctionReturn(PETSC_SUCCESS);
99: ctx = *_ctx;
100: if (ctx->list) PetscCall(PetscFree(ctx->list));
101: if (ctx->pcell_offsets) PetscCall(PetscFree(ctx->pcell_offsets));
102: PetscCall(PetscFree(ctx));
103: *_ctx = NULL;
104: PetscFunctionReturn(PETSC_SUCCESS);
105: }
107: /*@
108: DMSwarmSortGetNumberOfPointsPerCell - Returns the number of points in a cell
110: Not Collective
112: Input Parameters:
113: + dm - a `DMSWARM` objects
114: - e - the index of the cell
116: Output Parameter:
117: . npoints - the number of points in the cell
119: Level: advanced
121: Notes:
122: You must call `DMSwarmSortGetAccess()` before you can call `DMSwarmSortGetNumberOfPointsPerCell()`
124: .seealso: `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSortGetAccess()`, `DMSwarmSortGetPointsPerCell()`
125: @*/
126: PetscErrorCode DMSwarmSortGetNumberOfPointsPerCell(DM dm, PetscInt e, PetscInt *npoints)
127: {
128: DM_Swarm *swarm = (DM_Swarm *)dm->data;
129: PetscInt points_per_cell;
130: DMSwarmSort ctx;
132: PetscFunctionBegin;
133: ctx = swarm->sort_context;
134: PetscCheck(ctx, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "The DMSwarmSort context has not been created. Must call DMSwarmSortGetAccess() first");
135: PetscCheck(ctx->isvalid, PETSC_COMM_SELF, PETSC_ERR_USER, "SwarmPointSort container is not valid. Must call DMSwarmSortGetAccess() first");
136: PetscCheck(e < ctx->ncells, PETSC_COMM_SELF, PETSC_ERR_USER, "Cell index (%" PetscInt_FMT ") is greater than max number of local cells (%" PetscInt_FMT ")", e, ctx->ncells);
137: PetscCheck(e >= 0, PETSC_COMM_SELF, PETSC_ERR_USER, "Cell index (%" PetscInt_FMT ") cannot be negative", e);
138: points_per_cell = ctx->pcell_offsets[e + 1] - ctx->pcell_offsets[e];
139: *npoints = points_per_cell;
140: PetscFunctionReturn(PETSC_SUCCESS);
141: }
143: /*@C
144: DMSwarmSortGetPointsPerCell - Creates an array of point indices for all points in a cell
146: Not Collective
148: Input Parameters:
149: + dm - a `DMSWARM` object
150: . e - the index of the cell
151: . npoints - the number of points in the cell
152: - pidlist - array of the indices identifying all points in cell e
154: Level: advanced
156: Notes:
157: You must call `DMSwarmSortGetAccess()` before you can call `DMSwarmSortGetPointsPerCell()`
159: The array `pidlist` is internally created and must be free'd by the user
161: .seealso: `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSortGetAccess()`, `DMSwarmSortGetNumberOfPointsPerCell()`
162: @*/
163: PetscErrorCode DMSwarmSortGetPointsPerCell(DM dm, PetscInt e, PetscInt *npoints, PetscInt **pidlist)
164: {
165: DM_Swarm *swarm = (DM_Swarm *)dm->data;
166: PetscInt points_per_cell;
167: PetscInt p, pid, pid_unsorted;
168: PetscInt *plist;
169: DMSwarmSort ctx;
171: PetscFunctionBegin;
172: ctx = swarm->sort_context;
173: PetscCheck(ctx, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "The DMSwarmSort context has not been created. Must call DMSwarmSortGetAccess() first");
174: PetscCall(DMSwarmSortGetNumberOfPointsPerCell(dm, e, &points_per_cell));
175: PetscCall(PetscMalloc1(points_per_cell, &plist));
176: for (p = 0; p < points_per_cell; p++) {
177: pid = ctx->pcell_offsets[e] + p;
178: pid_unsorted = ctx->list[pid].point_index;
179: plist[p] = pid_unsorted;
180: }
181: *npoints = points_per_cell;
182: *pidlist = plist;
183: PetscFunctionReturn(PETSC_SUCCESS);
184: }
186: /*@
187: DMSwarmSortGetAccess - Setups up a `DMSWARM` point sort context for efficient traversal of points within a cell
189: Not Collective
191: Input Parameter:
192: . dm - a `DMSWARM` object
194: Level: advanced
196: Notes:
197: Calling `DMSwarmSortGetAccess()` creates a list which enables easy identification of all points contained in a
198: given cell. This method does not explicitly sort the data within the `DMSWARM` based on the cell index associated
199: with a `DMSWARM` point.
201: The sort context is valid only for the `DMSWARM` points defined at the time when `DMSwarmSortGetAccess()` was called.
202: For example, suppose the swarm contained NP points when `DMSwarmSortGetAccess()` was called. If the user subsequently
203: adds 10 additional points to the swarm, the sort context is still valid, but only for the first NP points.
204: The indices associated with the 10 new additional points will not be contained within the sort context.
205: This means that the user can still safely perform queries via `DMSwarmSortGetPointsPerCell()` and
206: `DMSwarmSortGetPointsPerCell()`, however the results return will be based on the first NP points.
208: If any` DMSWARM` re-sizing method is called after `DMSwarmSortGetAccess()` which modifies any of the first NP entries
209: in the `DMSWARM`, the sort context will become invalid. Currently there are no guards to prevent the user from
210: invalidating the sort context. For this reason, we highly recommend you do not use `DMSwarmRemovePointAtIndex()` in
211: between calls to `DMSwarmSortGetAccess()` and `DMSwarmSortRestoreAccess()`.
213: To facilitate safe removal of points using the sort context, we suggest a "two pass" strategy in which the
214: first pass "marks" points for removal, and the second pass actually removes the points from the `DMSWARM`
216: You must call `DMSwarmSortGetAccess()` before you can call `DMSwarmSortGetPointsPerCell()` or `DMSwarmSortGetNumberOfPointsPerCell()`
218: The sort context may become invalid if any re-sizing methods are applied which alter the first NP points
219: within swarm at the time `DMSwarmSortGetAccess()` was called.
221: You must call `DMSwarmSortRestoreAccess()` when you no longer need access to the sort context
223: .seealso: `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSortRestoreAccess()`
224: @*/
225: PetscErrorCode DMSwarmSortGetAccess(DM dm)
226: {
227: DM_Swarm *swarm = (DM_Swarm *)dm->data;
228: PetscInt ncells;
229: DM celldm;
230: PetscBool isda, isplex, isshell;
232: PetscFunctionBegin;
233: if (!swarm->sort_context) PetscCall(DMSwarmSortCreate(&swarm->sort_context));
235: /* get the number of cells */
236: PetscCall(DMSwarmGetCellDM(dm, &celldm));
237: PetscCall(PetscObjectTypeCompare((PetscObject)celldm, DMDA, &isda));
238: PetscCall(PetscObjectTypeCompare((PetscObject)celldm, DMPLEX, &isplex));
239: PetscCall(PetscObjectTypeCompare((PetscObject)celldm, DMSHELL, &isshell));
240: ncells = 0;
241: if (isda) {
242: PetscInt nel, npe;
243: const PetscInt *element;
245: PetscCall(DMDAGetElements(celldm, &nel, &npe, &element));
246: ncells = nel;
247: PetscCall(DMDARestoreElements(celldm, &nel, &npe, &element));
248: } else if (isplex) {
249: PetscInt ps, pe;
251: PetscCall(DMPlexGetHeightStratum(celldm, 0, &ps, &pe));
252: ncells = pe - ps;
253: } else if (isshell) {
254: PetscErrorCode (*method_DMShellGetNumberOfCells)(DM, PetscInt *);
256: PetscCall(PetscObjectQueryFunction((PetscObject)celldm, "DMGetNumberOfCells_C", &method_DMShellGetNumberOfCells));
257: if (method_DMShellGetNumberOfCells) {
258: PetscCall(method_DMShellGetNumberOfCells(celldm, &ncells));
259: } else
260: SETERRQ(PetscObjectComm((PetscObject)dm), 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);");
261: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Cannot determine the number of cells for a DM not of type DA, PLEX or SHELL");
263: /* setup */
264: PetscCall(DMSwarmSortSetup(swarm->sort_context, dm, ncells));
265: PetscFunctionReturn(PETSC_SUCCESS);
266: }
268: /*@
269: DMSwarmSortRestoreAccess - Invalidates the `DMSWARM` point sorting context previously computed with `DMSwarmSortGetAccess()`
271: Not Collective
273: Input Parameter:
274: . dm - a `DMSWARM` object
276: Level: advanced
278: Note:
279: You must call `DMSwarmSortGetAccess()` before calling `DMSwarmSortRestoreAccess()`
281: .seealso: `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSortGetAccess()`
282: @*/
283: PetscErrorCode DMSwarmSortRestoreAccess(DM dm)
284: {
285: DM_Swarm *swarm = (DM_Swarm *)dm->data;
287: PetscFunctionBegin;
288: if (!swarm->sort_context) PetscFunctionReturn(PETSC_SUCCESS);
289: PetscCheck(swarm->sort_context->isvalid, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "You must call DMSwarmSortGetAccess() before calling DMSwarmSortRestoreAccess()");
290: swarm->sort_context->isvalid = PETSC_FALSE;
291: PetscFunctionReturn(PETSC_SUCCESS);
292: }
294: /*@
295: DMSwarmSortGetIsValid - Gets the isvalid flag associated with a `DMSWARM` point sorting context
297: Not Collective
299: Input Parameter:
300: . dm - a `DMSWARM` object
302: Output Parameter:
303: . isvalid - flag indicating whether the sort context is up-to-date
305: Level: advanced
307: .seealso: `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSortGetAccess()`
308: @*/
309: PetscErrorCode DMSwarmSortGetIsValid(DM dm, PetscBool *isvalid)
310: {
311: DM_Swarm *swarm = (DM_Swarm *)dm->data;
313: PetscFunctionBegin;
314: if (!swarm->sort_context) {
315: *isvalid = PETSC_FALSE;
316: PetscFunctionReturn(PETSC_SUCCESS);
317: }
318: *isvalid = swarm->sort_context->isvalid;
319: PetscFunctionReturn(PETSC_SUCCESS);
320: }
322: /*@
323: DMSwarmSortGetSizes - Gets the sizes associated with a `DMSWARM` point sorting context
325: Not Collective
327: Input Parameter:
328: . dm - a `DMSWARM` object
330: Output Parameters:
331: + ncells - number of cells within the sort context (pass `NULL` to ignore)
332: - npoints - number of points used to create the sort context (pass `NULL` to ignore)
334: Level: advanced
336: .seealso: `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSortGetAccess()`
337: @*/
338: PetscErrorCode DMSwarmSortGetSizes(DM dm, PetscInt *ncells, PetscInt *npoints)
339: {
340: DM_Swarm *swarm = (DM_Swarm *)dm->data;
342: PetscFunctionBegin;
343: if (!swarm->sort_context) {
344: if (ncells) *ncells = 0;
345: if (npoints) *npoints = 0;
346: PetscFunctionReturn(PETSC_SUCCESS);
347: }
348: if (ncells) *ncells = swarm->sort_context->ncells;
349: if (npoints) *npoints = swarm->sort_context->npoints;
350: PetscFunctionReturn(PETSC_SUCCESS);
351: }