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++) {
76: ctx->pcell_offsets[ ctx->list[p].cell_index ]++;
77: }
79: /* create offset list */
80: count = 0;
81: for (c=0; c<ctx->ncells; c++) {
82: tmp = ctx->pcell_offsets[c];
83: ctx->pcell_offsets[c] = count;
84: count = count + tmp;
85: }
86: ctx->pcell_offsets[c] = count;
88: ctx->isvalid = PETSC_TRUE;
89: PetscLogEventEnd(DMSWARM_Sort,0,0,0,0);
90: return 0;
91: }
93: PetscErrorCode DMSwarmSortDestroy(DMSwarmSort *_ctx)
94: {
95: DMSwarmSort ctx;
97: if (!_ctx) return 0;
98: if (!*_ctx) return 0;
99: ctx = *_ctx;
100: if (ctx->list) {
101: PetscFree(ctx->list);
102: }
103: if (ctx->pcell_offsets) {
104: PetscFree(ctx->pcell_offsets);
105: }
106: PetscFree(ctx);
107: *_ctx = NULL;
108: return 0;
109: }
111: /*@C
112: DMSwarmSortGetNumberOfPointsPerCell - Returns the number of points in a cell
114: Not collective
116: Input parameters:
117: + dm - a DMSwarm objects
118: . e - the index of the cell
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: 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: ctx = swarm->sort_context;
139: points_per_cell = ctx->pcell_offsets[e+1] - ctx->pcell_offsets[e];
140: *npoints = points_per_cell;
141: return 0;
142: }
144: /*@C
145: DMSwarmSortGetPointsPerCell - Creates an array of point indices for all points in a cell
147: Not collective
149: Input parameters:
150: + dm - a DMSwarm object
151: . e - the index of the cell
152: . npoints - the number of points in the cell
153: - pidlist - array of the indices indentifying all points in cell e
155: Level: advanced
157: Notes:
158: You must call DMSwarmSortGetAccess() before you can call DMSwarmSortGetPointsPerCell()
160: The array pidlist is internally created and must be free'd by the user
162: .seealso: DMSwarmSetType(), DMSwarmSortGetAccess(), DMSwarmSortGetNumberOfPointsPerCell()
163: @*/
164: PETSC_EXTERN PetscErrorCode DMSwarmSortGetPointsPerCell(DM dm,PetscInt e,PetscInt *npoints,PetscInt **pidlist)
165: {
166: DM_Swarm *swarm = (DM_Swarm*)dm->data;
167: PetscInt points_per_cell;
168: PetscInt p,pid,pid_unsorted;
169: PetscInt *plist;
170: DMSwarmSort ctx;
172: ctx = swarm->sort_context;
174: DMSwarmSortGetNumberOfPointsPerCell(dm,e,&points_per_cell);
175: 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;
184: return 0;
185: }
187: /*@C
188: DMSwarmSortGetAccess - Setups up a DMSwarm point sort context for efficient traversal of points within a cell
190: Not collective
192: Input parameter:
193: . dm - a DMSwarm object
195: Calling DMSwarmSortGetAccess() creates a list which enables easy identification of all points contained in a
196: given cell. This method does not explicitly sort the data within the DMSwarm based on the cell index associated
197: with a DMSwarm point.
199: The sort context is valid only for the DMSwarm points defined at the time when DMSwarmSortGetAccess() was called.
200: For example, suppose the swarm contained NP points when DMSwarmSortGetAccess() was called. If the user subsequently
201: adds 10 additional points to the swarm, the sort context is still valid, but only for the first NP points.
202: The indices associated with the 10 new additional points will not be contained within the sort context.
203: This means that the user can still safely perform queries via DMSwarmSortGetPointsPerCell() and
204: DMSwarmSortGetPointsPerCell(), however the results return will be based on the first NP points.
206: If any DMSwam re-sizing method is called after DMSwarmSortGetAccess() which modifies any of the first NP entries
207: in the DMSwarm, the sort context will become invalid. Currently there are no guards to prevent the user from
208: invalidating the sort context. For this reason, we highly recommend you do not use DMSwarmRemovePointAtIndex() in
209: between calls to DMSwarmSortGetAccess() and DMSwarmSortRestoreAccess().
211: To facilitate safe removal of points using the sort context, we suggest a "two pass" strategy in which the
212: first pass "marks" points for removal, and the second pass actually removes the points from the DMSwarm.
214: Notes:
215: You must call DMSwarmSortGetAccess() before you can call DMSwarmSortGetPointsPerCell() or DMSwarmSortGetNumberOfPointsPerCell()
217: The sort context may become invalid if any re-sizing methods are applied which alter the first NP points
218: within swarm at the time DMSwarmSortGetAccess() was called.
220: You must call DMSwarmSortRestoreAccess() when you no longer need access to the sort context
222: Level: advanced
224: .seealso: DMSwarmSetType(), DMSwarmSortRestoreAccess()
225: @*/
226: PETSC_EXTERN PetscErrorCode DMSwarmSortGetAccess(DM dm)
227: {
228: DM_Swarm *swarm = (DM_Swarm*)dm->data;
229: PetscInt ncells;
230: DM celldm;
231: PetscBool isda,isplex,isshell;
233: if (!swarm->sort_context) {
234: DMSwarmSortCreate(&swarm->sort_context);
235: }
237: /* get the number of cells */
238: DMSwarmGetCellDM(dm,&celldm);
239: PetscObjectTypeCompare((PetscObject)celldm,DMDA,&isda);
240: PetscObjectTypeCompare((PetscObject)celldm,DMPLEX,&isplex);
241: PetscObjectTypeCompare((PetscObject)celldm,DMSHELL,&isshell);
242: ncells = 0;
243: if (isda) {
244: PetscInt nel,npe;
245: const PetscInt *element;
247: DMDAGetElements(celldm,&nel,&npe,&element);
248: ncells = nel;
249: DMDARestoreElements(celldm,&nel,&npe,&element);
250: } else if (isplex) {
251: PetscInt ps,pe;
253: DMPlexGetHeightStratum(celldm,0,&ps,&pe);
254: ncells = pe - ps;
255: } else if (isshell) {
256: PetscErrorCode (*method_DMShellGetNumberOfCells)(DM,PetscInt*);
258: PetscObjectQueryFunction((PetscObject)celldm,"DMGetNumberOfCells_C",&method_DMShellGetNumberOfCells);
259: if (method_DMShellGetNumberOfCells) {
260: method_DMShellGetNumberOfCells(celldm,&ncells);
261: } else 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);");
262: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Cannot determine the number of cells for a DM not of type DA, PLEX or SHELL");
264: /* setup */
265: DMSwarmSortSetup(swarm->sort_context,dm,ncells);
266: return 0;
267: }
269: /*@C
270: DMSwarmSortRestoreAccess - Invalidates the DMSwarm point sorting context
272: Not collective
274: Input parameter:
275: . dm - a DMSwarm object
277: Level: advanced
279: Note:
280: You must call DMSwarmSortGetAccess() before calling DMSwarmSortRestoreAccess()
282: .seealso: DMSwarmSetType(), DMSwarmSortGetAccess()
283: @*/
284: PETSC_EXTERN PetscErrorCode DMSwarmSortRestoreAccess(DM dm)
285: {
286: DM_Swarm *swarm = (DM_Swarm*)dm->data;
288: if (!swarm->sort_context) return 0;
290: swarm->sort_context->isvalid = PETSC_FALSE;
291: return 0;
292: }
294: /*@C
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: DMSwarmSetType(), DMSwarmSortGetAccess()
308: @*/
309: PETSC_EXTERN PetscErrorCode DMSwarmSortGetIsValid(DM dm,PetscBool *isvalid)
310: {
311: DM_Swarm *swarm = (DM_Swarm*)dm->data;
313: if (!swarm->sort_context) {
314: *isvalid = PETSC_FALSE;
315: return 0;
316: }
317: *isvalid = swarm->sort_context->isvalid;
318: return 0;
319: }
321: /*@C
322: DMSwarmSortGetSizes - Gets the sizes associated with a DMSwarm point sorting context
324: Not collective
326: Input parameter:
327: . dm - a DMSwarm object
329: Output parameters:
330: + ncells - number of cells within the sort context (pass NULL to ignore)
331: - npoints - number of points used to create the sort context (pass NULL to ignore)
333: Level: advanced
335: .seealso: DMSwarmSetType(), DMSwarmSortGetAccess()
336: @*/
337: PETSC_EXTERN PetscErrorCode DMSwarmSortGetSizes(DM dm,PetscInt *ncells,PetscInt *npoints)
338: {
339: DM_Swarm *swarm = (DM_Swarm*)dm->data;
341: if (!swarm->sort_context) {
342: if (ncells) *ncells = 0;
343: if (npoints) *npoints = 0;
344: return 0;
345: }
346: if (ncells) *ncells = swarm->sort_context->ncells;
347: if (npoints) *npoints = swarm->sort_context->npoints;
348: return 0;
349: }