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;
     10:   SwarmPoint *pointB;

     12:   pointA = (SwarmPoint*)dataA;
     13:   pointB = (SwarmPoint*)dataB;

     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;

     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;

     56:   if (!ctx) return(0);
     57:   if (ctx->isvalid) return(0);

     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);

     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);

     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);

     82:   DMSwarmSortApplyCellIndexSort(ctx);

     84:   /* sum points per cell */
     85:   for (p=0; p<ctx->npoints; p++) {
     86:     ctx->pcell_offsets[ ctx->list[p].cell_index ]++;
     87:   }

     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;

     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;

    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

    126:    Not collective

    128:    Input parameters:
    129: +  dm - a DMSwarm objects
    130: .  e - the index of the cell
    131: -  npoints - the number of points in the cell

    133:    Level: advanced

    135:    Notes:
    136:    You must call DMSwarmSortGetAccess() before you can call DMSwarmSortGetNumberOfPointsPerCell()

    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;

    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

    160:    Not collective

    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

    168:    Level: advanced

    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

    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;

    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;

    198:   return(0);
    199: }

    201: /*@C
    202:    DMSwarmSortGetAccess - Setups up a DMSwarm point sort context for efficient traversal of points within a cell

    204:    Not collective

    206:    Input parameter:
    207: .  dm - a DMSwarm object

    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.

    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.

    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().

    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.

    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

    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;

    247:   if (!swarm->sort_context) {
    248:     DMSwarmSortCreate(&swarm->sort_context);
    249:   }

    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;

    261:     DMDAGetElements(celldm,&nel,&npe,&element);
    262:     ncells = nel;
    263:     DMDARestoreElements(celldm,&nel,&npe,&element);
    264:   } else if (isplex) {
    265:     PetscInt ps,pe;

    267:     DMPlexGetHeightStratum(celldm,0,&ps,&pe);
    268:     ncells = pe - ps;
    269:   } else if (isshell) {
    270:     PetscErrorCode (*method_DMShellGetNumberOfCells)(DM,PetscInt*);

    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");

    278:   /* setup */
    279:   DMSwarmSortSetup(swarm->sort_context,dm,ncells);

    281:   return(0);
    282: }

    284: /*@C
    285:    DMSwarmSortRestoreAccess - Invalidates the DMSwarm point sorting context

    287:    Not collective

    289:    Input parameter:
    290: .  dm - a DMSwarm object

    292:    Level: advanced

    294:    Note:
    295:    You must call DMSwarmSortGetAccess() before calling DMSwarmSortRestoreAccess()

    297: .seealso: DMSwarmSetType(), DMSwarmSortGetAccess()
    298: @*/
    299: PETSC_EXTERN PetscErrorCode DMSwarmSortRestoreAccess(DM dm)
    300: {
    301:   DM_Swarm *swarm = (DM_Swarm*)dm->data;

    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

    313:    Not collective

    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

    323: .seealso: DMSwarmSetType(), DMSwarmSortGetAccess()
    324: @*/
    325: PETSC_EXTERN PetscErrorCode DMSwarmSortGetIsValid(DM dm,PetscBool *isvalid)
    326: {
    327:   DM_Swarm *swarm = (DM_Swarm*)dm->data;

    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

    341:    Not collective

    343:    Input parameter:
    344: .  dm - a DMSwarm object

    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)

    350:    Level: advanced

    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;

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