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:   DMSwarmCellDM celldm;
     44:   PetscInt     *swarm_cellid;
     45:   PetscInt      p, npoints;
     46:   PetscInt      tmp, c, count;
     47:   const char   *cellid;

     49:   PetscFunctionBegin;
     50:   if (!ctx) PetscFunctionReturn(PETSC_SUCCESS);
     51:   if (ctx->isvalid) PetscFunctionReturn(PETSC_SUCCESS);

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

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

     69:   PetscCall(DMSwarmGetCellDMActive(dm, &celldm));
     70:   PetscCall(DMSwarmCellDMGetCellID(celldm, &cellid));
     71:   PetscCall(DMSwarmGetField(dm, cellid, NULL, NULL, (void **)&swarm_cellid));
     72:   for (p = 0; p < ctx->npoints; p++) {
     73:     ctx->list[p].point_index = p;
     74:     ctx->list[p].cell_index  = swarm_cellid[p];
     75:   }
     76:   PetscCall(DMSwarmRestoreField(dm, cellid, NULL, NULL, (void **)&swarm_cellid));
     77:   PetscCall(DMSwarmSortApplyCellIndexSort(ctx));

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

     82:   /* create offset list */
     83:   count = 0;
     84:   for (c = 0; c < ctx->ncells; c++) {
     85:     tmp                   = ctx->pcell_offsets[c];
     86:     ctx->pcell_offsets[c] = count;
     87:     count                 = count + tmp;
     88:   }
     89:   ctx->pcell_offsets[c] = count;

     91:   ctx->isvalid = PETSC_TRUE;
     92:   PetscCall(PetscLogEventEnd(DMSWARM_Sort, 0, 0, 0, 0));
     93:   PetscFunctionReturn(PETSC_SUCCESS);
     94: }

     96: PetscErrorCode DMSwarmSortDestroy(DMSwarmSort *ctx)
     97: {
     98:   DMSwarmSort ictx;

    100:   PetscFunctionBegin;
    101:   if (!ctx) PetscFunctionReturn(PETSC_SUCCESS);
    102:   if (!*ctx) PetscFunctionReturn(PETSC_SUCCESS);
    103:   ictx = *ctx;
    104:   PetscCall(PetscFree(ictx->list));
    105:   PetscCall(PetscFree(ictx->pcell_offsets));
    106:   PetscCall(PetscFree(ictx));
    107:   *ctx = NULL;
    108:   PetscFunctionReturn(PETSC_SUCCESS);
    109: }

    111: /*@
    112:   DMSwarmSortGetNumberOfPointsPerCell - Returns the number of points in a cell

    114:   Not Collective

    116:   Input Parameters:
    117: + sw   - a `DMSWARM` objects
    118: - cell - the cell number in the cell `DM`

    120:   Output Parameter:
    121: . npoints - the number of points in the cell

    123:   Level: advanced

    125:   Notes:
    126:   You must call `DMSwarmSortGetAccess()` before you can call `DMSwarmSortGetNumberOfPointsPerCell()`

    128: .seealso: `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSortGetAccess()`, `DMSwarmSortGetPointsPerCell()`
    129: @*/
    130: PetscErrorCode DMSwarmSortGetNumberOfPointsPerCell(DM sw, PetscInt cell, PetscInt *npoints)
    131: {
    132:   DMSwarmCellDM celldm;
    133:   DMSwarmSort   ctx;

    135:   PetscFunctionBegin;
    136:   PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
    137:   PetscCall(DMSwarmCellDMGetSort(celldm, &ctx));
    138:   PetscCheck(ctx, PetscObjectComm((PetscObject)sw), PETSC_ERR_USER, "The DMSwarmSort context has not been created. Must call DMSwarmSortGetAccess() first");
    139:   PetscCheck(ctx->isvalid, PETSC_COMM_SELF, PETSC_ERR_USER, "SwarmPointSort container is not valid. Must call DMSwarmSortGetAccess() first");
    140:   PetscCheck(cell < ctx->ncells, PETSC_COMM_SELF, PETSC_ERR_USER, "Cell index (%" PetscInt_FMT ") is greater than max number of local cells (%" PetscInt_FMT ")", cell, ctx->ncells);
    141:   PetscCheck(cell >= 0, PETSC_COMM_SELF, PETSC_ERR_USER, "Cell index (%" PetscInt_FMT ") cannot be negative", cell);
    142:   *npoints = ctx->pcell_offsets[cell + 1] - ctx->pcell_offsets[cell];
    143:   PetscFunctionReturn(PETSC_SUCCESS);
    144: }

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

    149:   Not Collective

    151:   Input Parameters:
    152: + sw      - a `DMSWARM` object
    153: . cell    - the cell number in the cell `DM`
    154: . npoints - the number of points in the cell
    155: - pidlist - array of the indices identifying all points in cell e

    157:   Level: advanced

    159:   Note:
    160:   You must call `DMSwarmSortGetAccess()` before you can call `DMSwarmSortGetPointsPerCell()`, and call `DMSwarmRestorePointsPerCell()` afterwards

    162: .seealso: `DMSWARM`, `DMSwarmSetType()`, `DMSwarmRestorePointsPerCell()`, `DMSwarmSortGetAccess()`, `DMSwarmSortGetNumberOfPointsPerCell()`
    163: @*/
    164: PetscErrorCode DMSwarmSortGetPointsPerCell(DM sw, PetscInt cell, PetscInt *npoints, PetscInt **pidlist)
    165: {
    166:   DMSwarmCellDM celldm;
    167:   PetscInt      pid, pid_unsorted;
    168:   DMSwarmSort   ctx;

    170:   PetscFunctionBegin;
    171:   PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
    172:   PetscCall(DMSwarmCellDMGetSort(celldm, &ctx));
    173:   PetscCheck(ctx, PetscObjectComm((PetscObject)sw), PETSC_ERR_USER, "The DMSwarmSort context has not been created. Must call DMSwarmSortGetAccess() first");
    174:   PetscCall(DMSwarmSortGetNumberOfPointsPerCell(sw, cell, npoints));
    175:   PetscCall(DMGetWorkArray(sw, *npoints, MPIU_SCALAR, pidlist));
    176:   for (PetscInt p = 0; p < *npoints; ++p) {
    177:     pid           = ctx->pcell_offsets[cell] + p;
    178:     pid_unsorted  = ctx->list[pid].point_index;
    179:     (*pidlist)[p] = pid_unsorted;
    180:   }
    181:   PetscFunctionReturn(PETSC_SUCCESS);
    182: }

    184: /*@C
    185:   DMSwarmSortRestorePointsPerCell - Restores an array of point indices for all points in a cell

    187:   Not Collective

    189:   Input Parameters:
    190: + dm      - a `DMSWARM` object
    191: . e       - the index of the cell
    192: . npoints - the number of points in the cell
    193: - pidlist - array of the indices identifying all points in cell e

    195:   Level: advanced

    197:   Note:
    198:   You must call `DMSwarmSortGetAccess()` and `DMSwarmSortGetPointsPerCell()` before you can call `DMSwarmSortRestorePointsPerCell()`

    200: .seealso: `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSortGetPointsPerCell()`, `DMSwarmSortGetAccess()`, `DMSwarmSortGetNumberOfPointsPerCell()`
    201: @*/
    202: PetscErrorCode DMSwarmSortRestorePointsPerCell(DM dm, PetscInt e, PetscInt *npoints, PetscInt **pidlist)
    203: {
    204:   PetscFunctionBegin;
    205:   PetscCall(DMRestoreWorkArray(dm, *npoints, MPIU_SCALAR, pidlist));
    206:   PetscFunctionReturn(PETSC_SUCCESS);
    207: }

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

    212:   Not Collective

    214:   Input Parameter:
    215: . sw - a `DMSWARM` object

    217:   Level: advanced

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

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

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

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

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

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

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

    246: .seealso: `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSortRestoreAccess()`
    247: @*/
    248: PetscErrorCode DMSwarmSortGetAccess(DM sw)
    249: {
    250:   DM            dm;
    251:   DMSwarmCellDM celldm;
    252:   DMSwarmSort   ctx;
    253:   PetscInt      ncells = 0;
    254:   PetscBool     isda, isplex, isshell;

    256:   PetscFunctionBegin;
    257:   PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
    258:   PetscCall(DMSwarmCellDMGetSort(celldm, &ctx));
    259:   if (!ctx) {
    260:     PetscCall(DMSwarmSortCreate(&ctx));
    261:     PetscCall(DMSwarmCellDMSetSort(celldm, ctx));
    262:   }

    264:   /* get the number of cells */
    265:   PetscCall(DMSwarmGetCellDM(sw, &dm));
    266:   PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMDA, &isda));
    267:   PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMPLEX, &isplex));
    268:   PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMSHELL, &isshell));
    269:   if (isda) {
    270:     const PetscInt *element;
    271:     PetscInt        nel, npe;

    273:     PetscCall(DMDAGetElements(dm, &nel, &npe, &element));
    274:     ncells = nel;
    275:     PetscCall(DMDARestoreElements(dm, &nel, &npe, &element));
    276:   } else if (isplex) {
    277:     PetscInt ps, pe;

    279:     PetscCall(DMPlexGetHeightStratum(dm, 0, &ps, &pe));
    280:     ncells = pe - ps;
    281:   } else if (isshell) {
    282:     PetscErrorCode (*method_DMShellGetNumberOfCells)(DM, PetscInt *);

    284:     PetscCall(PetscObjectQueryFunction((PetscObject)dm, "DMGetNumberOfCells_C", &method_DMShellGetNumberOfCells));
    285:     if (method_DMShellGetNumberOfCells) {
    286:       PetscCall(method_DMShellGetNumberOfCells(dm, &ncells));
    287:     } else
    288:       SETERRQ(PetscObjectComm((PetscObject)sw), 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);");
    289:   } else SETERRQ(PetscObjectComm((PetscObject)sw), PETSC_ERR_SUP, "Cannot determine the number of cells for a DM not of type DA, PLEX or SHELL");

    291:   /* setup */
    292:   PetscCall(DMSwarmSortSetup(ctx, sw, ncells));
    293:   PetscFunctionReturn(PETSC_SUCCESS);
    294: }

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

    299:   Not Collective

    301:   Input Parameter:
    302: . sw - a `DMSWARM` object

    304:   Level: advanced

    306:   Note:
    307:   You must call `DMSwarmSortGetAccess()` before calling `DMSwarmSortRestoreAccess()`

    309: .seealso: `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSortGetAccess()`
    310: @*/
    311: PetscErrorCode DMSwarmSortRestoreAccess(DM sw)
    312: {
    313:   DMSwarmCellDM celldm;
    314:   DMSwarmSort   ctx;

    316:   PetscFunctionBegin;
    317:   PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
    318:   PetscCall(DMSwarmCellDMGetSort(celldm, &ctx));
    319:   if (!ctx) PetscFunctionReturn(PETSC_SUCCESS);
    320:   PetscCheck(ctx->isvalid, PetscObjectComm((PetscObject)sw), PETSC_ERR_SUP, "You must call DMSwarmSortGetAccess() before calling DMSwarmSortRestoreAccess()");
    321:   ctx->isvalid = PETSC_FALSE;
    322:   PetscFunctionReturn(PETSC_SUCCESS);
    323: }

    325: /*@
    326:   DMSwarmSortGetIsValid - Gets the isvalid flag associated with a `DMSWARM` point sorting context

    328:   Not Collective

    330:   Input Parameter:
    331: . sw - a `DMSWARM` object

    333:   Output Parameter:
    334: . isvalid - flag indicating whether the sort context is up-to-date

    336:   Level: advanced

    338: .seealso: `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSortGetAccess()`
    339: @*/
    340: PetscErrorCode DMSwarmSortGetIsValid(DM sw, PetscBool *isvalid)
    341: {
    342:   DMSwarmCellDM celldm;
    343:   DMSwarmSort   ctx;

    345:   PetscFunctionBegin;
    346:   PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
    347:   PetscCall(DMSwarmCellDMGetSort(celldm, &ctx));
    348:   if (!ctx) {
    349:     *isvalid = PETSC_FALSE;
    350:     PetscFunctionReturn(PETSC_SUCCESS);
    351:   }
    352:   *isvalid = ctx->isvalid;
    353:   PetscFunctionReturn(PETSC_SUCCESS);
    354: }

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

    359:   Not Collective

    361:   Input Parameter:
    362: . sw - a `DMSWARM` object

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

    368:   Level: advanced

    370: .seealso: `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSortGetAccess()`
    371: @*/
    372: PetscErrorCode DMSwarmSortGetSizes(DM sw, PetscInt *ncells, PetscInt *npoints)
    373: {
    374:   DMSwarmCellDM celldm;
    375:   DMSwarmSort   ctx;

    377:   PetscFunctionBegin;
    378:   PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
    379:   PetscCall(DMSwarmCellDMGetSort(celldm, &ctx));
    380:   if (!ctx) {
    381:     if (ncells) *ncells = 0;
    382:     if (npoints) *npoints = 0;
    383:     PetscFunctionReturn(PETSC_SUCCESS);
    384:   }
    385:   if (ncells) *ncells = ctx->ncells;
    386:   if (npoints) *npoints = ctx->npoints;
    387:   PetscFunctionReturn(PETSC_SUCCESS);
    388: }