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: }