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 <petscdt.h>
  7: #include "../src/dm/impls/swarm/data_bucket.h"

  9: #include <petsc/private/petscfeimpl.h>

 11: /*
 12:  Error checking to ensure the swarm type is correct and that a cell DM has been set
 13: */
 14: #define DMSWARMPICVALID(dm) \
 15: { \
 16:   DM_Swarm *_swarm = (DM_Swarm*)(dm)->data; \
 18:   else \
 20: }

 22: /* Coordinate insertition/addition API */
 23: /*@C
 24:    DMSwarmSetPointsUniformCoordinates - Set point coordinates in a DMSwarm on a regular (ijk) grid

 26:    Collective on dm

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

 35:    Level: beginner

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

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

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

 76:   for (b=0; b<bs; b++) {
 77:     if (npoints[b] > 1) {
 78:       dx[b] = (max[b] - min[b])/((PetscReal)(npoints[b]-1));
 79:     } else {
 80:       dx[b] = 0.0;
 81:     }
 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);
147:   PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell);
148:   n_found = 0;
149:   for (p=0; p<n_estimate; p++) {
150:     if (LA_sfcell[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) {
151:       n_found++;
152:     }
153:   }

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

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

185:   PetscSFDestroy(&sfcell);
186:   VecDestroy(&pos);
187:   return 0;
188: }

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

193:    Collective on dm

195:    Input parameters:
196: +  dm - the DMSwarm
197: .  npoints - the number of points to insert
198: .  coor - the coordinate values
199: .  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
200: -  mode - indicates whether to append points to the swarm (ADD_VALUES), or over-ride existing points (INSERT_VALUES)

202:    Level: beginner

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

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

230:   DMSWARMPICVALID(dm);
231:   PetscObjectGetComm((PetscObject)dm,&comm);
232:   MPI_Comm_rank(comm,&rank);

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

248:   /* broadcast points from rank 0 if requested */
249:   if (redundant) {
250:     my_npoints = npoints;
251:     MPI_Bcast(&my_npoints,1,MPIU_INT,0,comm);

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

264:   /* determine the number of points living in the bounding box */
265:   n_estimate = 0;
266:   for (i=0; i<my_npoints; i++) {
267:     PetscBool point_inside = PETSC_TRUE;

269:     for (b=0; b<bs; b++) {
270:       if (my_coor[bs*i+b] < gmin[b]) { point_inside = PETSC_FALSE; }
271:       if (my_coor[bs*i+b] > gmax[b]) { point_inside = PETSC_FALSE; }
272:     }
273:     if (point_inside) { n_estimate++; }
274:   }

276:   /* create candidate list */
277:   VecCreate(PETSC_COMM_SELF,&pos);
278:   VecSetSizes(pos,bs*n_estimate,PETSC_DECIDE);
279:   VecSetBlockSize(pos,bs);
280:   VecSetFromOptions(pos);
281:   VecGetArray(pos,&_pos);

283:   n_estimate = 0;
284:   for (i=0; i<my_npoints; i++) {
285:     PetscBool point_inside = PETSC_TRUE;

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

300:   /* locate points */
301:   DMLocatePoints(celldm,pos,DM_POINTLOCATION_NONE,&sfcell);

303:   PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell);
304:   n_found = 0;
305:   for (p=0; p<n_estimate; p++) {
306:     if (LA_sfcell[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) {
307:       n_found++;
308:     }
309:   }

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

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

341:   if (redundant) {
342:     if (rank > 0) {
343:       PetscFree(my_coor);
344:     }
345:   }
346:   PetscSFDestroy(&sfcell);
347:   VecDestroy(&pos);
348:   return 0;
349: }

351: extern PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_DA(DM,DM,DMSwarmPICLayoutType,PetscInt);
352: extern PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX(DM,DM,DMSwarmPICLayoutType,PetscInt);

354: /*@C
355:    DMSwarmInsertPointsUsingCellDM - Insert point coordinates within each cell

357:    Not collective

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

364:    Level: beginner

366:    Notes:

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

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

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

377: .seealso: DMSwarmPICLayoutType, DMSwarmSetType(), DMSwarmSetCellDM(), DMSwarmType
378: @*/
379: PETSC_EXTERN PetscErrorCode DMSwarmInsertPointsUsingCellDM(DM dm,DMSwarmPICLayoutType layout_type,PetscInt fill_param)
380: {
381:   DM             celldm;
382:   PetscBool      isDA,isPLEX;

384:   DMSWARMPICVALID(dm);
385:   DMSwarmGetCellDM(dm,&celldm);
386:   PetscObjectTypeCompare((PetscObject)celldm,DMDA,&isDA);
387:   PetscObjectTypeCompare((PetscObject)celldm,DMPLEX,&isPLEX);
388:   if (isDA) {
389:     private_DMSwarmInsertPointsUsingCellDM_DA(dm,celldm,layout_type,fill_param);
390:   } else if (isPLEX) {
391:     private_DMSwarmInsertPointsUsingCellDM_PLEX(dm,celldm,layout_type,fill_param);
392:   } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only supported for cell DMs of type DMDA and DMPLEX");
393:   return 0;
394: }

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

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

401:    Not collective

403:    Input parameters:
404: +  dm - the DMSwarm
405: .  celldm - the cell DM
406: .  npoints - the number of points to insert in each cell
407: -  xi - the coordinates (defined in the local coordinate system for each cell) to insert

409:  Level: beginner

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

416: $    PetscReal *coor;
417: $    DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&coor);
418: $    // user code to define the coordinates here
419: $    DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&coor);

421: .seealso: DMSwarmSetCellDM(), DMSwarmInsertPointsUsingCellDM()
422: @*/
423: PETSC_EXTERN PetscErrorCode DMSwarmSetPointCoordinatesCellwise(DM dm,PetscInt npoints,PetscReal xi[])
424: {
425:   DM             celldm;
426:   PetscBool      isDA,isPLEX;

428:   DMSWARMPICVALID(dm);
429:   DMSwarmGetCellDM(dm,&celldm);
430:   PetscObjectTypeCompare((PetscObject)celldm,DMDA,&isDA);
431:   PetscObjectTypeCompare((PetscObject)celldm,DMPLEX,&isPLEX);
433:   else if (isPLEX) {
434:     private_DMSwarmSetPointCoordinatesCellwise_PLEX(dm,celldm,npoints,xi);
435:   } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only supported for cell DMs of type DMDA and DMPLEX");
436:   return 0;
437: }

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

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

446:    Collective on dm

448:    Input parameters:
449: +  dm - the DMSwarm
450: .  nfields - the number of swarm fields to project
451: .  fieldnames - the textual names of the swarm fields to project
452: .  fields - an array of Vec's of length nfields
453: -  reuse - flag indicating whether the array and contents of fields should be re-used or internally allocated

455:    Currently, the only available projection method consists of
456:      phi_i = \sum_{p=0}^{np} N_i(x_p) phi_p dJ / \sum_{p=0}^{np} N_i(x_p) dJ
457:    where phi_p is the swarm field at point p,
458:      N_i() is the cell DM basis function at vertex i,
459:      dJ is the determinant of the cell Jacobian and
460:      phi_i is the projected vertex value of the field phi.

462:    Level: beginner

464:    Notes:

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

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

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

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

475: .seealso: DMSwarmSetType(), DMSwarmSetCellDM(), DMSwarmType
476: @*/
477: PETSC_EXTERN PetscErrorCode DMSwarmProjectFields(DM dm,PetscInt nfields,const char *fieldnames[],Vec **fields,PetscBool reuse)
478: {
479:   DM_Swarm         *swarm = (DM_Swarm*)dm->data;
480:   DMSwarmDataField *gfield;
481:   DM               celldm;
482:   PetscBool        isDA,isPLEX;
483:   Vec              *vecs;
484:   PetscInt         f,nvecs;
485:   PetscInt         project_type = 0;

487:   DMSWARMPICVALID(dm);
488:   DMSwarmGetCellDM(dm,&celldm);
489:   PetscMalloc1(nfields,&gfield);
490:   nvecs = 0;
491:   for (f=0; f<nfields; f++) {
492:     DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,fieldnames[f],&gfield[f]);
495:     nvecs += gfield[f]->bs;
496:   }
497:   if (!reuse) {
498:     PetscMalloc1(nvecs,&vecs);
499:     for (f=0; f<nvecs; f++) {
500:       DMCreateGlobalVector(celldm,&vecs[f]);
501:       PetscObjectSetName((PetscObject)vecs[f],gfield[f]->name);
502:     }
503:   } else {
504:     vecs = *fields;
505:   }

507:   PetscObjectTypeCompare((PetscObject)celldm,DMDA,&isDA);
508:   PetscObjectTypeCompare((PetscObject)celldm,DMPLEX,&isPLEX);
509:   if (isDA) {
510:     private_DMSwarmProjectFields_DA(dm,celldm,project_type,nfields,gfield,vecs);
511:   } else if (isPLEX) {
512:     private_DMSwarmProjectFields_PLEX(dm,celldm,project_type,nfields,gfield,vecs);
513:   } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only supported for cell DMs of type DMDA and DMPLEX");

515:   PetscFree(gfield);
516:   if (!reuse) *fields = vecs;
517:   return 0;
518: }

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

523:    Not collective

525:    Input parameter:
526: .  dm - the DMSwarm

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

532:    Level: beginner

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

537: .seealso: DMSwarmSetType(), DMSwarmSetCellDM(), DMSwarmType
538: @*/
539: PETSC_EXTERN PetscErrorCode DMSwarmCreatePointPerCellCount(DM dm,PetscInt *ncells,PetscInt **count)
540: {
541:   PetscBool      isvalid;
542:   PetscInt       nel;
543:   PetscInt       *sum;

545:   DMSwarmSortGetIsValid(dm,&isvalid);
546:   nel = 0;
547:   if (isvalid) {
548:     PetscInt e;

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

552:     PetscMalloc1(nel,&sum);
553:     for (e=0; e<nel; e++) {
554:       DMSwarmSortGetNumberOfPointsPerCell(dm,e,&sum[e]);
555:     }
556:   } else {
557:     DM        celldm;
558:     PetscBool isda,isplex,isshell;
559:     PetscInt  p,npoints;
560:     PetscInt *swarm_cellid;

562:     /* get the number of cells */
563:     DMSwarmGetCellDM(dm,&celldm);
564:     PetscObjectTypeCompare((PetscObject)celldm,DMDA,&isda);
565:     PetscObjectTypeCompare((PetscObject)celldm,DMPLEX,&isplex);
566:     PetscObjectTypeCompare((PetscObject)celldm,DMSHELL,&isshell);
567:     if (isda) {
568:       PetscInt       _nel,_npe;
569:       const PetscInt *_element;

571:       DMDAGetElements(celldm,&_nel,&_npe,&_element);
572:       nel = _nel;
573:       DMDARestoreElements(celldm,&_nel,&_npe,&_element);
574:     } else if (isplex) {
575:       PetscInt ps,pe;

577:       DMPlexGetHeightStratum(celldm,0,&ps,&pe);
578:       nel = pe - ps;
579:     } else if (isshell) {
580:       PetscErrorCode (*method_DMShellGetNumberOfCells)(DM,PetscInt*);

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

588:     PetscMalloc1(nel,&sum);
589:     PetscArrayzero(sum,nel);
590:     DMSwarmGetLocalSize(dm,&npoints);
591:     DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
592:     for (p=0; p<npoints; p++) {
593:       if (swarm_cellid[p] != DMLOCATEPOINT_POINT_NOT_FOUND) {
594:         sum[ swarm_cellid[p] ]++;
595:       }
596:     }
597:     DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
598:   }
599:   if (ncells) { *ncells = nel; }
600:   *count  = sum;
601:   return 0;
602: }

604: /*@
605:   DMSwarmGetNumSpecies - Get the number of particle species

607:   Not collective

609:   Input parameter:
610: . dm - the DMSwarm

612:   Output parameters:
613: . Ns - the number of species

615:   Level: intermediate

617: .seealso: DMSwarmSetNumSpecies(), DMSwarmSetType(), DMSwarmType
618: @*/
619: PetscErrorCode DMSwarmGetNumSpecies(DM sw, PetscInt *Ns)
620: {
621:   DM_Swarm *swarm = (DM_Swarm *) sw->data;

623:   *Ns = swarm->Ns;
624:   return 0;
625: }

627: /*@
628:   DMSwarmSetNumSpecies - Set the number of particle species

630:   Not collective

632:   Input parameter:
633: + dm - the DMSwarm
634: - Ns - the number of species

636:   Level: intermediate

638: .seealso: DMSwarmGetNumSpecies(), DMSwarmSetType(), DMSwarmType
639: @*/
640: PetscErrorCode DMSwarmSetNumSpecies(DM sw, PetscInt Ns)
641: {
642:   DM_Swarm *swarm = (DM_Swarm *) sw->data;

644:   swarm->Ns =  Ns;
645:   return 0;
646: }

648: /*@C
649:   DMSwarmComputeLocalSize - Compute the local number and distribution of particles based upon a density function

651:   Not collective

653:   Input Parameters:
654: + sw      - The DMSwarm
655: . N       - The target number of particles
656: - density - The density field for the particle layout, normalized to unity

658:   Note: One particle will be created for each species.

660:   Level: advanced

662: .seealso: DMSwarmComputeLocalSizeFromOptions()
663: @*/
664: PetscErrorCode DMSwarmComputeLocalSize(DM sw, PetscInt N, PetscProbFunc density)
665: {
666:   DM               dm;
667:   PetscQuadrature  quad;
668:   const PetscReal *xq, *wq;
669:   PetscInt        *npc, *cellid;
670:   PetscReal        xi0[3], scale[1] = {.01};
671:   PetscInt         Ns, cStart, cEnd, c, dim, d, Nq, q, Np = 0, p;
672:   PetscBool        simplex;

674:   DMSwarmGetNumSpecies(sw, &Ns);
675:   DMSwarmGetCellDM(sw, &dm);
676:   DMGetDimension(dm, &dim);
677:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
678:   DMPlexIsSimplex(dm, &simplex);
679:   if (simplex) PetscDTStroudConicalQuadrature(dim, 1, 5, -1.0, 1.0, &quad);
680:   else         PetscDTGaussTensorQuadrature(dim, 1, 5, -1.0, 1.0, &quad);
681:   PetscQuadratureGetData(quad, NULL, NULL, &Nq, &xq, &wq);
682:   PetscMalloc1(cEnd-cStart, &npc);
683:   /* Integrate the density function to get the number of particles in each cell */
684:   for (d = 0; d < dim; ++d) xi0[d] = -1.0;
685:   for (c = 0; c < cEnd-cStart; ++c) {
686:     const PetscInt cell = c + cStart;
687:     PetscReal v0[3], J[9], invJ[9], detJ;
688:     PetscReal n_int = 0.;

690:     DMPlexComputeCellGeometryFEM(dm, cell, NULL, v0, J, invJ, &detJ);
691:     for (q = 0; q < Nq; ++q) {
692:       PetscReal xr[3], den[3];

694:       CoordinatesRefToReal(dim, dim, xi0, v0, J, &xq[q*dim], xr);
695:       density(xr, scale, den);
696:       n_int += N*den[0]*wq[q];
697:     }
698:     npc[c]  = (PetscInt) n_int;
699:     npc[c] *= Ns;
700:     Np     += npc[c];
701:   }
702:   PetscQuadratureDestroy(&quad);
703:   DMSwarmSetLocalSizes(sw, Np, 0);

705:   DMSwarmGetField(sw, DMSwarmPICField_cellid, NULL, NULL, (void **) &cellid);
706:   for (c = 0, p = 0; c < cEnd-cStart; ++c) {
707:     for (q = 0; q < npc[c]; ++q, ++p) cellid[p] = c;
708:   }
709:   DMSwarmRestoreField(sw, DMSwarmPICField_cellid, NULL, NULL, (void **) &cellid);
710:   PetscFree(npc);
711:   return 0;
712: }

714: /*@
715:   DMSwarmComputeLocalSizeFromOptions - Compute the local number and distribution of particles based upon a density function determined by options

717:   Not collective

719:   Input Parameters:
720: , sw - The DMSwarm

722:   Level: advanced

724: .seealso: DMSwarmComputeLocalSize()
725: @*/
726: PetscErrorCode DMSwarmComputeLocalSizeFromOptions(DM sw)
727: {
728:   DTProbDensityType den = DTPROB_DENSITY_CONSTANT;
729:   PetscProbFunc     pdf;
730:   PetscInt          N, Ns, dim;
731:   PetscBool         flg;
732:   const char       *prefix;
733:   PetscErrorCode    ierr;

735:   PetscOptionsBegin(PetscObjectComm((PetscObject) sw), "", "DMSwarm Options", "DMSWARM");
736:   PetscOptionsInt("-dm_swarm_num_particles", "The target number of particles", "", N, &N, NULL);
737:   PetscOptionsInt("-dm_swarm_num_species", "The number of species", "DMSwarmSetNumSpecies", Ns, &Ns, &flg);
738:   if (flg) DMSwarmSetNumSpecies(sw, Ns);
739:   PetscOptionsEnum("-dm_swarm_density", "Method to compute particle density <constant, gaussian>", "", DTProbDensityTypes, (PetscEnum) den, (PetscEnum *) &den, NULL);
740:   PetscOptionsEnd();

742:   DMGetDimension(sw, &dim);
743:   PetscObjectGetOptionsPrefix((PetscObject) sw, &prefix);
744:   PetscProbCreateFromOptions(dim, prefix, "-dm_swarm_coordinate_density", &pdf, NULL, NULL);
745:   DMSwarmComputeLocalSize(sw, N, pdf);
746:   return 0;
747: }

749: /*@
750:   DMSwarmInitializeCoordinates - Determine the initial coordinates of particles for a PIC method

752:   Not collective

754:   Input Parameters:
755: , sw - The DMSwarm

757:   Note: Currently, we randomly place particles in their assigned cell

759:   Level: advanced

761: .seealso: DMSwarmComputeLocalSize(), DMSwarmInitializeVelocities()
762: @*/
763: PetscErrorCode DMSwarmInitializeCoordinates(DM sw)
764: {
765:   DM             dm;
766:   PetscRandom    rnd;
767:   PetscScalar   *weight;
768:   PetscReal     *x, xi0[3];
769:   PetscInt      *species;
770:   PetscBool      removePoints = PETSC_TRUE;
771:   PetscDataType  dtype;
772:   PetscInt       Ns, cStart, cEnd, c, dim, d, s, bs;

775:   DMSwarmGetNumSpecies(sw, &Ns);
776:   DMSwarmGetCellDM(sw, &dm);
777:   DMGetDimension(dm, &dim);
778:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);

780:   /* Set particle position randomly in cell, set weights to 1 */
781:   PetscRandomCreate(PetscObjectComm((PetscObject) dm), &rnd);
782:   PetscRandomSetInterval(rnd, -1.0, 1.0);
783:   PetscRandomSetFromOptions(rnd);
784:   DMSwarmGetField(sw, "w_q", &bs, &dtype, (void **) &weight);
785:   DMSwarmGetField(sw, "DMSwarmPIC_coor", &bs, &dtype, (void **) &x);
786:   DMSwarmGetField(sw, "species", NULL, NULL, (void **) &species);
787:   DMSwarmSortGetAccess(sw);
788:   for (d = 0; d < dim; ++d) xi0[d] = -1.0;
789:   for (c = cStart; c < cEnd; ++c) {
790:     PetscReal v0[3], J[9], invJ[9], detJ;
791:     PetscInt *pidx, Npc, q;

793:     DMSwarmSortGetPointsPerCell(sw, c, &Npc, &pidx);
794:     DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
795:     for (q = 0; q < Npc; ++q) {
796:       const PetscInt p = pidx[q];
797:       PetscReal      xref[3];

799:       for (d = 0; d < dim; ++d) PetscRandomGetValueReal(rnd, &xref[d]);
800:       CoordinatesRefToReal(dim, dim, xi0, v0, J, xref, &x[p*dim]);

802:       weight[p] = 1.0;
803:       for (s = 0; s < Ns; ++s) species[p] = p % Ns;
804:     }
805:     PetscFree(pidx);
806:   }
807:   PetscRandomDestroy(&rnd);
808:   DMSwarmSortRestoreAccess(sw);
809:   DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **) &weight);
810:   DMSwarmRestoreField(sw, DMSwarmPICField_coor, NULL, NULL, (void **) &x);
811:   DMSwarmRestoreField(sw, "species", NULL, NULL, (void **) &species);
812:   DMSwarmMigrate(sw, removePoints);
813:   DMLocalizeCoordinates(sw);
814:   return 0;
815: }

817: /*@C
818:   DMSwarmInitializeVelocities - Set the initial velocities of particles using a distribution.

820:   Collective on dm

822:   Input Parameters:
823: + sw      - The DMSwarm object
824: . sampler - A function which uniformly samples the velocity PDF
825: - v0      - The velocity scale for nondimensionalization for each species

827:   Level: advanced

829: .seealso: DMSwarmComputeLocalSize(), DMSwarmInitializeCoordinates(), DMSwarmInitializeVelocitiesFromOptions()
830: @*/
831: PetscErrorCode DMSwarmInitializeVelocities(DM sw, PetscProbFunc sampler, const PetscReal v0[])
832: {
833:   PetscRandom  rnd;
834:   PetscReal   *v;
835:   PetscInt    *species;
836:   PetscInt     dim, Np, p;

838:   PetscRandomCreate(PetscObjectComm((PetscObject) sw), &rnd);
839:   PetscRandomSetInterval(rnd, 0, 1.);
840:   PetscRandomSetFromOptions(rnd);

842:   DMGetDimension(sw, &dim);
843:   DMSwarmGetLocalSize(sw, &Np);
844:   DMSwarmGetField(sw, "velocity", NULL, NULL, (void **) &v);
845:   DMSwarmGetField(sw, "species", NULL, NULL, (void **) &species);
846:   for (p = 0; p < Np; ++p) {
847:     PetscInt  s = species[p], d;
848:     PetscReal a[3], vel[3];

850:     for (d = 0; d < dim; ++d) PetscRandomGetValueReal(rnd, &a[d]);
851:     sampler(a, NULL, vel);
852:     for (d = 0; d < dim; ++d) {v[p*dim+d] = (v0[s] / v0[0]) * vel[d];}
853:   }
854:   DMSwarmRestoreField(sw, "velocity", NULL, NULL, (void **) &v);
855:   DMSwarmRestoreField(sw, "species", NULL, NULL, (void **) &species);
856:   PetscRandomDestroy(&rnd);
857:   return 0;
858: }

860: /*@
861:   DMSwarmInitializeVelocitiesFromOptions - Set the initial velocities of particles using a distribution determined from options.

863:   Collective on dm

865:   Input Parameters:
866: + sw      - The DMSwarm object
867: - v0      - The velocity scale for nondimensionalization for each species

869:   Level: advanced

871: .seealso: DMSwarmComputeLocalSize(), DMSwarmInitializeCoordinates(), DMSwarmInitializeVelocities()
872: @*/
873: PetscErrorCode DMSwarmInitializeVelocitiesFromOptions(DM sw, const PetscReal v0[])
874: {
875:   PetscProbFunc  sampler;
876:   PetscInt       dim;
877:   const char    *prefix;

879:   DMGetDimension(sw, &dim);
880:   PetscObjectGetOptionsPrefix((PetscObject) sw, &prefix);
881:   PetscProbCreateFromOptions(dim, prefix, "-dm_swarm_velocity_density", NULL, NULL, &sampler);
882:   DMSwarmInitializeVelocities(sw, sampler, v0);
883:   return 0;
884: }