Actual source code: swarmpic_sort.c
petsc-3.13.6 2020-09-29
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;
10: SwarmPoint *pointB;
11:
12: pointA = (SwarmPoint*)dataA;
13: pointB = (SwarmPoint*)dataB;
14:
15: if (pointA->cell_index < pointB->cell_index) {
16: return -1;
17: } else if (pointA->cell_index > pointB->cell_index) {
18: return 1;
19: } else {
20: return 0;
21: }
22: }
24: PetscErrorCode DMSwarmSortApplyCellIndexSort(DMSwarmSort ctx)
25: {
27: qsort(ctx->list,ctx->npoints,sizeof(SwarmPoint),sort_CompareSwarmPoint);
28: return(0);
29: }
31: PetscErrorCode DMSwarmSortCreate(DMSwarmSort *_ctx)
32: {
34: DMSwarmSort ctx;
35:
37: PetscMalloc1(1,&ctx);
38: PetscMemzero(ctx,sizeof(struct _p_DMSwarmSort));
39: ctx->isvalid = PETSC_FALSE;
40: ctx->ncells = 0;
41: ctx->npoints = 0;
42: PetscMalloc1(1,&ctx->pcell_offsets);
43: PetscMalloc1(1,&ctx->list);
44: *_ctx = ctx;
45: return(0);
46: }
48: PetscErrorCode DMSwarmSortSetup(DMSwarmSort ctx,DM dm,PetscInt ncells)
49: {
50: PetscInt *swarm_cellid;
51: PetscInt p,npoints;
52: PetscInt tmp,c,count;
53: PetscErrorCode ierr;
54:
56: if (!ctx) return(0);
57: if (ctx->isvalid) return(0);
58:
59: PetscLogEventBegin(DMSWARM_Sort,0,0,0,0);
60: /* check the number of cells */
61: if (ncells != ctx->ncells) {
62: PetscRealloc(sizeof(PetscInt)*(ncells + 1),&ctx->pcell_offsets);
63: ctx->ncells = ncells;
64: }
65: PetscArrayzero(ctx->pcell_offsets,ctx->ncells + 1);
66:
67: /* get the number of points */
68: DMSwarmGetLocalSize(dm,&npoints);
69: if (npoints != ctx->npoints) {
70: PetscRealloc(sizeof(SwarmPoint)*npoints,&ctx->list);
71: ctx->npoints = npoints;
72: }
73: PetscArrayzero(ctx->list,npoints);
74:
75: DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
76: for (p=0; p<ctx->npoints; p++) {
77: ctx->list[p].point_index = p;
78: ctx->list[p].cell_index = swarm_cellid[p];
79: }
80: DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
81:
82: DMSwarmSortApplyCellIndexSort(ctx);
83:
84: /* sum points per cell */
85: for (p=0; p<ctx->npoints; p++) {
86: ctx->pcell_offsets[ ctx->list[p].cell_index ]++;
87: }
88:
89: /* create offset list */
90: count = 0;
91: for (c=0; c<ctx->ncells; c++) {
92: tmp = ctx->pcell_offsets[c];
93: ctx->pcell_offsets[c] = count;
94: count = count + tmp;
95: }
96: ctx->pcell_offsets[c] = count;
97:
98: ctx->isvalid = PETSC_TRUE;
99: PetscLogEventEnd(DMSWARM_Sort,0,0,0,0);
100: return(0);
101: }
103: PetscErrorCode DMSwarmSortDestroy(DMSwarmSort *_ctx)
104: {
105: DMSwarmSort ctx;
106: PetscErrorCode ierr;
107:
109: if (!_ctx) return(0);
110: if (!*_ctx) return(0);
111: ctx = *_ctx;
112: if (ctx->list) {
113: PetscFree(ctx->list);
114: }
115: if (ctx->pcell_offsets) {
116: PetscFree(ctx->pcell_offsets);
117: }
118: PetscFree(ctx);
119: *_ctx = NULL;
120: return(0);
121: }
123: /*@C
124: DMSwarmSortGetNumberOfPointsPerCell - Returns the number of points in a cell
125:
126: Not collective
127:
128: Input parameters:
129: + dm - a DMSwarm objects
130: . e - the index of the cell
131: - npoints - the number of points in the cell
132:
133: Level: advanced
134:
135: Notes:
136: You must call DMSwarmSortGetAccess() before you can call DMSwarmSortGetNumberOfPointsPerCell()
137:
138: .seealso: DMSwarmSetType(), DMSwarmSortGetAccess(), DMSwarmSortGetPointsPerCell()
139: @*/
140: PetscErrorCode DMSwarmSortGetNumberOfPointsPerCell(DM dm,PetscInt e,PetscInt *npoints)
141: {
142: DM_Swarm *swarm = (DM_Swarm*)dm->data;
143: PetscInt points_per_cell;
144: DMSwarmSort ctx;
145:
147: ctx = swarm->sort_context;
148: if (!ctx) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"The DMSwarmSort context has not been created. Must call DMSwarmSortGetAccess() first");
149: if (!ctx->isvalid) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"SwarmPointSort container is not valid. Must call DMSwarmSortGetAccess() first");
150: if (e >= ctx->ncells) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Cell index (%D) is greater than max number of local cells (%D)",e,ctx->ncells);
151: if (e < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Cell index (%D) cannot be negative",e);
152: points_per_cell = ctx->pcell_offsets[e+1] - ctx->pcell_offsets[e];
153: *npoints = points_per_cell;
154: return(0);
155: }
157: /*@C
158: DMSwarmSortGetPointsPerCell - Creates an array of point indices for all points in a cell
159:
160: Not collective
161:
162: Input parameters:
163: + dm - a DMSwarm object
164: . e - the index of the cell
165: . npoints - the number of points in the cell
166: - pidlist - array of the indices indentifying all points in cell e
167:
168: Level: advanced
169:
170: Notes:
171: - You must call DMSwarmSortGetAccess() before you can call DMSwarmSortGetPointsPerCell()
172: - The array pidlist is internally created and must be free'd by the user
173:
174: .seealso: DMSwarmSetType(), DMSwarmSortGetAccess(), DMSwarmSortGetNumberOfPointsPerCell()
175: @*/
176: PETSC_EXTERN PetscErrorCode DMSwarmSortGetPointsPerCell(DM dm,PetscInt e,PetscInt *npoints,PetscInt **pidlist)
177: {
178: DM_Swarm *swarm = (DM_Swarm*)dm->data;
180: PetscInt points_per_cell;
181: PetscInt p,pid,pid_unsorted;
182: PetscInt *plist;
183: DMSwarmSort ctx;
184:
186: ctx = swarm->sort_context;
187: if (!ctx) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"The DMSwarmSort context has not been created. Must call DMSwarmSortGetAccess() first");
188: DMSwarmSortGetNumberOfPointsPerCell(dm,e,&points_per_cell);
189: PetscMalloc1(points_per_cell,&plist);
190: for (p=0; p<points_per_cell; p++) {
191: pid = ctx->pcell_offsets[e] + p;
192: pid_unsorted = ctx->list[pid].point_index;
193: plist[p] = pid_unsorted;
194: }
195: *npoints = points_per_cell;
196: *pidlist = plist;
197:
198: return(0);
199: }
201: /*@C
202: DMSwarmSortGetAccess - Setups up a DMSwarm point sort context for efficient traversal of points within a cell
203:
204: Not collective
205:
206: Input parameter:
207: . dm - a DMSwarm object
208:
209: Calling DMSwarmSortGetAccess() creates a list which enables easy identification of all points contained in a
210: given cell. This method does not explicitly sort the data within the DMSwarm based on the cell index associated
211: with a DMSwarm point.
212:
213: The sort context is valid only for the DMSwarm points defined at the time when DMSwarmSortGetAccess() was called.
214: For example, suppose the swarm contained NP points when DMSwarmSortGetAccess() was called. If the user subsequently
215: adds 10 additional points to the swarm, the sort context is still valid, but only for the first NP points.
216: The indices associated with the 10 new additional points will not be contained within the sort context.
217: This means that the user can still safely perform queries via DMSwarmSortGetPointsPerCell() and
218: DMSwarmSortGetPointsPerCell(), however the results return will be based on the first NP points.
219:
220: If any DMSwam re-sizing method is called after DMSwarmSortGetAccess() which modifies any of the first NP entries
221: in the DMSwarm, the sort context will become invalid. Currently there are no guards to prevent the user from
222: invalidating the sort context. For this reason, we highly recommend you do not use DMSwarmRemovePointAtIndex() in
223: between calls to DMSwarmSortGetAccess() and DMSwarmSortRestoreAccess().
224:
225: To facilitate safe removal of points using the sort context, we suggest a "two pass" strategy in which the
226: first pass "marks" points for removal, and the second pass actually removes the points from the DMSwarm.
227:
228: Notes:
229: - You must call DMSwarmSortGetAccess() before you can call DMSwarmSortGetPointsPerCell() or DMSwarmSortGetNumberOfPointsPerCell()
230: - The sort context may become invalid if any re-sizing methods are applied which alter the first NP points
231: within swarm at the time DMSwarmSortGetAccess() was called.
232: - You must call DMSwarmSortRestoreAccess() when you no longer need access to the sort context
234: Level: advanced
235:
236: .seealso: DMSwarmSetType(), DMSwarmSortRestoreAccess()
237: @*/
238: PETSC_EXTERN PetscErrorCode DMSwarmSortGetAccess(DM dm)
239: {
240: DM_Swarm *swarm = (DM_Swarm*)dm->data;
241: PetscErrorCode ierr;
242: PetscInt ncells;
243: DM celldm;
244: PetscBool isda,isplex,isshell;
245:
247: if (!swarm->sort_context) {
248: DMSwarmSortCreate(&swarm->sort_context);
249: }
250:
251: /* get the number of cells */
252: DMSwarmGetCellDM(dm,&celldm);
253: PetscObjectTypeCompare((PetscObject)celldm,DMDA,&isda);
254: PetscObjectTypeCompare((PetscObject)celldm,DMPLEX,&isplex);
255: PetscObjectTypeCompare((PetscObject)celldm,DMSHELL,&isshell);
256: ncells = 0;
257: if (isda) {
258: PetscInt nel,npe;
259: const PetscInt *element;
260:
261: DMDAGetElements(celldm,&nel,&npe,&element);
262: ncells = nel;
263: DMDARestoreElements(celldm,&nel,&npe,&element);
264: } else if (isplex) {
265: PetscInt ps,pe;
266:
267: DMPlexGetHeightStratum(celldm,0,&ps,&pe);
268: ncells = pe - ps;
269: } else if (isshell) {
270: PetscErrorCode (*method_DMShellGetNumberOfCells)(DM,PetscInt*);
271:
272: PetscObjectQueryFunction((PetscObject)celldm,"DMGetNumberOfCells_C",&method_DMShellGetNumberOfCells);
273: if (method_DMShellGetNumberOfCells) {
274: method_DMShellGetNumberOfCells(celldm,&ncells);
275: } 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 );");
276: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Cannot determine the number of cells for a DM not of type DA, PLEX or SHELL");
277:
278: /* setup */
279: DMSwarmSortSetup(swarm->sort_context,dm,ncells);
280:
281: return(0);
282: }
284: /*@C
285: DMSwarmSortRestoreAccess - Invalidates the DMSwarm point sorting context
286:
287: Not collective
288:
289: Input parameter:
290: . dm - a DMSwarm object
291:
292: Level: advanced
293:
294: Note:
295: You must call DMSwarmSortGetAccess() before calling DMSwarmSortRestoreAccess()
296:
297: .seealso: DMSwarmSetType(), DMSwarmSortGetAccess()
298: @*/
299: PETSC_EXTERN PetscErrorCode DMSwarmSortRestoreAccess(DM dm)
300: {
301: DM_Swarm *swarm = (DM_Swarm*)dm->data;
302:
304: if (!swarm->sort_context) return(0);
305: if (!swarm->sort_context->isvalid) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"You must call DMSwarmSortGetAccess() before calling DMSwarmSortRestoreAccess()");
306: swarm->sort_context->isvalid = PETSC_FALSE;
307: return(0);
308: }
310: /*@C
311: DMSwarmSortGetIsValid - Gets the isvalid flag associated with a DMSwarm point sorting context
312:
313: Not collective
314:
315: Input parameter:
316: . dm - a DMSwarm object
318: Output parameter:
319: . isvalid - flag indicating whether the sort context is up-to-date
321: Level: advanced
322:
323: .seealso: DMSwarmSetType(), DMSwarmSortGetAccess()
324: @*/
325: PETSC_EXTERN PetscErrorCode DMSwarmSortGetIsValid(DM dm,PetscBool *isvalid)
326: {
327: DM_Swarm *swarm = (DM_Swarm*)dm->data;
328:
330: if (!swarm->sort_context) {
331: *isvalid = PETSC_FALSE;
332: return(0);
333: }
334: *isvalid = swarm->sort_context->isvalid;
335: return(0);
336: }
338: /*@C
339: DMSwarmSortGetSizes - Gets the sizes associated with a DMSwarm point sorting context
340:
341: Not collective
342:
343: Input parameter:
344: . dm - a DMSwarm object
345:
346: Output parameters:
347: + ncells - number of cells within the sort context (pass NULL to ignore)
348: - npoints - number of points used to create the sort context (pass NULL to ignore)
349:
350: Level: advanced
351:
352: .seealso: DMSwarmSetType(), DMSwarmSortGetAccess()
353: @*/
354: PETSC_EXTERN PetscErrorCode DMSwarmSortGetSizes(DM dm,PetscInt *ncells,PetscInt *npoints)
355: {
356: DM_Swarm *swarm = (DM_Swarm*)dm->data;
357:
359: if (!swarm->sort_context) {
360: if (ncells) { *ncells = 0; }
361: if (npoints) { *npoints = 0; }
362: return(0);
363: }
364: if (ncells) { *ncells = swarm->sort_context->ncells; }
365: if (npoints) { *npoints = swarm->sort_context->npoints; }
366: return(0);
367: }