Actual source code: swarmpic.c

petsc-3.11.4 2019-09-28
Report Typos and Errors
  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
 22:  
 23:    Collective on DM
 24:  
 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)
 31:  
 32:    Level: beginner
 33:  
 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
 38:  
 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;
 58: 
 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:     }
 81: 
 82:     _npoints[b] = npoints[b];
 83:   }
 84: 
 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;
 93: 
 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:   }
108: 
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);
115: 
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;
123: 
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);
144: 
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:   }
167: 
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);
188: 
189:   return(0);
190: }

192: /*@C
193:    DMSwarmSetPointCoordinates - Set point coordinates in a DMSwarm from a user defined list
194:  
195:    Collective on DM
196:  
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)
203:  
204:    Level: beginner
205:  
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.
210:  
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;
232: 
234:   DMSWARMPICVALID(dm);
235:   PetscObjectGetComm((PetscObject)dm,&comm);
236:   MPI_Comm_rank(comm,&rank);
237: 
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);
251: 
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:   }
267: 
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;
272: 
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:   }
279: 
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);
286: 
287:   n_estimate = 0;
288:   for (i=0; i<my_npoints; i++) {
289:     PetscBool point_inside = PETSC_TRUE;
290: 
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);
303: 
304:   /* locate points */
305:   DMLocatePoints(celldm,pos,DM_POINTLOCATION_NONE,&sfcell);
306: 
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:   }
314: 
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:   }
326: 
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);
344: 
345:   if (redundant) {
346:     if (rank > 0) {
347:       PetscFree(my_coor);
348:     }
349:   }
350:   PetscSFDestroy(&sfcell);
351:   VecDestroy(&pos);
352: 
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
361:  
362:    Not collective
363:  
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");
400: 
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
424:  
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;
437: 
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");
447: 
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
458:  
459:    Collective on Vec
460:  
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
467:  
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.
474:  
475:    Level: beginner
476:  
477:    Notes:
478:  
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.
481:  
482:    Only swarm fields registered with data type = PETSC_REAL can be projected onto the cell DM.
483:  
484:    Only swarm fields of block size = 1 can currently be projected.
485:  
486:    The only projection methods currently only support the DA (2D) and PLEX (triangles 2D).
487:  
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;
500: 
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:   }
521: 
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");
529: 
530:   PetscFree(gfield);
531:   if (!reuse) {
532:     *fields = vecs;
533:   }
534: 
535:   return(0);
536: }

538: /*@C
539:    DMSwarmCreatePointPerCellCount - Count the number of points within all cells in the cell DM
540:  
541:    Not collective
542:  
543:    Input parameter:
544: .  dm - the DMSwarm
545:  
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
549:  
550:    Level: beginner
551:  
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;
563: 
565:   DMSwarmSortGetIsValid(dm,&isvalid);
566:   nel = 0;
567:   if (isvalid) {
568:     PetscInt e;
569: 
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;
590: 
591:       DMDAGetElements(celldm,&_nel,&_npe,&_element);
592:       nel = _nel;
593:       DMDARestoreElements(celldm,&_nel,&_npe,&_element);
594:     } else if (isplex) {
595:       PetscInt ps,pe;
596: 
597:       DMPlexGetHeightStratum(celldm,0,&ps,&pe);
598:       nel = pe - ps;
599:     } else if (isshell) {
600:       PetscErrorCode (*method_DMShellGetNumberOfCells)(DM,PetscInt*);
601: 
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");
607: 
608:     PetscMalloc1(nel,&sum);
609:     PetscMemzero(sum,sizeof(PetscInt)*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: }