Actual source code: swarmpic.c
petsc-3.14.6 2021-03-30
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: }