Actual source code: swarmpic.c

      1: #define PETSCDM_DLL
      2: #include <petsc/private/dmswarmimpl.h>
      3: #include <petscsf.h>
      4: #include <petscdmda.h>
      5: #include <petscdmplex.h>
      6: #include "../src/dm/impls/swarm/data_bucket.h"

      8: /*
      9:  Error chceking macto to ensure the swarm type is correct and that a cell DM has been set
     10: */
     11: #define DMSWARMPICVALID(dm) \
     12: { \
     13:   DM_Swarm *_swarm = (DM_Swarm*)(dm)->data; \
     14:   if (_swarm->swarm_type != DMSWARM_PIC) SETERRQ(PetscObjectComm((PetscObject)(dm)),PETSC_ERR_SUP,"Only valid for DMSwarm-PIC. You must call DMSwarmSetType(dm,DMSWARM_PIC)"); \
     15:   else \
     16:     if (!_swarm->dmcell) SETERRQ(PetscObjectComm((PetscObject)(dm)),PETSC_ERR_SUP,"Only valid for DMSwarmPIC if the cell DM is set. You must call DMSwarmSetCellDM(dm,celldm)"); \
     17: }

     19: /* Coordinate insertition/addition API */
     20: /*@C
     21:    DMSwarmSetPointsUniformCoordinates - Set point coordinates in a DMSwarm on a regular (ijk) grid

     23:    Collective on dm

     25:    Input parameters:
     26: +  dm - the DMSwarm
     27: .  min - minimum coordinate values in the x, y, z directions (array of length dim)
     28: .  max - maximum coordinate values in the x, y, z directions (array of length dim)
     29: .  npoints - number of points in each spatial direction (array of length dim)
     30: -  mode - indicates whether to append points to the swarm (ADD_VALUES), or over-ride existing points (INSERT_VALUES)

     32:    Level: beginner

     34:    Notes:
     35:    When using mode = INSERT_VALUES, this method will reset the number of particles in the DMSwarm
     36:    to be npoints[0]*npoints[1] (2D) or npoints[0]*npoints[1]*npoints[2] (3D). When using mode = ADD_VALUES,
     37:    new points will be appended to any already existing in the DMSwarm

     39: .seealso: DMSwarmSetType(), DMSwarmSetCellDM(), DMSwarmType
     40: @*/
     41: PETSC_EXTERN PetscErrorCode DMSwarmSetPointsUniformCoordinates(DM dm,PetscReal min[],PetscReal max[],PetscInt npoints[],InsertMode mode)
     42: {
     43:   PetscErrorCode    ierr;
     44:   PetscReal         gmin[] = {PETSC_MAX_REAL ,PETSC_MAX_REAL, PETSC_MAX_REAL};
     45:   PetscReal         gmax[] = {PETSC_MIN_REAL, PETSC_MIN_REAL, PETSC_MIN_REAL};
     46:   PetscInt          i,j,k,N,bs,b,n_estimate,n_curr,n_new_est,p,n_found;
     47:   Vec               coorlocal;
     48:   const PetscScalar *_coor;
     49:   DM                celldm;
     50:   PetscReal         dx[3];
     51:   PetscInt          _npoints[] = { 0, 0, 1 };
     52:   Vec               pos;
     53:   PetscScalar       *_pos;
     54:   PetscReal         *swarm_coor;
     55:   PetscInt          *swarm_cellid;
     56:   PetscSF           sfcell = NULL;
     57:   const PetscSFNode *LA_sfcell;

     60:   DMSWARMPICVALID(dm);
     61:   DMSwarmGetCellDM(dm,&celldm);
     62:   DMGetCoordinatesLocal(celldm,&coorlocal);
     63:   VecGetSize(coorlocal,&N);
     64:   VecGetBlockSize(coorlocal,&bs);
     65:   N = N / bs;
     66:   VecGetArrayRead(coorlocal,&_coor);
     67:   for (i=0; i<N; i++) {
     68:     for (b=0; b<bs; b++) {
     69:       gmin[b] = PetscMin(gmin[b],PetscRealPart(_coor[bs*i+b]));
     70:       gmax[b] = PetscMax(gmax[b],PetscRealPart(_coor[bs*i+b]));
     71:     }
     72:   }
     73:   VecRestoreArrayRead(coorlocal,&_coor);

     75:   for (b=0; b<bs; b++) {
     76:     if (npoints[b] > 1) {
     77:       dx[b] = (max[b] - min[b])/((PetscReal)(npoints[b]-1));
     78:     } else {
     79:       dx[b] = 0.0;
     80:     }

     82:     _npoints[b] = npoints[b];
     83:   }

     85:   /* determine number of points living in the bounding box */
     86:   n_estimate = 0;
     87:   for (k=0; k<_npoints[2]; k++) {
     88:     for (j=0; j<_npoints[1]; j++) {
     89:       for (i=0; i<_npoints[0]; i++) {
     90:         PetscReal xp[] = {0.0,0.0,0.0};
     91:         PetscInt ijk[3];
     92:         PetscBool point_inside = PETSC_TRUE;

     94:         ijk[0] = i;
     95:         ijk[1] = j;
     96:         ijk[2] = k;
     97:         for (b=0; b<bs; b++) {
     98:           xp[b] = min[b] + ijk[b] * dx[b];
     99:         }
    100:         for (b=0; b<bs; b++) {
    101:           if (xp[b] < gmin[b]) { point_inside = PETSC_FALSE; }
    102:           if (xp[b] > gmax[b]) { point_inside = PETSC_FALSE; }
    103:         }
    104:         if (point_inside) { n_estimate++; }
    105:       }
    106:     }
    107:   }

    109:   /* create candidate list */
    110:   VecCreate(PETSC_COMM_SELF,&pos);
    111:   VecSetSizes(pos,bs*n_estimate,PETSC_DECIDE);
    112:   VecSetBlockSize(pos,bs);
    113:   VecSetFromOptions(pos);
    114:   VecGetArray(pos,&_pos);

    116:   n_estimate = 0;
    117:   for (k=0; k<_npoints[2]; k++) {
    118:     for (j=0; j<_npoints[1]; j++) {
    119:       for (i=0; i<_npoints[0]; i++) {
    120:         PetscReal xp[] = {0.0,0.0,0.0};
    121:         PetscInt  ijk[3];
    122:         PetscBool point_inside = PETSC_TRUE;

    124:         ijk[0] = i;
    125:         ijk[1] = j;
    126:         ijk[2] = k;
    127:         for (b=0; b<bs; b++) {
    128:           xp[b] = min[b] + ijk[b] * dx[b];
    129:         }
    130:         for (b=0; b<bs; b++) {
    131:           if (xp[b] < gmin[b]) { point_inside = PETSC_FALSE; }
    132:           if (xp[b] > gmax[b]) { point_inside = PETSC_FALSE; }
    133:         }
    134:         if (point_inside) {
    135:           for (b=0; b<bs; b++) {
    136:             _pos[bs*n_estimate+b] = xp[b];
    137:           }
    138:           n_estimate++;
    139:         }
    140:       }
    141:     }
    142:   }
    143:   VecRestoreArray(pos,&_pos);

    145:   /* locate points */
    146:   DMLocatePoints(celldm,pos,DM_POINTLOCATION_NONE,&sfcell);

    148:   PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell);
    149:   n_found = 0;
    150:   for (p=0; p<n_estimate; p++) {
    151:     if (LA_sfcell[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) {
    152:       n_found++;
    153:     }
    154:   }

    156:   /* adjust size */
    157:   if (mode == ADD_VALUES) {
    158:     DMSwarmGetLocalSize(dm,&n_curr);
    159:     n_new_est = n_curr + n_found;
    160:     DMSwarmSetLocalSizes(dm,n_new_est,-1);
    161:   }
    162:   if (mode == INSERT_VALUES) {
    163:     n_curr = 0;
    164:     n_new_est = n_found;
    165:     DMSwarmSetLocalSizes(dm,n_new_est,-1);
    166:   }

    168:   /* initialize new coords, cell owners, pid */
    169:   VecGetArrayRead(pos,&_coor);
    170:   DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
    171:   DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
    172:   n_found = 0;
    173:   for (p=0; p<n_estimate; p++) {
    174:     if (LA_sfcell[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) {
    175:       for (b=0; b<bs; b++) {
    176:         swarm_coor[bs*(n_curr + n_found) + b] = PetscRealPart(_coor[bs*p+b]);
    177:       }
    178:       swarm_cellid[n_curr + n_found] = LA_sfcell[p].index;
    179:       n_found++;
    180:     }
    181:   }
    182:   DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
    183:   DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
    184:   VecRestoreArrayRead(pos,&_coor);

    186:   PetscSFDestroy(&sfcell);
    187:   VecDestroy(&pos);

    189:   return(0);
    190: }

    192: /*@C
    193:    DMSwarmSetPointCoordinates - Set point coordinates in a DMSwarm from a user defined list

    195:    Collective on dm

    197:    Input parameters:
    198: +  dm - the DMSwarm
    199: .  npoints - the number of points to insert
    200: .  coor - the coordinate values
    201: .  redundant - if set to PETSC_TRUE, it is assumed that npoints and coor[] are only valid on rank 0 and should be broadcast to other ranks
    202: -  mode - indicates whether to append points to the swarm (ADD_VALUES), or over-ride existing points (INSERT_VALUES)

    204:    Level: beginner

    206:    Notes:
    207:    If the user has specified redundant = PETSC_FALSE, the cell DM will attempt to locate the coordinates provided by coor[] within
    208:    its sub-domain. If they any values within coor[] are not located in the sub-domain, they will be ignored and will not get
    209:    added to the DMSwarm.

    211: .seealso: DMSwarmSetType(), DMSwarmSetCellDM(), DMSwarmType, DMSwarmSetPointsUniformCoordinates()
    212: @*/
    213: PETSC_EXTERN PetscErrorCode DMSwarmSetPointCoordinates(DM dm,PetscInt npoints,PetscReal coor[],PetscBool redundant,InsertMode mode)
    214: {
    215:   PetscErrorCode    ierr;
    216:   PetscReal         gmin[] = {PETSC_MAX_REAL ,PETSC_MAX_REAL, PETSC_MAX_REAL};
    217:   PetscReal         gmax[] = {PETSC_MIN_REAL, PETSC_MIN_REAL, PETSC_MIN_REAL};
    218:   PetscInt          i,N,bs,b,n_estimate,n_curr,n_new_est,p,n_found;
    219:   Vec               coorlocal;
    220:   const PetscScalar *_coor;
    221:   DM                celldm;
    222:   Vec               pos;
    223:   PetscScalar       *_pos;
    224:   PetscReal         *swarm_coor;
    225:   PetscInt          *swarm_cellid;
    226:   PetscSF           sfcell = NULL;
    227:   const PetscSFNode *LA_sfcell;
    228:   PetscReal         *my_coor;
    229:   PetscInt          my_npoints;
    230:   PetscMPIInt       rank;
    231:   MPI_Comm          comm;

    234:   DMSWARMPICVALID(dm);
    235:   PetscObjectGetComm((PetscObject)dm,&comm);
    236:   MPI_Comm_rank(comm,&rank);

    238:   DMSwarmGetCellDM(dm,&celldm);
    239:   DMGetCoordinatesLocal(celldm,&coorlocal);
    240:   VecGetSize(coorlocal,&N);
    241:   VecGetBlockSize(coorlocal,&bs);
    242:   N = N / bs;
    243:   VecGetArrayRead(coorlocal,&_coor);
    244:   for (i=0; i<N; i++) {
    245:     for (b=0; b<bs; b++) {
    246:       gmin[b] = PetscMin(gmin[b],PetscRealPart(_coor[bs*i+b]));
    247:       gmax[b] = PetscMax(gmax[b],PetscRealPart(_coor[bs*i+b]));
    248:     }
    249:   }
    250:   VecRestoreArrayRead(coorlocal,&_coor);

    252:   /* broadcast points from rank 0 if requested */
    253:   if (redundant) {
    254:     my_npoints = npoints;
    255:     MPI_Bcast(&my_npoints,1,MPIU_INT,0,comm);

    257:     if (rank > 0) { /* allocate space */
    258:       PetscMalloc1(bs*my_npoints,&my_coor);
    259:     } else {
    260:       my_coor = coor;
    261:     }
    262:     MPI_Bcast(my_coor,bs*my_npoints,MPIU_REAL,0,comm);
    263:   } else {
    264:     my_npoints = npoints;
    265:     my_coor = coor;
    266:   }

    268:   /* determine the number of points living in the bounding box */
    269:   n_estimate = 0;
    270:   for (i=0; i<my_npoints; i++) {
    271:     PetscBool point_inside = PETSC_TRUE;

    273:     for (b=0; b<bs; b++) {
    274:       if (my_coor[bs*i+b] < gmin[b]) { point_inside = PETSC_FALSE; }
    275:       if (my_coor[bs*i+b] > gmax[b]) { point_inside = PETSC_FALSE; }
    276:     }
    277:     if (point_inside) { n_estimate++; }
    278:   }

    280:   /* create candidate list */
    281:   VecCreate(PETSC_COMM_SELF,&pos);
    282:   VecSetSizes(pos,bs*n_estimate,PETSC_DECIDE);
    283:   VecSetBlockSize(pos,bs);
    284:   VecSetFromOptions(pos);
    285:   VecGetArray(pos,&_pos);

    287:   n_estimate = 0;
    288:   for (i=0; i<my_npoints; i++) {
    289:     PetscBool point_inside = PETSC_TRUE;

    291:     for (b=0; b<bs; b++) {
    292:       if (my_coor[bs*i+b] < gmin[b]) { point_inside = PETSC_FALSE; }
    293:       if (my_coor[bs*i+b] > gmax[b]) { point_inside = PETSC_FALSE; }
    294:     }
    295:     if (point_inside) {
    296:       for (b=0; b<bs; b++) {
    297:         _pos[bs*n_estimate+b] = my_coor[bs*i+b];
    298:       }
    299:       n_estimate++;
    300:     }
    301:   }
    302:   VecRestoreArray(pos,&_pos);

    304:   /* locate points */
    305:   DMLocatePoints(celldm,pos,DM_POINTLOCATION_NONE,&sfcell);

    307:   PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell);
    308:   n_found = 0;
    309:   for (p=0; p<n_estimate; p++) {
    310:     if (LA_sfcell[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) {
    311:       n_found++;
    312:     }
    313:   }

    315:   /* adjust size */
    316:   if (mode == ADD_VALUES) {
    317:     DMSwarmGetLocalSize(dm,&n_curr);
    318:     n_new_est = n_curr + n_found;
    319:     DMSwarmSetLocalSizes(dm,n_new_est,-1);
    320:   }
    321:   if (mode == INSERT_VALUES) {
    322:     n_curr = 0;
    323:     n_new_est = n_found;
    324:     DMSwarmSetLocalSizes(dm,n_new_est,-1);
    325:   }

    327:   /* initialize new coords, cell owners, pid */
    328:   VecGetArrayRead(pos,&_coor);
    329:   DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
    330:   DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
    331:   n_found = 0;
    332:   for (p=0; p<n_estimate; p++) {
    333:     if (LA_sfcell[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) {
    334:       for (b=0; b<bs; b++) {
    335:         swarm_coor[bs*(n_curr + n_found) + b] = PetscRealPart(_coor[bs*p+b]);
    336:       }
    337:       swarm_cellid[n_curr + n_found] = LA_sfcell[p].index;
    338:       n_found++;
    339:     }
    340:   }
    341:   DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
    342:   DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
    343:   VecRestoreArrayRead(pos,&_coor);

    345:   if (redundant) {
    346:     if (rank > 0) {
    347:       PetscFree(my_coor);
    348:     }
    349:   }
    350:   PetscSFDestroy(&sfcell);
    351:   VecDestroy(&pos);

    353:   return(0);
    354: }

    356: extern PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_DA(DM,DM,DMSwarmPICLayoutType,PetscInt);
    357: extern PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX(DM,DM,DMSwarmPICLayoutType,PetscInt);

    359: /*@C
    360:    DMSwarmInsertPointsUsingCellDM - Insert point coordinates within each cell

    362:    Not collective

    364:    Input parameters:
    365: +  dm - the DMSwarm
    366: .  layout_type - method used to fill each cell with the cell DM
    367: -  fill_param - parameter controlling how many points per cell are added (the meaning of this parameter is dependent on the layout type)

    369:    Level: beginner

    371:    Notes:

    373:    The insert method will reset any previous defined points within the DMSwarm.

    375:    When using a DMDA both 2D and 3D are supported for all layout types provided you are using DMDA_ELEMENT_Q1.

    377:    When using a DMPLEX the following case are supported:
    378:    (i) DMSWARMPIC_LAYOUT_REGULAR: 2D (triangle),
    379:    (ii) DMSWARMPIC_LAYOUT_GAUSS: 2D and 3D provided the cell is a tri/tet or a quad/hex,
    380:    (iii) DMSWARMPIC_LAYOUT_SUBDIVISION: 2D and 3D for quad/hex and 2D tri.

    382: .seealso: DMSwarmPICLayoutType, DMSwarmSetType(), DMSwarmSetCellDM(), DMSwarmType
    383: @*/
    384: PETSC_EXTERN PetscErrorCode DMSwarmInsertPointsUsingCellDM(DM dm,DMSwarmPICLayoutType layout_type,PetscInt fill_param)
    385: {
    387:   DM             celldm;
    388:   PetscBool      isDA,isPLEX;

    391:   DMSWARMPICVALID(dm);
    392:   DMSwarmGetCellDM(dm,&celldm);
    393:   PetscObjectTypeCompare((PetscObject)celldm,DMDA,&isDA);
    394:   PetscObjectTypeCompare((PetscObject)celldm,DMPLEX,&isPLEX);
    395:   if (isDA) {
    396:     private_DMSwarmInsertPointsUsingCellDM_DA(dm,celldm,layout_type,fill_param);
    397:   } else if (isPLEX) {
    398:     private_DMSwarmInsertPointsUsingCellDM_PLEX(dm,celldm,layout_type,fill_param);
    399:   } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only supported for cell DMs of type DMDA and DMPLEX");

    401:   return(0);
    402: }

    405: extern PetscErrorCode private_DMSwarmSetPointCoordinatesCellwise_PLEX(DM,DM,PetscInt,PetscReal*);

    407: /*@C
    408:    DMSwarmSetPointCoordinatesCellwise - Insert point coordinates (defined over the reference cell) within each cell

    410:    Not collective

    412:    Input parameters:
    413: +  dm - the DMSwarm
    414: .  celldm - the cell DM
    415: .  npoints - the number of points to insert in each cell
    416: -  xi - the coordinates (defined in the local coordinate system for each cell) to insert

    418:  Level: beginner

    420:  Notes:
    421:  The method will reset any previous defined points within the DMSwarm.
    422:  Only supported for DMPLEX. If you are using a DMDA it is recommended to either use
    423:  DMSwarmInsertPointsUsingCellDM(), or extract and set the coordinates yourself the following code

    425: $    PetscReal *coor;
    426: $    DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&coor);
    427: $    // user code to define the coordinates here
    428: $    DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&coor);

    430: .seealso: DMSwarmSetCellDM(), DMSwarmInsertPointsUsingCellDM()
    431: @*/
    432: PETSC_EXTERN PetscErrorCode DMSwarmSetPointCoordinatesCellwise(DM dm,PetscInt npoints,PetscReal xi[])
    433: {
    435:   DM             celldm;
    436:   PetscBool      isDA,isPLEX;

    439:   DMSWARMPICVALID(dm);
    440:   DMSwarmGetCellDM(dm,&celldm);
    441:   PetscObjectTypeCompare((PetscObject)celldm,DMDA,&isDA);
    442:   PetscObjectTypeCompare((PetscObject)celldm,DMPLEX,&isPLEX);
    443:   if (isDA) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only supported for cell DMs of type DMPLEX. Recommended you use DMSwarmInsertPointsUsingCellDM()");
    444:   else if (isPLEX) {
    445:     private_DMSwarmSetPointCoordinatesCellwise_PLEX(dm,celldm,npoints,xi);
    446:   } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only supported for cell DMs of type DMDA and DMPLEX");

    448:   return(0);
    449: }

    452: /* Field projection API */
    453: extern PetscErrorCode private_DMSwarmProjectFields_DA(DM swarm,DM celldm,PetscInt project_type,PetscInt nfields,DMSwarmDataField dfield[],Vec vecs[]);
    454: extern PetscErrorCode private_DMSwarmProjectFields_PLEX(DM swarm,DM celldm,PetscInt project_type,PetscInt nfields,DMSwarmDataField dfield[],Vec vecs[]);

    456: /*@C
    457:    DMSwarmProjectFields - Project a set of swarm fields onto the cell DM

    459:    Collective on dm

    461:    Input parameters:
    462: +  dm - the DMSwarm
    463: .  nfields - the number of swarm fields to project
    464: .  fieldnames - the textual names of the swarm fields to project
    465: .  fields - an array of Vec's of length nfields
    466: -  reuse - flag indicating whether the array and contents of fields should be re-used or internally allocated

    468:    Currently, the only available projection method consists of
    469:      phi_i = \sum_{p=0}^{np} N_i(x_p) phi_p dJ / \sum_{p=0}^{np} N_i(x_p) dJ
    470:    where phi_p is the swarm field at point p,
    471:      N_i() is the cell DM basis function at vertex i,
    472:      dJ is the determinant of the cell Jacobian and
    473:      phi_i is the projected vertex value of the field phi.

    475:    Level: beginner

    477:    Notes:

    479:    If reuse = PETSC_FALSE, this function will allocate the array of Vec's, and each individual Vec.
    480:      The user is responsible for destroying both the array and the individual Vec objects.

    482:    Only swarm fields registered with data type = PETSC_REAL can be projected onto the cell DM.

    484:    Only swarm fields of block size = 1 can currently be projected.

    486:    The only projection methods currently only support the DA (2D) and PLEX (triangles 2D).

    488: .seealso: DMSwarmSetType(), DMSwarmSetCellDM(), DMSwarmType
    489: @*/
    490: PETSC_EXTERN PetscErrorCode DMSwarmProjectFields(DM dm,PetscInt nfields,const char *fieldnames[],Vec **fields,PetscBool reuse)
    491: {
    492:   DM_Swarm         *swarm = (DM_Swarm*)dm->data;
    493:   DMSwarmDataField *gfield;
    494:   DM               celldm;
    495:   PetscBool        isDA,isPLEX;
    496:   Vec              *vecs;
    497:   PetscInt         f,nvecs;
    498:   PetscInt         project_type = 0;

    502:   DMSWARMPICVALID(dm);
    503:   DMSwarmGetCellDM(dm,&celldm);
    504:   PetscMalloc1(nfields,&gfield);
    505:   nvecs = 0;
    506:   for (f=0; f<nfields; f++) {
    507:     DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,fieldnames[f],&gfield[f]);
    508:     if (gfield[f]->petsc_type != PETSC_REAL) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Projection only valid for fields using a data type = PETSC_REAL");
    509:     if (gfield[f]->bs != 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Projection only valid for fields with block size = 1");
    510:     nvecs += gfield[f]->bs;
    511:   }
    512:   if (!reuse) {
    513:     PetscMalloc1(nvecs,&vecs);
    514:     for (f=0; f<nvecs; f++) {
    515:       DMCreateGlobalVector(celldm,&vecs[f]);
    516:       PetscObjectSetName((PetscObject)vecs[f],gfield[f]->name);
    517:     }
    518:   } else {
    519:     vecs = *fields;
    520:   }

    522:   PetscObjectTypeCompare((PetscObject)celldm,DMDA,&isDA);
    523:   PetscObjectTypeCompare((PetscObject)celldm,DMPLEX,&isPLEX);
    524:   if (isDA) {
    525:     private_DMSwarmProjectFields_DA(dm,celldm,project_type,nfields,gfield,vecs);
    526:   } else if (isPLEX) {
    527:     private_DMSwarmProjectFields_PLEX(dm,celldm,project_type,nfields,gfield,vecs);
    528:   } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only supported for cell DMs of type DMDA and DMPLEX");

    530:   PetscFree(gfield);
    531:   if (!reuse) {
    532:     *fields = vecs;
    533:   }

    535:   return(0);
    536: }

    538: /*@C
    539:    DMSwarmCreatePointPerCellCount - Count the number of points within all cells in the cell DM

    541:    Not collective

    543:    Input parameter:
    544: .  dm - the DMSwarm

    546:    Output parameters:
    547: +  ncells - the number of cells in the cell DM (optional argument, pass NULL to ignore)
    548: -  count - array of length ncells containing the number of points per cell

    550:    Level: beginner

    552:    Notes:
    553:    The array count is allocated internally and must be free'd by the user.

    555: .seealso: DMSwarmSetType(), DMSwarmSetCellDM(), DMSwarmType
    556: @*/
    557: PETSC_EXTERN PetscErrorCode DMSwarmCreatePointPerCellCount(DM dm,PetscInt *ncells,PetscInt **count)
    558: {
    560:   PetscBool      isvalid;
    561:   PetscInt       nel;
    562:   PetscInt       *sum;

    565:   DMSwarmSortGetIsValid(dm,&isvalid);
    566:   nel = 0;
    567:   if (isvalid) {
    568:     PetscInt e;

    570:     DMSwarmSortGetSizes(dm,&nel,NULL);

    572:     PetscMalloc1(nel,&sum);
    573:     for (e=0; e<nel; e++) {
    574:       DMSwarmSortGetNumberOfPointsPerCell(dm,e,&sum[e]);
    575:     }
    576:   } else {
    577:     DM        celldm;
    578:     PetscBool isda,isplex,isshell;
    579:     PetscInt  p,npoints;
    580:     PetscInt *swarm_cellid;

    582:     /* get the number of cells */
    583:     DMSwarmGetCellDM(dm,&celldm);
    584:     PetscObjectTypeCompare((PetscObject)celldm,DMDA,&isda);
    585:     PetscObjectTypeCompare((PetscObject)celldm,DMPLEX,&isplex);
    586:     PetscObjectTypeCompare((PetscObject)celldm,DMSHELL,&isshell);
    587:     if (isda) {
    588:       PetscInt       _nel,_npe;
    589:       const PetscInt *_element;

    591:       DMDAGetElements(celldm,&_nel,&_npe,&_element);
    592:       nel = _nel;
    593:       DMDARestoreElements(celldm,&_nel,&_npe,&_element);
    594:     } else if (isplex) {
    595:       PetscInt ps,pe;

    597:       DMPlexGetHeightStratum(celldm,0,&ps,&pe);
    598:       nel = pe - ps;
    599:     } else if (isshell) {
    600:       PetscErrorCode (*method_DMShellGetNumberOfCells)(DM,PetscInt*);

    602:       PetscObjectQueryFunction((PetscObject)celldm,"DMGetNumberOfCells_C",&method_DMShellGetNumberOfCells);
    603:       if (method_DMShellGetNumberOfCells) {
    604:         method_DMShellGetNumberOfCells(celldm,&nel);
    605:       } 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);");
    606:     } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Cannot determine the number of cells for a DM not of type DA, PLEX or SHELL");

    608:     PetscMalloc1(nel,&sum);
    609:     PetscArrayzero(sum,nel);
    610:     DMSwarmGetLocalSize(dm,&npoints);
    611:     DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
    612:     for (p=0; p<npoints; p++) {
    613:       if (swarm_cellid[p] != DMLOCATEPOINT_POINT_NOT_FOUND) {
    614:         sum[ swarm_cellid[p] ]++;
    615:       }
    616:     }
    617:     DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
    618:   }
    619:   if (ncells) { *ncells = nel; }
    620:   *count  = sum;
    621:   return(0);
    622: }