Actual source code: swarmpic_sort.c
1: #include <petscdm.h>
2: #include <petscdmda.h>
3: #include <petscdmplex.h>
4: #include <petscdmswarm.h>
5: #include <petsc/private/dmswarmimpl.h>
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 *_ctx)
29: {
30: DMSwarmSort ctx;
32: PetscFunctionBegin;
33: PetscCall(PetscNew(&ctx));
34: ctx->isvalid = PETSC_FALSE;
35: ctx->ncells = 0;
36: ctx->npoints = 0;
37: PetscCall(PetscMalloc1(1, &ctx->pcell_offsets));
38: PetscCall(PetscMalloc1(1, &ctx->list));
39: *_ctx = ctx;
40: PetscFunctionReturn(PETSC_SUCCESS);
41: }
43: static PetscErrorCode DMSwarmSortSetup(DMSwarmSort ctx, DM dm, PetscInt ncells)
44: {
45: PetscInt *swarm_cellid;
46: PetscInt p, npoints;
47: PetscInt tmp, c, count;
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(DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid));
70: for (p = 0; p < ctx->npoints; p++) {
71: ctx->list[p].point_index = p;
72: ctx->list[p].cell_index = swarm_cellid[p];
73: }
74: PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid));
75: PetscCall(DMSwarmSortApplyCellIndexSort(ctx));
77: /* sum points per cell */
78: for (p = 0; p < ctx->npoints; p++) ctx->pcell_offsets[ctx->list[p].cell_index]++;
80: /* create offset list */
81: count = 0;
82: for (c = 0; c < ctx->ncells; c++) {
83: tmp = ctx->pcell_offsets[c];
84: ctx->pcell_offsets[c] = count;
85: count = count + tmp;
86: }
87: ctx->pcell_offsets[c] = count;
89: ctx->isvalid = PETSC_TRUE;
90: PetscCall(PetscLogEventEnd(DMSWARM_Sort, 0, 0, 0, 0));
91: PetscFunctionReturn(PETSC_SUCCESS);
92: }
94: PetscErrorCode DMSwarmSortDestroy(DMSwarmSort *_ctx)
95: {
96: DMSwarmSort ctx;
98: PetscFunctionBegin;
99: if (!_ctx) PetscFunctionReturn(PETSC_SUCCESS);
100: if (!*_ctx) PetscFunctionReturn(PETSC_SUCCESS);
101: ctx = *_ctx;
102: if (ctx->list) PetscCall(PetscFree(ctx->list));
103: if (ctx->pcell_offsets) PetscCall(PetscFree(ctx->pcell_offsets));
104: PetscCall(PetscFree(ctx));
105: *_ctx = NULL;
106: PetscFunctionReturn(PETSC_SUCCESS);
107: }
109: /*@C
110: DMSwarmSortGetNumberOfPointsPerCell - Returns the number of points in a cell
112: Not Collective
114: Input Parameters:
115: + dm - a `DMSWARM` objects
116: - e - the index of the cell
118: Output Parameter:
119: . npoints - the number of points in the cell
121: Level: advanced
123: Notes:
124: You must call `DMSwarmSortGetAccess()` before you can call `DMSwarmSortGetNumberOfPointsPerCell()`
126: .seealso: `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSortGetAccess()`, `DMSwarmSortGetPointsPerCell()`
127: @*/
128: PetscErrorCode DMSwarmSortGetNumberOfPointsPerCell(DM dm, PetscInt e, PetscInt *npoints)
129: {
130: DM_Swarm *swarm = (DM_Swarm *)dm->data;
131: PetscInt points_per_cell;
132: DMSwarmSort ctx;
134: PetscFunctionBegin;
135: ctx = swarm->sort_context;
136: PetscCheck(ctx, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "The DMSwarmSort context has not been created. Must call DMSwarmSortGetAccess() first");
137: PetscCheck(ctx->isvalid, PETSC_COMM_SELF, PETSC_ERR_USER, "SwarmPointSort container is not valid. Must call DMSwarmSortGetAccess() first");
138: 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);
139: PetscCheck(e >= 0, PETSC_COMM_SELF, PETSC_ERR_USER, "Cell index (%" PetscInt_FMT ") cannot be negative", e);
140: points_per_cell = ctx->pcell_offsets[e + 1] - ctx->pcell_offsets[e];
141: *npoints = points_per_cell;
142: PetscFunctionReturn(PETSC_SUCCESS);
143: }
145: /*@C
146: DMSwarmSortGetPointsPerCell - Creates an array of point indices for all points in a cell
148: Not Collective
150: Input Parameters:
151: + dm - a `DMSWARM` object
152: . e - the index of the cell
153: . npoints - the number of points in the cell
154: - pidlist - array of the indices identifying all points in cell e
156: Level: advanced
158: Notes:
159: You must call `DMSwarmSortGetAccess()` before you can call `DMSwarmSortGetPointsPerCell()`
161: The array `pidlist` is internally created and must be free'd by the user
163: .seealso: `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSortGetAccess()`, `DMSwarmSortGetNumberOfPointsPerCell()`
164: @*/
165: PETSC_EXTERN PetscErrorCode DMSwarmSortGetPointsPerCell(DM dm, PetscInt e, PetscInt *npoints, PetscInt **pidlist)
166: {
167: DM_Swarm *swarm = (DM_Swarm *)dm->data;
168: PetscInt points_per_cell;
169: PetscInt p, pid, pid_unsorted;
170: PetscInt *plist;
171: DMSwarmSort ctx;
173: PetscFunctionBegin;
174: ctx = swarm->sort_context;
175: PetscCheck(ctx, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "The DMSwarmSort context has not been created. Must call DMSwarmSortGetAccess() first");
176: PetscCall(DMSwarmSortGetNumberOfPointsPerCell(dm, e, &points_per_cell));
177: PetscCall(PetscMalloc1(points_per_cell, &plist));
178: for (p = 0; p < points_per_cell; p++) {
179: pid = ctx->pcell_offsets[e] + p;
180: pid_unsorted = ctx->list[pid].point_index;
181: plist[p] = pid_unsorted;
182: }
183: *npoints = points_per_cell;
184: *pidlist = plist;
185: PetscFunctionReturn(PETSC_SUCCESS);
186: }
188: /*@C
189: DMSwarmSortGetAccess - Setups up a `DMSWARM` point sort context for efficient traversal of points within a cell
191: Not Collective
193: Input Parameter:
194: . dm - a `DMSWARM` object
196: Level: advanced
198: Notes:
199: Calling `DMSwarmSortGetAccess()` creates a list which enables easy identification of all points contained in a
200: given cell. This method does not explicitly sort the data within the `DMSWARM` based on the cell index associated
201: with a `DMSWARM` point.
203: The sort context is valid only for the `DMSWARM` points defined at the time when `DMSwarmSortGetAccess()` was called.
204: For example, suppose the swarm contained NP points when `DMSwarmSortGetAccess()` was called. If the user subsequently
205: adds 10 additional points to the swarm, the sort context is still valid, but only for the first NP points.
206: The indices associated with the 10 new additional points will not be contained within the sort context.
207: This means that the user can still safely perform queries via `DMSwarmSortGetPointsPerCell()` and
208: `DMSwarmSortGetPointsPerCell()`, however the results return will be based on the first NP points.
210: If any` DMSWARM` re-sizing method is called after `DMSwarmSortGetAccess()` which modifies any of the first NP entries
211: in the `DMSWARM`, the sort context will become invalid. Currently there are no guards to prevent the user from
212: invalidating the sort context. For this reason, we highly recommend you do not use `DMSwarmRemovePointAtIndex()` in
213: between calls to `DMSwarmSortGetAccess()` and `DMSwarmSortRestoreAccess()`.
215: To facilitate safe removal of points using the sort context, we suggest a "two pass" strategy in which the
216: first pass "marks" points for removal, and the second pass actually removes the points from the `DMSWARM`
218: You must call `DMSwarmSortGetAccess()` before you can call `DMSwarmSortGetPointsPerCell()` or `DMSwarmSortGetNumberOfPointsPerCell()`
220: The sort context may become invalid if any re-sizing methods are applied which alter the first NP points
221: within swarm at the time `DMSwarmSortGetAccess()` was called.
223: You must call `DMSwarmSortRestoreAccess()` when you no longer need access to the sort context
225: .seealso: `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSortRestoreAccess()`
226: @*/
227: PETSC_EXTERN PetscErrorCode DMSwarmSortGetAccess(DM dm)
228: {
229: DM_Swarm *swarm = (DM_Swarm *)dm->data;
230: PetscInt ncells;
231: DM celldm;
232: PetscBool isda, isplex, isshell;
234: PetscFunctionBegin;
235: if (!swarm->sort_context) PetscCall(DMSwarmSortCreate(&swarm->sort_context));
237: /* get the number of cells */
238: PetscCall(DMSwarmGetCellDM(dm, &celldm));
239: PetscCall(PetscObjectTypeCompare((PetscObject)celldm, DMDA, &isda));
240: PetscCall(PetscObjectTypeCompare((PetscObject)celldm, DMPLEX, &isplex));
241: PetscCall(PetscObjectTypeCompare((PetscObject)celldm, DMSHELL, &isshell));
242: ncells = 0;
243: if (isda) {
244: PetscInt nel, npe;
245: const PetscInt *element;
247: PetscCall(DMDAGetElements(celldm, &nel, &npe, &element));
248: ncells = nel;
249: PetscCall(DMDARestoreElements(celldm, &nel, &npe, &element));
250: } else if (isplex) {
251: PetscInt ps, pe;
253: PetscCall(DMPlexGetHeightStratum(celldm, 0, &ps, &pe));
254: ncells = pe - ps;
255: } else if (isshell) {
256: PetscErrorCode (*method_DMShellGetNumberOfCells)(DM, PetscInt *);
258: PetscCall(PetscObjectQueryFunction((PetscObject)celldm, "DMGetNumberOfCells_C", &method_DMShellGetNumberOfCells));
259: if (method_DMShellGetNumberOfCells) {
260: PetscCall(method_DMShellGetNumberOfCells(celldm, &ncells));
261: } else
262: 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);");
263: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Cannot determine the number of cells for a DM not of type DA, PLEX or SHELL");
265: /* setup */
266: PetscCall(DMSwarmSortSetup(swarm->sort_context, dm, ncells));
267: PetscFunctionReturn(PETSC_SUCCESS);
268: }
270: /*@C
271: DMSwarmSortRestoreAccess - Invalidates the `DMSWARM` point sorting context
273: Not Collective
275: Input Parameter:
276: . dm - a `DMSWARM` object
278: Level: advanced
280: Note:
281: You must call `DMSwarmSortGetAccess()` before calling `DMSwarmSortRestoreAccess()`
283: .seealso: `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSortGetAccess()`
284: @*/
285: PETSC_EXTERN PetscErrorCode DMSwarmSortRestoreAccess(DM dm)
286: {
287: DM_Swarm *swarm = (DM_Swarm *)dm->data;
289: PetscFunctionBegin;
290: if (!swarm->sort_context) PetscFunctionReturn(PETSC_SUCCESS);
291: PetscCheck(swarm->sort_context->isvalid, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "You must call DMSwarmSortGetAccess() before calling DMSwarmSortRestoreAccess()");
292: swarm->sort_context->isvalid = PETSC_FALSE;
293: PetscFunctionReturn(PETSC_SUCCESS);
294: }
296: /*@C
297: DMSwarmSortGetIsValid - Gets the isvalid flag associated with a `DMSWARM` point sorting context
299: Not Collective
301: Input Parameter:
302: . dm - a `DMSWARM` object
304: Output Parameter:
305: . isvalid - flag indicating whether the sort context is up-to-date
307: Level: advanced
309: .seealso: `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSortGetAccess()`
310: @*/
311: PETSC_EXTERN PetscErrorCode DMSwarmSortGetIsValid(DM dm, PetscBool *isvalid)
312: {
313: DM_Swarm *swarm = (DM_Swarm *)dm->data;
315: PetscFunctionBegin;
316: if (!swarm->sort_context) {
317: *isvalid = PETSC_FALSE;
318: PetscFunctionReturn(PETSC_SUCCESS);
319: }
320: *isvalid = swarm->sort_context->isvalid;
321: PetscFunctionReturn(PETSC_SUCCESS);
322: }
324: /*@C
325: DMSwarmSortGetSizes - Gets the sizes associated with a `DMSWARM` point sorting context
327: Not Collective
329: Input Parameter:
330: . dm - a `DMSWARM` object
332: Output Parameters:
333: + ncells - number of cells within the sort context (pass `NULL` to ignore)
334: - npoints - number of points used to create the sort context (pass `NULL` to ignore)
336: Level: advanced
338: .seealso: `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSortGetAccess()`
339: @*/
340: PETSC_EXTERN PetscErrorCode DMSwarmSortGetSizes(DM dm, PetscInt *ncells, PetscInt *npoints)
341: {
342: DM_Swarm *swarm = (DM_Swarm *)dm->data;
344: PetscFunctionBegin;
345: if (!swarm->sort_context) {
346: if (ncells) *ncells = 0;
347: if (npoints) *npoints = 0;
348: PetscFunctionReturn(PETSC_SUCCESS);
349: }
350: if (ncells) *ncells = swarm->sort_context->ncells;
351: if (npoints) *npoints = swarm->sort_context->npoints;
352: PetscFunctionReturn(PETSC_SUCCESS);
353: }