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: 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: PetscErrorCode DMSwarmSortApplyCellIndexSort(DMSwarmSort ctx)
22: {
23: qsort(ctx->list, ctx->npoints, sizeof(SwarmPoint), sort_CompareSwarmPoint);
24: return 0;
25: }
27: PetscErrorCode DMSwarmSortCreate(DMSwarmSort *_ctx)
28: {
29: DMSwarmSort ctx;
31: PetscNew(&ctx);
32: ctx->isvalid = PETSC_FALSE;
33: ctx->ncells = 0;
34: ctx->npoints = 0;
35: PetscMalloc1(1, &ctx->pcell_offsets);
36: PetscMalloc1(1, &ctx->list);
37: *_ctx = ctx;
38: return 0;
39: }
41: PetscErrorCode DMSwarmSortSetup(DMSwarmSort ctx, DM dm, PetscInt ncells)
42: {
43: PetscInt *swarm_cellid;
44: PetscInt p, npoints;
45: PetscInt tmp, c, count;
47: if (!ctx) return 0;
48: if (ctx->isvalid) return 0;
50: PetscLogEventBegin(DMSWARM_Sort, 0, 0, 0, 0);
51: /* check the number of cells */
52: if (ncells != ctx->ncells) {
53: PetscRealloc(sizeof(PetscInt) * (ncells + 1), &ctx->pcell_offsets);
54: ctx->ncells = ncells;
55: }
56: PetscArrayzero(ctx->pcell_offsets, ctx->ncells + 1);
58: /* get the number of points */
59: DMSwarmGetLocalSize(dm, &npoints);
60: if (npoints != ctx->npoints) {
61: PetscRealloc(sizeof(SwarmPoint) * npoints, &ctx->list);
62: ctx->npoints = npoints;
63: }
64: PetscArrayzero(ctx->list, npoints);
66: DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid);
67: for (p = 0; p < ctx->npoints; p++) {
68: ctx->list[p].point_index = p;
69: ctx->list[p].cell_index = swarm_cellid[p];
70: }
71: DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid);
72: DMSwarmSortApplyCellIndexSort(ctx);
74: /* sum points per cell */
75: for (p = 0; p < ctx->npoints; p++) ctx->pcell_offsets[ctx->list[p].cell_index]++;
77: /* create offset list */
78: count = 0;
79: for (c = 0; c < ctx->ncells; c++) {
80: tmp = ctx->pcell_offsets[c];
81: ctx->pcell_offsets[c] = count;
82: count = count + tmp;
83: }
84: ctx->pcell_offsets[c] = count;
86: ctx->isvalid = PETSC_TRUE;
87: PetscLogEventEnd(DMSWARM_Sort, 0, 0, 0, 0);
88: return 0;
89: }
91: PetscErrorCode DMSwarmSortDestroy(DMSwarmSort *_ctx)
92: {
93: DMSwarmSort ctx;
95: if (!_ctx) return 0;
96: if (!*_ctx) return 0;
97: ctx = *_ctx;
98: if (ctx->list) PetscFree(ctx->list);
99: if (ctx->pcell_offsets) PetscFree(ctx->pcell_offsets);
100: PetscFree(ctx);
101: *_ctx = NULL;
102: return 0;
103: }
105: /*@C
106: DMSwarmSortGetNumberOfPointsPerCell - Returns the number of points in a cell
108: Not collective
110: Input parameters:
111: + dm - a DMSwarm objects
112: . e - the index of the cell
113: - npoints - the number of points in the cell
115: Level: advanced
117: Notes:
118: You must call DMSwarmSortGetAccess() before you can call DMSwarmSortGetNumberOfPointsPerCell()
120: .seealso: `DMSwarmSetType()`, `DMSwarmSortGetAccess()`, `DMSwarmSortGetPointsPerCell()`
121: @*/
122: PetscErrorCode DMSwarmSortGetNumberOfPointsPerCell(DM dm, PetscInt e, PetscInt *npoints)
123: {
124: DM_Swarm *swarm = (DM_Swarm *)dm->data;
125: PetscInt points_per_cell;
126: DMSwarmSort ctx;
128: ctx = swarm->sort_context;
133: points_per_cell = ctx->pcell_offsets[e + 1] - ctx->pcell_offsets[e];
134: *npoints = points_per_cell;
135: return 0;
136: }
138: /*@C
139: DMSwarmSortGetPointsPerCell - Creates an array of point indices for all points in a cell
141: Not collective
143: Input parameters:
144: + dm - a DMSwarm object
145: . e - the index of the cell
146: . npoints - the number of points in the cell
147: - pidlist - array of the indices identifying all points in cell e
149: Level: advanced
151: Notes:
152: You must call DMSwarmSortGetAccess() before you can call DMSwarmSortGetPointsPerCell()
154: The array pidlist is internally created and must be free'd by the user
156: .seealso: `DMSwarmSetType()`, `DMSwarmSortGetAccess()`, `DMSwarmSortGetNumberOfPointsPerCell()`
157: @*/
158: PETSC_EXTERN PetscErrorCode DMSwarmSortGetPointsPerCell(DM dm, PetscInt e, PetscInt *npoints, PetscInt **pidlist)
159: {
160: DM_Swarm *swarm = (DM_Swarm *)dm->data;
161: PetscInt points_per_cell;
162: PetscInt p, pid, pid_unsorted;
163: PetscInt *plist;
164: DMSwarmSort ctx;
166: ctx = swarm->sort_context;
168: DMSwarmSortGetNumberOfPointsPerCell(dm, e, &points_per_cell);
169: PetscMalloc1(points_per_cell, &plist);
170: for (p = 0; p < points_per_cell; p++) {
171: pid = ctx->pcell_offsets[e] + p;
172: pid_unsorted = ctx->list[pid].point_index;
173: plist[p] = pid_unsorted;
174: }
175: *npoints = points_per_cell;
176: *pidlist = plist;
178: return 0;
179: }
181: /*@C
182: DMSwarmSortGetAccess - Setups up a DMSwarm point sort context for efficient traversal of points within a cell
184: Not collective
186: Input parameter:
187: . dm - a DMSwarm object
189: Calling DMSwarmSortGetAccess() creates a list which enables easy identification of all points contained in a
190: given cell. This method does not explicitly sort the data within the DMSwarm based on the cell index associated
191: with a DMSwarm point.
193: The sort context is valid only for the DMSwarm points defined at the time when DMSwarmSortGetAccess() was called.
194: For example, suppose the swarm contained NP points when DMSwarmSortGetAccess() was called. If the user subsequently
195: adds 10 additional points to the swarm, the sort context is still valid, but only for the first NP points.
196: The indices associated with the 10 new additional points will not be contained within the sort context.
197: This means that the user can still safely perform queries via DMSwarmSortGetPointsPerCell() and
198: DMSwarmSortGetPointsPerCell(), however the results return will be based on the first NP points.
200: If any DMSwam re-sizing method is called after DMSwarmSortGetAccess() which modifies any of the first NP entries
201: in the DMSwarm, the sort context will become invalid. Currently there are no guards to prevent the user from
202: invalidating the sort context. For this reason, we highly recommend you do not use DMSwarmRemovePointAtIndex() in
203: between calls to DMSwarmSortGetAccess() and DMSwarmSortRestoreAccess().
205: To facilitate safe removal of points using the sort context, we suggest a "two pass" strategy in which the
206: first pass "marks" points for removal, and the second pass actually removes the points from the DMSwarm.
208: Notes:
209: You must call DMSwarmSortGetAccess() before you can call DMSwarmSortGetPointsPerCell() or DMSwarmSortGetNumberOfPointsPerCell()
211: The sort context may become invalid if any re-sizing methods are applied which alter the first NP points
212: within swarm at the time DMSwarmSortGetAccess() was called.
214: You must call DMSwarmSortRestoreAccess() when you no longer need access to the sort context
216: Level: advanced
218: .seealso: `DMSwarmSetType()`, `DMSwarmSortRestoreAccess()`
219: @*/
220: PETSC_EXTERN PetscErrorCode DMSwarmSortGetAccess(DM dm)
221: {
222: DM_Swarm *swarm = (DM_Swarm *)dm->data;
223: PetscInt ncells;
224: DM celldm;
225: PetscBool isda, isplex, isshell;
227: if (!swarm->sort_context) DMSwarmSortCreate(&swarm->sort_context);
229: /* get the number of cells */
230: DMSwarmGetCellDM(dm, &celldm);
231: PetscObjectTypeCompare((PetscObject)celldm, DMDA, &isda);
232: PetscObjectTypeCompare((PetscObject)celldm, DMPLEX, &isplex);
233: PetscObjectTypeCompare((PetscObject)celldm, DMSHELL, &isshell);
234: ncells = 0;
235: if (isda) {
236: PetscInt nel, npe;
237: const PetscInt *element;
239: DMDAGetElements(celldm, &nel, &npe, &element);
240: ncells = nel;
241: DMDARestoreElements(celldm, &nel, &npe, &element);
242: } else if (isplex) {
243: PetscInt ps, pe;
245: DMPlexGetHeightStratum(celldm, 0, &ps, &pe);
246: ncells = pe - ps;
247: } else if (isshell) {
248: PetscErrorCode (*method_DMShellGetNumberOfCells)(DM, PetscInt *);
250: PetscObjectQueryFunction((PetscObject)celldm, "DMGetNumberOfCells_C", &method_DMShellGetNumberOfCells);
251: if (method_DMShellGetNumberOfCells) {
252: method_DMShellGetNumberOfCells(celldm, &ncells);
253: } else
254: 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);");
255: } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Cannot determine the number of cells for a DM not of type DA, PLEX or SHELL");
257: /* setup */
258: DMSwarmSortSetup(swarm->sort_context, dm, ncells);
259: return 0;
260: }
262: /*@C
263: DMSwarmSortRestoreAccess - Invalidates the DMSwarm point sorting context
265: Not collective
267: Input parameter:
268: . dm - a DMSwarm object
270: Level: advanced
272: Note:
273: You must call DMSwarmSortGetAccess() before calling DMSwarmSortRestoreAccess()
275: .seealso: `DMSwarmSetType()`, `DMSwarmSortGetAccess()`
276: @*/
277: PETSC_EXTERN PetscErrorCode DMSwarmSortRestoreAccess(DM dm)
278: {
279: DM_Swarm *swarm = (DM_Swarm *)dm->data;
281: if (!swarm->sort_context) return 0;
283: swarm->sort_context->isvalid = PETSC_FALSE;
284: return 0;
285: }
287: /*@C
288: DMSwarmSortGetIsValid - Gets the isvalid flag associated with a DMSwarm point sorting context
290: Not collective
292: Input parameter:
293: . dm - a DMSwarm object
295: Output parameter:
296: . isvalid - flag indicating whether the sort context is up-to-date
298: Level: advanced
300: .seealso: `DMSwarmSetType()`, `DMSwarmSortGetAccess()`
301: @*/
302: PETSC_EXTERN PetscErrorCode DMSwarmSortGetIsValid(DM dm, PetscBool *isvalid)
303: {
304: DM_Swarm *swarm = (DM_Swarm *)dm->data;
306: if (!swarm->sort_context) {
307: *isvalid = PETSC_FALSE;
308: return 0;
309: }
310: *isvalid = swarm->sort_context->isvalid;
311: return 0;
312: }
314: /*@C
315: DMSwarmSortGetSizes - Gets the sizes associated with a DMSwarm point sorting context
317: Not collective
319: Input parameter:
320: . dm - a DMSwarm object
322: Output parameters:
323: + ncells - number of cells within the sort context (pass NULL to ignore)
324: - npoints - number of points used to create the sort context (pass NULL to ignore)
326: Level: advanced
328: .seealso: `DMSwarmSetType()`, `DMSwarmSortGetAccess()`
329: @*/
330: PETSC_EXTERN PetscErrorCode DMSwarmSortGetSizes(DM dm, PetscInt *ncells, PetscInt *npoints)
331: {
332: DM_Swarm *swarm = (DM_Swarm *)dm->data;
334: if (!swarm->sort_context) {
335: if (ncells) *ncells = 0;
336: if (npoints) *npoints = 0;
337: return 0;
338: }
339: if (ncells) *ncells = swarm->sort_context->ncells;
340: if (npoints) *npoints = swarm->sort_context->npoints;
341: return 0;
342: }