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