Actual source code: swarmpic_sort.c

      1: #include <petscdmda.h>
      2: #include <petscdmplex.h>
      3: #include <petsc/private/dmswarmimpl.h>

      5: static int sort_CompareSwarmPoint(const void *dataA, const void *dataB)
      6: {
      7:   SwarmPoint *pointA = (SwarmPoint *)dataA;
      8:   SwarmPoint *pointB = (SwarmPoint *)dataB;

     10:   if (pointA->cell_index < pointB->cell_index) {
     11:     return -1;
     12:   } else if (pointA->cell_index > pointB->cell_index) {
     13:     return 1;
     14:   } else {
     15:     return 0;
     16:   }
     17: }

     19: static PetscErrorCode DMSwarmSortApplyCellIndexSort(DMSwarmSort ctx)
     20: {
     21:   PetscFunctionBegin;
     22:   if (ctx->list) qsort(ctx->list, ctx->npoints, sizeof(SwarmPoint), sort_CompareSwarmPoint);
     23:   PetscFunctionReturn(PETSC_SUCCESS);
     24: }

     26: static PetscErrorCode DMSwarmSortCreate(DMSwarmSort *_ctx)
     27: {
     28:   DMSwarmSort ctx;

     30:   PetscFunctionBegin;
     31:   PetscCall(PetscNew(&ctx));
     32:   ctx->isvalid = PETSC_FALSE;
     33:   ctx->ncells  = 0;
     34:   ctx->npoints = 0;
     35:   PetscCall(PetscMalloc1(1, &ctx->pcell_offsets));
     36:   PetscCall(PetscMalloc1(1, &ctx->list));
     37:   *_ctx = ctx;
     38:   PetscFunctionReturn(PETSC_SUCCESS);
     39: }

     41: static PetscErrorCode DMSwarmSortSetup(DMSwarmSort ctx, DM dm, PetscInt ncells)
     42: {
     43:   PetscInt *swarm_cellid;
     44:   PetscInt  p, npoints;
     45:   PetscInt  tmp, c, count;

     47:   PetscFunctionBegin;
     48:   if (!ctx) PetscFunctionReturn(PETSC_SUCCESS);
     49:   if (ctx->isvalid) PetscFunctionReturn(PETSC_SUCCESS);

     51:   PetscCall(PetscLogEventBegin(DMSWARM_Sort, 0, 0, 0, 0));
     52:   /* check the number of cells */
     53:   if (ncells != ctx->ncells) {
     54:     PetscCall(PetscRealloc(sizeof(PetscInt) * (ncells + 1), &ctx->pcell_offsets));
     55:     ctx->ncells = ncells;
     56:   }
     57:   PetscCall(PetscArrayzero(ctx->pcell_offsets, ctx->ncells + 1));

     59:   /* get the number of points */
     60:   PetscCall(DMSwarmGetLocalSize(dm, &npoints));
     61:   if (npoints != ctx->npoints) {
     62:     PetscCall(PetscRealloc(sizeof(SwarmPoint) * npoints, &ctx->list));
     63:     ctx->npoints = npoints;
     64:   }
     65:   PetscCall(PetscArrayzero(ctx->list, npoints));

     67:   PetscCall(DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid));
     68:   for (p = 0; p < ctx->npoints; p++) {
     69:     ctx->list[p].point_index = p;
     70:     ctx->list[p].cell_index  = swarm_cellid[p];
     71:   }
     72:   PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid));
     73:   PetscCall(DMSwarmSortApplyCellIndexSort(ctx));

     75:   /* sum points per cell */
     76:   for (p = 0; p < ctx->npoints; p++) ctx->pcell_offsets[ctx->list[p].cell_index]++;

     78:   /* create offset list */
     79:   count = 0;
     80:   for (c = 0; c < ctx->ncells; c++) {
     81:     tmp                   = ctx->pcell_offsets[c];
     82:     ctx->pcell_offsets[c] = count;
     83:     count                 = count + tmp;
     84:   }
     85:   ctx->pcell_offsets[c] = count;

     87:   ctx->isvalid = PETSC_TRUE;
     88:   PetscCall(PetscLogEventEnd(DMSWARM_Sort, 0, 0, 0, 0));
     89:   PetscFunctionReturn(PETSC_SUCCESS);
     90: }

     92: PetscErrorCode DMSwarmSortDestroy(DMSwarmSort *_ctx)
     93: {
     94:   DMSwarmSort ctx;

     96:   PetscFunctionBegin;
     97:   if (!_ctx) PetscFunctionReturn(PETSC_SUCCESS);
     98:   if (!*_ctx) PetscFunctionReturn(PETSC_SUCCESS);
     99:   ctx = *_ctx;
    100:   if (ctx->list) PetscCall(PetscFree(ctx->list));
    101:   if (ctx->pcell_offsets) PetscCall(PetscFree(ctx->pcell_offsets));
    102:   PetscCall(PetscFree(ctx));
    103:   *_ctx = NULL;
    104:   PetscFunctionReturn(PETSC_SUCCESS);
    105: }

    107: /*@
    108:   DMSwarmSortGetNumberOfPointsPerCell - Returns the number of points in a cell

    110:   Not Collective

    112:   Input Parameters:
    113: + dm - a `DMSWARM` objects
    114: - e  - the index of the cell

    116:   Output Parameter:
    117: . npoints - the number of points in the cell

    119:   Level: advanced

    121:   Notes:
    122:   You must call `DMSwarmSortGetAccess()` before you can call `DMSwarmSortGetNumberOfPointsPerCell()`

    124: .seealso: `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSortGetAccess()`, `DMSwarmSortGetPointsPerCell()`
    125: @*/
    126: PetscErrorCode DMSwarmSortGetNumberOfPointsPerCell(DM dm, PetscInt e, PetscInt *npoints)
    127: {
    128:   DM_Swarm   *swarm = (DM_Swarm *)dm->data;
    129:   PetscInt    points_per_cell;
    130:   DMSwarmSort ctx;

    132:   PetscFunctionBegin;
    133:   ctx = swarm->sort_context;
    134:   PetscCheck(ctx, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "The DMSwarmSort context has not been created. Must call DMSwarmSortGetAccess() first");
    135:   PetscCheck(ctx->isvalid, PETSC_COMM_SELF, PETSC_ERR_USER, "SwarmPointSort container is not valid. Must call DMSwarmSortGetAccess() first");
    136:   PetscCheck(e < ctx->ncells, PETSC_COMM_SELF, PETSC_ERR_USER, "Cell index (%" PetscInt_FMT ") is greater than max number of local cells (%" PetscInt_FMT ")", e, ctx->ncells);
    137:   PetscCheck(e >= 0, PETSC_COMM_SELF, PETSC_ERR_USER, "Cell index (%" PetscInt_FMT ") cannot be negative", e);
    138:   points_per_cell = ctx->pcell_offsets[e + 1] - ctx->pcell_offsets[e];
    139:   *npoints        = points_per_cell;
    140:   PetscFunctionReturn(PETSC_SUCCESS);
    141: }

    143: /*@C
    144:   DMSwarmSortGetPointsPerCell - Creates an array of point indices for all points in a cell

    146:   Not Collective

    148:   Input Parameters:
    149: + dm      - a `DMSWARM` object
    150: . e       - the index of the cell
    151: . npoints - the number of points in the cell
    152: - pidlist - array of the indices identifying all points in cell e

    154:   Level: advanced

    156:   Notes:
    157:   You must call `DMSwarmSortGetAccess()` before you can call `DMSwarmSortGetPointsPerCell()`

    159:   The array `pidlist` is internally created and must be free'd by the user

    161: .seealso: `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSortGetAccess()`, `DMSwarmSortGetNumberOfPointsPerCell()`
    162: @*/
    163: PetscErrorCode DMSwarmSortGetPointsPerCell(DM dm, PetscInt e, PetscInt *npoints, PetscInt **pidlist)
    164: {
    165:   DM_Swarm   *swarm = (DM_Swarm *)dm->data;
    166:   PetscInt    points_per_cell;
    167:   PetscInt    p, pid, pid_unsorted;
    168:   PetscInt   *plist;
    169:   DMSwarmSort ctx;

    171:   PetscFunctionBegin;
    172:   ctx = swarm->sort_context;
    173:   PetscCheck(ctx, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "The DMSwarmSort context has not been created. Must call DMSwarmSortGetAccess() first");
    174:   PetscCall(DMSwarmSortGetNumberOfPointsPerCell(dm, e, &points_per_cell));
    175:   PetscCall(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;
    183:   PetscFunctionReturn(PETSC_SUCCESS);
    184: }

    186: /*@
    187:   DMSwarmSortGetAccess - Setups up a `DMSWARM` point sort context for efficient traversal of points within a cell

    189:   Not Collective

    191:   Input Parameter:
    192: . dm - a `DMSWARM` object

    194:   Level: advanced

    196:   Notes:
    197:   Calling `DMSwarmSortGetAccess()` creates a list which enables easy identification of all points contained in a
    198:   given cell. This method does not explicitly sort the data within the `DMSWARM` based on the cell index associated
    199:   with a `DMSWARM` point.

    201:   The sort context is valid only for the `DMSWARM` points defined at the time when `DMSwarmSortGetAccess()` was called.
    202:   For example, suppose the swarm contained NP points when `DMSwarmSortGetAccess()` was called. If the user subsequently
    203:   adds 10 additional points to the swarm, the sort context is still valid, but only for the first NP points.
    204:   The indices associated with the 10 new additional points will not be contained within the sort context.
    205:   This means that the user can still safely perform queries via `DMSwarmSortGetPointsPerCell()` and
    206:   `DMSwarmSortGetPointsPerCell()`, however the results return will be based on the first NP points.

    208:   If any` DMSWARM` re-sizing method is called after `DMSwarmSortGetAccess()` which modifies any of the first NP entries
    209:   in the `DMSWARM`, the sort context will become invalid. Currently there are no guards to prevent the user from
    210:   invalidating the sort context. For this reason, we highly recommend you do not use `DMSwarmRemovePointAtIndex()` in
    211:   between calls to `DMSwarmSortGetAccess()` and `DMSwarmSortRestoreAccess()`.

    213:   To facilitate safe removal of points using the sort context, we suggest a "two pass" strategy in which the
    214:   first pass "marks" points for removal, and the second pass actually removes the points from the `DMSWARM`

    216:   You must call `DMSwarmSortGetAccess()` before you can call `DMSwarmSortGetPointsPerCell()` or `DMSwarmSortGetNumberOfPointsPerCell()`

    218:   The sort context may become invalid if any re-sizing methods are applied which alter the first NP points
    219:   within swarm at the time `DMSwarmSortGetAccess()` was called.

    221:   You must call `DMSwarmSortRestoreAccess()` when you no longer need access to the sort context

    223: .seealso: `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSortRestoreAccess()`
    224: @*/
    225: PetscErrorCode DMSwarmSortGetAccess(DM dm)
    226: {
    227:   DM_Swarm *swarm = (DM_Swarm *)dm->data;
    228:   PetscInt  ncells;
    229:   DM        celldm;
    230:   PetscBool isda, isplex, isshell;

    232:   PetscFunctionBegin;
    233:   if (!swarm->sort_context) PetscCall(DMSwarmSortCreate(&swarm->sort_context));

    235:   /* get the number of cells */
    236:   PetscCall(DMSwarmGetCellDM(dm, &celldm));
    237:   PetscCall(PetscObjectTypeCompare((PetscObject)celldm, DMDA, &isda));
    238:   PetscCall(PetscObjectTypeCompare((PetscObject)celldm, DMPLEX, &isplex));
    239:   PetscCall(PetscObjectTypeCompare((PetscObject)celldm, DMSHELL, &isshell));
    240:   ncells = 0;
    241:   if (isda) {
    242:     PetscInt        nel, npe;
    243:     const PetscInt *element;

    245:     PetscCall(DMDAGetElements(celldm, &nel, &npe, &element));
    246:     ncells = nel;
    247:     PetscCall(DMDARestoreElements(celldm, &nel, &npe, &element));
    248:   } else if (isplex) {
    249:     PetscInt ps, pe;

    251:     PetscCall(DMPlexGetHeightStratum(celldm, 0, &ps, &pe));
    252:     ncells = pe - ps;
    253:   } else if (isshell) {
    254:     PetscErrorCode (*method_DMShellGetNumberOfCells)(DM, PetscInt *);

    256:     PetscCall(PetscObjectQueryFunction((PetscObject)celldm, "DMGetNumberOfCells_C", &method_DMShellGetNumberOfCells));
    257:     if (method_DMShellGetNumberOfCells) {
    258:       PetscCall(method_DMShellGetNumberOfCells(celldm, &ncells));
    259:     } else
    260:       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);");
    261:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Cannot determine the number of cells for a DM not of type DA, PLEX or SHELL");

    263:   /* setup */
    264:   PetscCall(DMSwarmSortSetup(swarm->sort_context, dm, ncells));
    265:   PetscFunctionReturn(PETSC_SUCCESS);
    266: }

    268: /*@
    269:   DMSwarmSortRestoreAccess - Invalidates the `DMSWARM` point sorting context previously computed with `DMSwarmSortGetAccess()`

    271:   Not Collective

    273:   Input Parameter:
    274: . dm - a `DMSWARM` object

    276:   Level: advanced

    278:   Note:
    279:   You must call `DMSwarmSortGetAccess()` before calling `DMSwarmSortRestoreAccess()`

    281: .seealso: `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSortGetAccess()`
    282: @*/
    283: PetscErrorCode DMSwarmSortRestoreAccess(DM dm)
    284: {
    285:   DM_Swarm *swarm = (DM_Swarm *)dm->data;

    287:   PetscFunctionBegin;
    288:   if (!swarm->sort_context) PetscFunctionReturn(PETSC_SUCCESS);
    289:   PetscCheck(swarm->sort_context->isvalid, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "You must call DMSwarmSortGetAccess() before calling DMSwarmSortRestoreAccess()");
    290:   swarm->sort_context->isvalid = PETSC_FALSE;
    291:   PetscFunctionReturn(PETSC_SUCCESS);
    292: }

    294: /*@
    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: `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSortGetAccess()`
    308: @*/
    309: PetscErrorCode DMSwarmSortGetIsValid(DM dm, PetscBool *isvalid)
    310: {
    311:   DM_Swarm *swarm = (DM_Swarm *)dm->data;

    313:   PetscFunctionBegin;
    314:   if (!swarm->sort_context) {
    315:     *isvalid = PETSC_FALSE;
    316:     PetscFunctionReturn(PETSC_SUCCESS);
    317:   }
    318:   *isvalid = swarm->sort_context->isvalid;
    319:   PetscFunctionReturn(PETSC_SUCCESS);
    320: }

    322: /*@
    323:   DMSwarmSortGetSizes - Gets the sizes associated with a `DMSWARM` point sorting context

    325:   Not Collective

    327:   Input Parameter:
    328: . dm - a `DMSWARM` object

    330:   Output Parameters:
    331: + ncells  - number of cells within the sort context (pass `NULL` to ignore)
    332: - npoints - number of points used to create the sort context (pass `NULL` to ignore)

    334:   Level: advanced

    336: .seealso: `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSortGetAccess()`
    337: @*/
    338: PetscErrorCode DMSwarmSortGetSizes(DM dm, PetscInt *ncells, PetscInt *npoints)
    339: {
    340:   DM_Swarm *swarm = (DM_Swarm *)dm->data;

    342:   PetscFunctionBegin;
    343:   if (!swarm->sort_context) {
    344:     if (ncells) *ncells = 0;
    345:     if (npoints) *npoints = 0;
    346:     PetscFunctionReturn(PETSC_SUCCESS);
    347:   }
    348:   if (ncells) *ncells = swarm->sort_context->ncells;
    349:   if (npoints) *npoints = swarm->sort_context->npoints;
    350:   PetscFunctionReturn(PETSC_SUCCESS);
    351: }