1 #include <petscdm.h>
2 #include <petscdmda.h>
3 #include <petscdmplex.h>
4 #include <petscdmswarm.h>
5 #include <petsc/private/dmswarmimpl.h>
6
sort_CompareSwarmPoint(const void * dataA,const void * dataB)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 }
23
DMSwarmSortApplyCellIndexSort(DMSwarmSort ctx)24 PetscErrorCode DMSwarmSortApplyCellIndexSort(DMSwarmSort ctx)
25 {
26 PetscFunctionBegin;
27 qsort(ctx->list,ctx->npoints,sizeof(SwarmPoint),sort_CompareSwarmPoint);
28 PetscFunctionReturn(0);
29 }
30
DMSwarmSortCreate(DMSwarmSort * _ctx)31 PetscErrorCode DMSwarmSortCreate(DMSwarmSort *_ctx)
32 {
33 PetscErrorCode ierr;
34 DMSwarmSort ctx;
35
36 PetscFunctionBegin;
37 ierr = PetscMalloc1(1,&ctx);CHKERRQ(ierr);
38 ierr = PetscMemzero(ctx,sizeof(struct _p_DMSwarmSort));CHKERRQ(ierr);
39 ctx->isvalid = PETSC_FALSE;
40 ctx->ncells = 0;
41 ctx->npoints = 0;
42 ierr = PetscMalloc1(1,&ctx->pcell_offsets);CHKERRQ(ierr);
43 ierr = PetscMalloc1(1,&ctx->list);CHKERRQ(ierr);
44 *_ctx = ctx;
45 PetscFunctionReturn(0);
46 }
47
DMSwarmSortSetup(DMSwarmSort ctx,DM dm,PetscInt ncells)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
55 PetscFunctionBegin;
56 if (!ctx) PetscFunctionReturn(0);
57 if (ctx->isvalid) PetscFunctionReturn(0);
58
59 ierr = PetscLogEventBegin(DMSWARM_Sort,0,0,0,0);CHKERRQ(ierr);
60 /* check the number of cells */
61 if (ncells != ctx->ncells) {
62 ierr = PetscRealloc(sizeof(PetscInt)*(ncells + 1),&ctx->pcell_offsets);CHKERRQ(ierr);
63 ctx->ncells = ncells;
64 }
65 ierr = PetscArrayzero(ctx->pcell_offsets,ctx->ncells + 1);CHKERRQ(ierr);
66
67 /* get the number of points */
68 ierr = DMSwarmGetLocalSize(dm,&npoints);CHKERRQ(ierr);
69 if (npoints != ctx->npoints) {
70 ierr = PetscRealloc(sizeof(SwarmPoint)*npoints,&ctx->list);CHKERRQ(ierr);
71 ctx->npoints = npoints;
72 }
73 ierr = PetscArrayzero(ctx->list,npoints);CHKERRQ(ierr);
74
75 ierr = DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr);
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 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr);
81
82 ierr = DMSwarmSortApplyCellIndexSort(ctx);CHKERRQ(ierr);
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 ierr = PetscLogEventEnd(DMSWARM_Sort,0,0,0,0);CHKERRQ(ierr);
100 PetscFunctionReturn(0);
101 }
102
DMSwarmSortDestroy(DMSwarmSort * _ctx)103 PetscErrorCode DMSwarmSortDestroy(DMSwarmSort *_ctx)
104 {
105 DMSwarmSort ctx;
106 PetscErrorCode ierr;
107
108 PetscFunctionBegin;
109 if (!_ctx) PetscFunctionReturn(0);
110 if (!*_ctx) PetscFunctionReturn(0);
111 ctx = *_ctx;
112 if (ctx->list) {
113 ierr = PetscFree(ctx->list);CHKERRQ(ierr);
114 }
115 if (ctx->pcell_offsets) {
116 ierr = PetscFree(ctx->pcell_offsets);CHKERRQ(ierr);
117 }
118 ierr = PetscFree(ctx);CHKERRQ(ierr);
119 *_ctx = NULL;
120 PetscFunctionReturn(0);
121 }
122
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 @*/
DMSwarmSortGetNumberOfPointsPerCell(DM dm,PetscInt e,PetscInt * npoints)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
146 PetscFunctionBegin;
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 PetscFunctionReturn(0);
155 }
156
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 @*/
DMSwarmSortGetPointsPerCell(DM dm,PetscInt e,PetscInt * npoints,PetscInt ** pidlist)176 PETSC_EXTERN PetscErrorCode DMSwarmSortGetPointsPerCell(DM dm,PetscInt e,PetscInt *npoints,PetscInt **pidlist)
177 {
178 DM_Swarm *swarm = (DM_Swarm*)dm->data;
179 PetscErrorCode ierr;
180 PetscInt points_per_cell;
181 PetscInt p,pid,pid_unsorted;
182 PetscInt *plist;
183 DMSwarmSort ctx;
184
185 PetscFunctionBegin;
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 ierr = DMSwarmSortGetNumberOfPointsPerCell(dm,e,&points_per_cell);CHKERRQ(ierr);
189 ierr = PetscMalloc1(points_per_cell,&plist);CHKERRQ(ierr);
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 PetscFunctionReturn(0);
199 }
200
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
233
234 Level: advanced
235
236 .seealso: DMSwarmSetType(), DMSwarmSortRestoreAccess()
237 @*/
DMSwarmSortGetAccess(DM dm)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
246 PetscFunctionBegin;
247 if (!swarm->sort_context) {
248 ierr = DMSwarmSortCreate(&swarm->sort_context);CHKERRQ(ierr);
249 }
250
251 /* get the number of cells */
252 ierr = DMSwarmGetCellDM(dm,&celldm);CHKERRQ(ierr);
253 ierr = PetscObjectTypeCompare((PetscObject)celldm,DMDA,&isda);CHKERRQ(ierr);
254 ierr = PetscObjectTypeCompare((PetscObject)celldm,DMPLEX,&isplex);CHKERRQ(ierr);
255 ierr = PetscObjectTypeCompare((PetscObject)celldm,DMSHELL,&isshell);CHKERRQ(ierr);
256 ncells = 0;
257 if (isda) {
258 PetscInt nel,npe;
259 const PetscInt *element;
260
261 ierr = DMDAGetElements(celldm,&nel,&npe,&element);CHKERRQ(ierr);
262 ncells = nel;
263 ierr = DMDARestoreElements(celldm,&nel,&npe,&element);CHKERRQ(ierr);
264 } else if (isplex) {
265 PetscInt ps,pe;
266
267 ierr = DMPlexGetHeightStratum(celldm,0,&ps,&pe);CHKERRQ(ierr);
268 ncells = pe - ps;
269 } else if (isshell) {
270 PetscErrorCode (*method_DMShellGetNumberOfCells)(DM,PetscInt*);
271
272 ierr = PetscObjectQueryFunction((PetscObject)celldm,"DMGetNumberOfCells_C",&method_DMShellGetNumberOfCells);CHKERRQ(ierr);
273 if (method_DMShellGetNumberOfCells) {
274 ierr = method_DMShellGetNumberOfCells(celldm,&ncells);CHKERRQ(ierr);
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 ierr = DMSwarmSortSetup(swarm->sort_context,dm,ncells);CHKERRQ(ierr);
280
281 PetscFunctionReturn(0);
282 }
283
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 @*/
DMSwarmSortRestoreAccess(DM dm)299 PETSC_EXTERN PetscErrorCode DMSwarmSortRestoreAccess(DM dm)
300 {
301 DM_Swarm *swarm = (DM_Swarm*)dm->data;
302
303 PetscFunctionBegin;
304 if (!swarm->sort_context) PetscFunctionReturn(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 PetscFunctionReturn(0);
308 }
309
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
317
318 Output parameter:
319 . isvalid - flag indicating whether the sort context is up-to-date
320
321 Level: advanced
322
323 .seealso: DMSwarmSetType(), DMSwarmSortGetAccess()
324 @*/
DMSwarmSortGetIsValid(DM dm,PetscBool * isvalid)325 PETSC_EXTERN PetscErrorCode DMSwarmSortGetIsValid(DM dm,PetscBool *isvalid)
326 {
327 DM_Swarm *swarm = (DM_Swarm*)dm->data;
328
329 PetscFunctionBegin;
330 if (!swarm->sort_context) {
331 *isvalid = PETSC_FALSE;
332 PetscFunctionReturn(0);
333 }
334 *isvalid = swarm->sort_context->isvalid;
335 PetscFunctionReturn(0);
336 }
337
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 @*/
DMSwarmSortGetSizes(DM dm,PetscInt * ncells,PetscInt * npoints)354 PETSC_EXTERN PetscErrorCode DMSwarmSortGetSizes(DM dm,PetscInt *ncells,PetscInt *npoints)
355 {
356 DM_Swarm *swarm = (DM_Swarm*)dm->data;
357
358 PetscFunctionBegin;
359 if (!swarm->sort_context) {
360 if (ncells) { *ncells = 0; }
361 if (npoints) { *npoints = 0; }
362 PetscFunctionReturn(0);
363 }
364 if (ncells) { *ncells = swarm->sort_context->ncells; }
365 if (npoints) { *npoints = swarm->sort_context->npoints; }
366 PetscFunctionReturn(0);
367 }
368