Actual source code: swarmpic.c
1: #define PETSCDM_DLL
2: #include <petsc/private/dmswarmimpl.h>
3: #include <petscsf.h>
4: #include <petscdmda.h>
5: #include <petscdmplex.h>
6: #include "../src/dm/impls/swarm/data_bucket.h"
8: /*
9: Error chceking macto to ensure the swarm type is correct and that a cell DM has been set
10: */
11: #define DMSWARMPICVALID(dm) \
12: { \
13: DM_Swarm *_swarm = (DM_Swarm*)(dm)->data; \
14: if (_swarm->swarm_type != DMSWARM_PIC) SETERRQ(PetscObjectComm((PetscObject)(dm)),PETSC_ERR_SUP,"Only valid for DMSwarm-PIC. You must call DMSwarmSetType(dm,DMSWARM_PIC)"); \
15: else \
16: if (!_swarm->dmcell) SETERRQ(PetscObjectComm((PetscObject)(dm)),PETSC_ERR_SUP,"Only valid for DMSwarmPIC if the cell DM is set. You must call DMSwarmSetCellDM(dm,celldm)"); \
17: }
19: /* Coordinate insertition/addition API */
20: /*@C
21: DMSwarmSetPointsUniformCoordinates - Set point coordinates in a DMSwarm on a regular (ijk) grid
23: Collective on dm
25: Input parameters:
26: + dm - the DMSwarm
27: . min - minimum coordinate values in the x, y, z directions (array of length dim)
28: . max - maximum coordinate values in the x, y, z directions (array of length dim)
29: . npoints - number of points in each spatial direction (array of length dim)
30: - mode - indicates whether to append points to the swarm (ADD_VALUES), or over-ride existing points (INSERT_VALUES)
32: Level: beginner
34: Notes:
35: When using mode = INSERT_VALUES, this method will reset the number of particles in the DMSwarm
36: to be npoints[0]*npoints[1] (2D) or npoints[0]*npoints[1]*npoints[2] (3D). When using mode = ADD_VALUES,
37: new points will be appended to any already existing in the DMSwarm
39: .seealso: DMSwarmSetType(), DMSwarmSetCellDM(), DMSwarmType
40: @*/
41: PETSC_EXTERN PetscErrorCode DMSwarmSetPointsUniformCoordinates(DM dm,PetscReal min[],PetscReal max[],PetscInt npoints[],InsertMode mode)
42: {
43: PetscErrorCode ierr;
44: PetscReal gmin[] = {PETSC_MAX_REAL ,PETSC_MAX_REAL, PETSC_MAX_REAL};
45: PetscReal gmax[] = {PETSC_MIN_REAL, PETSC_MIN_REAL, PETSC_MIN_REAL};
46: PetscInt i,j,k,N,bs,b,n_estimate,n_curr,n_new_est,p,n_found;
47: Vec coorlocal;
48: const PetscScalar *_coor;
49: DM celldm;
50: PetscReal dx[3];
51: PetscInt _npoints[] = { 0, 0, 1 };
52: Vec pos;
53: PetscScalar *_pos;
54: PetscReal *swarm_coor;
55: PetscInt *swarm_cellid;
56: PetscSF sfcell = NULL;
57: const PetscSFNode *LA_sfcell;
60: DMSWARMPICVALID(dm);
61: DMSwarmGetCellDM(dm,&celldm);
62: DMGetCoordinatesLocal(celldm,&coorlocal);
63: VecGetSize(coorlocal,&N);
64: VecGetBlockSize(coorlocal,&bs);
65: N = N / bs;
66: VecGetArrayRead(coorlocal,&_coor);
67: for (i=0; i<N; i++) {
68: for (b=0; b<bs; b++) {
69: gmin[b] = PetscMin(gmin[b],PetscRealPart(_coor[bs*i+b]));
70: gmax[b] = PetscMax(gmax[b],PetscRealPart(_coor[bs*i+b]));
71: }
72: }
73: VecRestoreArrayRead(coorlocal,&_coor);
75: for (b=0; b<bs; b++) {
76: if (npoints[b] > 1) {
77: dx[b] = (max[b] - min[b])/((PetscReal)(npoints[b]-1));
78: } else {
79: dx[b] = 0.0;
80: }
82: _npoints[b] = npoints[b];
83: }
85: /* determine number of points living in the bounding box */
86: n_estimate = 0;
87: for (k=0; k<_npoints[2]; k++) {
88: for (j=0; j<_npoints[1]; j++) {
89: for (i=0; i<_npoints[0]; i++) {
90: PetscReal xp[] = {0.0,0.0,0.0};
91: PetscInt ijk[3];
92: PetscBool point_inside = PETSC_TRUE;
94: ijk[0] = i;
95: ijk[1] = j;
96: ijk[2] = k;
97: for (b=0; b<bs; b++) {
98: xp[b] = min[b] + ijk[b] * dx[b];
99: }
100: for (b=0; b<bs; b++) {
101: if (xp[b] < gmin[b]) { point_inside = PETSC_FALSE; }
102: if (xp[b] > gmax[b]) { point_inside = PETSC_FALSE; }
103: }
104: if (point_inside) { n_estimate++; }
105: }
106: }
107: }
109: /* create candidate list */
110: VecCreate(PETSC_COMM_SELF,&pos);
111: VecSetSizes(pos,bs*n_estimate,PETSC_DECIDE);
112: VecSetBlockSize(pos,bs);
113: VecSetFromOptions(pos);
114: VecGetArray(pos,&_pos);
116: n_estimate = 0;
117: for (k=0; k<_npoints[2]; k++) {
118: for (j=0; j<_npoints[1]; j++) {
119: for (i=0; i<_npoints[0]; i++) {
120: PetscReal xp[] = {0.0,0.0,0.0};
121: PetscInt ijk[3];
122: PetscBool point_inside = PETSC_TRUE;
124: ijk[0] = i;
125: ijk[1] = j;
126: ijk[2] = k;
127: for (b=0; b<bs; b++) {
128: xp[b] = min[b] + ijk[b] * dx[b];
129: }
130: for (b=0; b<bs; b++) {
131: if (xp[b] < gmin[b]) { point_inside = PETSC_FALSE; }
132: if (xp[b] > gmax[b]) { point_inside = PETSC_FALSE; }
133: }
134: if (point_inside) {
135: for (b=0; b<bs; b++) {
136: _pos[bs*n_estimate+b] = xp[b];
137: }
138: n_estimate++;
139: }
140: }
141: }
142: }
143: VecRestoreArray(pos,&_pos);
145: /* locate points */
146: DMLocatePoints(celldm,pos,DM_POINTLOCATION_NONE,&sfcell);
148: PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell);
149: n_found = 0;
150: for (p=0; p<n_estimate; p++) {
151: if (LA_sfcell[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) {
152: n_found++;
153: }
154: }
156: /* adjust size */
157: if (mode == ADD_VALUES) {
158: DMSwarmGetLocalSize(dm,&n_curr);
159: n_new_est = n_curr + n_found;
160: DMSwarmSetLocalSizes(dm,n_new_est,-1);
161: }
162: if (mode == INSERT_VALUES) {
163: n_curr = 0;
164: n_new_est = n_found;
165: DMSwarmSetLocalSizes(dm,n_new_est,-1);
166: }
168: /* initialize new coords, cell owners, pid */
169: VecGetArrayRead(pos,&_coor);
170: DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
171: DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
172: n_found = 0;
173: for (p=0; p<n_estimate; p++) {
174: if (LA_sfcell[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) {
175: for (b=0; b<bs; b++) {
176: swarm_coor[bs*(n_curr + n_found) + b] = PetscRealPart(_coor[bs*p+b]);
177: }
178: swarm_cellid[n_curr + n_found] = LA_sfcell[p].index;
179: n_found++;
180: }
181: }
182: DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
183: DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
184: VecRestoreArrayRead(pos,&_coor);
186: PetscSFDestroy(&sfcell);
187: VecDestroy(&pos);
189: return(0);
190: }
192: /*@C
193: DMSwarmSetPointCoordinates - Set point coordinates in a DMSwarm from a user defined list
195: Collective on dm
197: Input parameters:
198: + dm - the DMSwarm
199: . npoints - the number of points to insert
200: . coor - the coordinate values
201: . redundant - if set to PETSC_TRUE, it is assumed that npoints and coor[] are only valid on rank 0 and should be broadcast to other ranks
202: - mode - indicates whether to append points to the swarm (ADD_VALUES), or over-ride existing points (INSERT_VALUES)
204: Level: beginner
206: Notes:
207: If the user has specified redundant = PETSC_FALSE, the cell DM will attempt to locate the coordinates provided by coor[] within
208: its sub-domain. If they any values within coor[] are not located in the sub-domain, they will be ignored and will not get
209: added to the DMSwarm.
211: .seealso: DMSwarmSetType(), DMSwarmSetCellDM(), DMSwarmType, DMSwarmSetPointsUniformCoordinates()
212: @*/
213: PETSC_EXTERN PetscErrorCode DMSwarmSetPointCoordinates(DM dm,PetscInt npoints,PetscReal coor[],PetscBool redundant,InsertMode mode)
214: {
215: PetscErrorCode ierr;
216: PetscReal gmin[] = {PETSC_MAX_REAL ,PETSC_MAX_REAL, PETSC_MAX_REAL};
217: PetscReal gmax[] = {PETSC_MIN_REAL, PETSC_MIN_REAL, PETSC_MIN_REAL};
218: PetscInt i,N,bs,b,n_estimate,n_curr,n_new_est,p,n_found;
219: Vec coorlocal;
220: const PetscScalar *_coor;
221: DM celldm;
222: Vec pos;
223: PetscScalar *_pos;
224: PetscReal *swarm_coor;
225: PetscInt *swarm_cellid;
226: PetscSF sfcell = NULL;
227: const PetscSFNode *LA_sfcell;
228: PetscReal *my_coor;
229: PetscInt my_npoints;
230: PetscMPIInt rank;
231: MPI_Comm comm;
234: DMSWARMPICVALID(dm);
235: PetscObjectGetComm((PetscObject)dm,&comm);
236: MPI_Comm_rank(comm,&rank);
238: DMSwarmGetCellDM(dm,&celldm);
239: DMGetCoordinatesLocal(celldm,&coorlocal);
240: VecGetSize(coorlocal,&N);
241: VecGetBlockSize(coorlocal,&bs);
242: N = N / bs;
243: VecGetArrayRead(coorlocal,&_coor);
244: for (i=0; i<N; i++) {
245: for (b=0; b<bs; b++) {
246: gmin[b] = PetscMin(gmin[b],PetscRealPart(_coor[bs*i+b]));
247: gmax[b] = PetscMax(gmax[b],PetscRealPart(_coor[bs*i+b]));
248: }
249: }
250: VecRestoreArrayRead(coorlocal,&_coor);
252: /* broadcast points from rank 0 if requested */
253: if (redundant) {
254: my_npoints = npoints;
255: MPI_Bcast(&my_npoints,1,MPIU_INT,0,comm);
257: if (rank > 0) { /* allocate space */
258: PetscMalloc1(bs*my_npoints,&my_coor);
259: } else {
260: my_coor = coor;
261: }
262: MPI_Bcast(my_coor,bs*my_npoints,MPIU_REAL,0,comm);
263: } else {
264: my_npoints = npoints;
265: my_coor = coor;
266: }
268: /* determine the number of points living in the bounding box */
269: n_estimate = 0;
270: for (i=0; i<my_npoints; i++) {
271: PetscBool point_inside = PETSC_TRUE;
273: for (b=0; b<bs; b++) {
274: if (my_coor[bs*i+b] < gmin[b]) { point_inside = PETSC_FALSE; }
275: if (my_coor[bs*i+b] > gmax[b]) { point_inside = PETSC_FALSE; }
276: }
277: if (point_inside) { n_estimate++; }
278: }
280: /* create candidate list */
281: VecCreate(PETSC_COMM_SELF,&pos);
282: VecSetSizes(pos,bs*n_estimate,PETSC_DECIDE);
283: VecSetBlockSize(pos,bs);
284: VecSetFromOptions(pos);
285: VecGetArray(pos,&_pos);
287: n_estimate = 0;
288: for (i=0; i<my_npoints; i++) {
289: PetscBool point_inside = PETSC_TRUE;
291: for (b=0; b<bs; b++) {
292: if (my_coor[bs*i+b] < gmin[b]) { point_inside = PETSC_FALSE; }
293: if (my_coor[bs*i+b] > gmax[b]) { point_inside = PETSC_FALSE; }
294: }
295: if (point_inside) {
296: for (b=0; b<bs; b++) {
297: _pos[bs*n_estimate+b] = my_coor[bs*i+b];
298: }
299: n_estimate++;
300: }
301: }
302: VecRestoreArray(pos,&_pos);
304: /* locate points */
305: DMLocatePoints(celldm,pos,DM_POINTLOCATION_NONE,&sfcell);
307: PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell);
308: n_found = 0;
309: for (p=0; p<n_estimate; p++) {
310: if (LA_sfcell[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) {
311: n_found++;
312: }
313: }
315: /* adjust size */
316: if (mode == ADD_VALUES) {
317: DMSwarmGetLocalSize(dm,&n_curr);
318: n_new_est = n_curr + n_found;
319: DMSwarmSetLocalSizes(dm,n_new_est,-1);
320: }
321: if (mode == INSERT_VALUES) {
322: n_curr = 0;
323: n_new_est = n_found;
324: DMSwarmSetLocalSizes(dm,n_new_est,-1);
325: }
327: /* initialize new coords, cell owners, pid */
328: VecGetArrayRead(pos,&_coor);
329: DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
330: DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
331: n_found = 0;
332: for (p=0; p<n_estimate; p++) {
333: if (LA_sfcell[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) {
334: for (b=0; b<bs; b++) {
335: swarm_coor[bs*(n_curr + n_found) + b] = PetscRealPart(_coor[bs*p+b]);
336: }
337: swarm_cellid[n_curr + n_found] = LA_sfcell[p].index;
338: n_found++;
339: }
340: }
341: DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
342: DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
343: VecRestoreArrayRead(pos,&_coor);
345: if (redundant) {
346: if (rank > 0) {
347: PetscFree(my_coor);
348: }
349: }
350: PetscSFDestroy(&sfcell);
351: VecDestroy(&pos);
353: return(0);
354: }
356: extern PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_DA(DM,DM,DMSwarmPICLayoutType,PetscInt);
357: extern PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX(DM,DM,DMSwarmPICLayoutType,PetscInt);
359: /*@C
360: DMSwarmInsertPointsUsingCellDM - Insert point coordinates within each cell
362: Not collective
364: Input parameters:
365: + dm - the DMSwarm
366: . layout_type - method used to fill each cell with the cell DM
367: - fill_param - parameter controlling how many points per cell are added (the meaning of this parameter is dependent on the layout type)
369: Level: beginner
371: Notes:
373: The insert method will reset any previous defined points within the DMSwarm.
375: When using a DMDA both 2D and 3D are supported for all layout types provided you are using DMDA_ELEMENT_Q1.
377: When using a DMPLEX the following case are supported:
378: (i) DMSWARMPIC_LAYOUT_REGULAR: 2D (triangle),
379: (ii) DMSWARMPIC_LAYOUT_GAUSS: 2D and 3D provided the cell is a tri/tet or a quad/hex,
380: (iii) DMSWARMPIC_LAYOUT_SUBDIVISION: 2D and 3D for quad/hex and 2D tri.
382: .seealso: DMSwarmPICLayoutType, DMSwarmSetType(), DMSwarmSetCellDM(), DMSwarmType
383: @*/
384: PETSC_EXTERN PetscErrorCode DMSwarmInsertPointsUsingCellDM(DM dm,DMSwarmPICLayoutType layout_type,PetscInt fill_param)
385: {
387: DM celldm;
388: PetscBool isDA,isPLEX;
391: DMSWARMPICVALID(dm);
392: DMSwarmGetCellDM(dm,&celldm);
393: PetscObjectTypeCompare((PetscObject)celldm,DMDA,&isDA);
394: PetscObjectTypeCompare((PetscObject)celldm,DMPLEX,&isPLEX);
395: if (isDA) {
396: private_DMSwarmInsertPointsUsingCellDM_DA(dm,celldm,layout_type,fill_param);
397: } else if (isPLEX) {
398: private_DMSwarmInsertPointsUsingCellDM_PLEX(dm,celldm,layout_type,fill_param);
399: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only supported for cell DMs of type DMDA and DMPLEX");
401: return(0);
402: }
404: extern PetscErrorCode private_DMSwarmSetPointCoordinatesCellwise_PLEX(DM,DM,PetscInt,PetscReal*);
406: /*@C
407: DMSwarmSetPointCoordinatesCellwise - Insert point coordinates (defined over the reference cell) within each cell
409: Not collective
411: Input parameters:
412: + dm - the DMSwarm
413: . celldm - the cell DM
414: . npoints - the number of points to insert in each cell
415: - xi - the coordinates (defined in the local coordinate system for each cell) to insert
417: Level: beginner
419: Notes:
420: The method will reset any previous defined points within the DMSwarm.
421: Only supported for DMPLEX. If you are using a DMDA it is recommended to either use
422: DMSwarmInsertPointsUsingCellDM(), or extract and set the coordinates yourself the following code
424: $ PetscReal *coor;
425: $ DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&coor);
426: $ // user code to define the coordinates here
427: $ DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&coor);
429: .seealso: DMSwarmSetCellDM(), DMSwarmInsertPointsUsingCellDM()
430: @*/
431: PETSC_EXTERN PetscErrorCode DMSwarmSetPointCoordinatesCellwise(DM dm,PetscInt npoints,PetscReal xi[])
432: {
434: DM celldm;
435: PetscBool isDA,isPLEX;
438: DMSWARMPICVALID(dm);
439: DMSwarmGetCellDM(dm,&celldm);
440: PetscObjectTypeCompare((PetscObject)celldm,DMDA,&isDA);
441: PetscObjectTypeCompare((PetscObject)celldm,DMPLEX,&isPLEX);
442: if (isDA) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only supported for cell DMs of type DMPLEX. Recommended you use DMSwarmInsertPointsUsingCellDM()");
443: else if (isPLEX) {
444: private_DMSwarmSetPointCoordinatesCellwise_PLEX(dm,celldm,npoints,xi);
445: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only supported for cell DMs of type DMDA and DMPLEX");
447: return(0);
448: }
450: /* Field projection API */
451: extern PetscErrorCode private_DMSwarmProjectFields_DA(DM swarm,DM celldm,PetscInt project_type,PetscInt nfields,DMSwarmDataField dfield[],Vec vecs[]);
452: extern PetscErrorCode private_DMSwarmProjectFields_PLEX(DM swarm,DM celldm,PetscInt project_type,PetscInt nfields,DMSwarmDataField dfield[],Vec vecs[]);
454: /*@C
455: DMSwarmProjectFields - Project a set of swarm fields onto the cell DM
457: Collective on dm
459: Input parameters:
460: + dm - the DMSwarm
461: . nfields - the number of swarm fields to project
462: . fieldnames - the textual names of the swarm fields to project
463: . fields - an array of Vec's of length nfields
464: - reuse - flag indicating whether the array and contents of fields should be re-used or internally allocated
466: Currently, the only available projection method consists of
467: phi_i = \sum_{p=0}^{np} N_i(x_p) phi_p dJ / \sum_{p=0}^{np} N_i(x_p) dJ
468: where phi_p is the swarm field at point p,
469: N_i() is the cell DM basis function at vertex i,
470: dJ is the determinant of the cell Jacobian and
471: phi_i is the projected vertex value of the field phi.
473: Level: beginner
475: Notes:
477: If reuse = PETSC_FALSE, this function will allocate the array of Vec's, and each individual Vec.
478: The user is responsible for destroying both the array and the individual Vec objects.
480: Only swarm fields registered with data type = PETSC_REAL can be projected onto the cell DM.
482: Only swarm fields of block size = 1 can currently be projected.
484: The only projection methods currently only support the DA (2D) and PLEX (triangles 2D).
486: .seealso: DMSwarmSetType(), DMSwarmSetCellDM(), DMSwarmType
487: @*/
488: PETSC_EXTERN PetscErrorCode DMSwarmProjectFields(DM dm,PetscInt nfields,const char *fieldnames[],Vec **fields,PetscBool reuse)
489: {
490: DM_Swarm *swarm = (DM_Swarm*)dm->data;
491: DMSwarmDataField *gfield;
492: DM celldm;
493: PetscBool isDA,isPLEX;
494: Vec *vecs;
495: PetscInt f,nvecs;
496: PetscInt project_type = 0;
500: DMSWARMPICVALID(dm);
501: DMSwarmGetCellDM(dm,&celldm);
502: PetscMalloc1(nfields,&gfield);
503: nvecs = 0;
504: for (f=0; f<nfields; f++) {
505: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,fieldnames[f],&gfield[f]);
506: 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");
507: if (gfield[f]->bs != 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Projection only valid for fields with block size = 1");
508: nvecs += gfield[f]->bs;
509: }
510: if (!reuse) {
511: PetscMalloc1(nvecs,&vecs);
512: for (f=0; f<nvecs; f++) {
513: DMCreateGlobalVector(celldm,&vecs[f]);
514: PetscObjectSetName((PetscObject)vecs[f],gfield[f]->name);
515: }
516: } else {
517: vecs = *fields;
518: }
520: PetscObjectTypeCompare((PetscObject)celldm,DMDA,&isDA);
521: PetscObjectTypeCompare((PetscObject)celldm,DMPLEX,&isPLEX);
522: if (isDA) {
523: private_DMSwarmProjectFields_DA(dm,celldm,project_type,nfields,gfield,vecs);
524: } else if (isPLEX) {
525: private_DMSwarmProjectFields_PLEX(dm,celldm,project_type,nfields,gfield,vecs);
526: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only supported for cell DMs of type DMDA and DMPLEX");
528: PetscFree(gfield);
529: if (!reuse) {
530: *fields = vecs;
531: }
533: return(0);
534: }
536: /*@C
537: DMSwarmCreatePointPerCellCount - Count the number of points within all cells in the cell DM
539: Not collective
541: Input parameter:
542: . dm - the DMSwarm
544: Output parameters:
545: + ncells - the number of cells in the cell DM (optional argument, pass NULL to ignore)
546: - count - array of length ncells containing the number of points per cell
548: Level: beginner
550: Notes:
551: The array count is allocated internally and must be free'd by the user.
553: .seealso: DMSwarmSetType(), DMSwarmSetCellDM(), DMSwarmType
554: @*/
555: PETSC_EXTERN PetscErrorCode DMSwarmCreatePointPerCellCount(DM dm,PetscInt *ncells,PetscInt **count)
556: {
558: PetscBool isvalid;
559: PetscInt nel;
560: PetscInt *sum;
563: DMSwarmSortGetIsValid(dm,&isvalid);
564: nel = 0;
565: if (isvalid) {
566: PetscInt e;
568: DMSwarmSortGetSizes(dm,&nel,NULL);
570: PetscMalloc1(nel,&sum);
571: for (e=0; e<nel; e++) {
572: DMSwarmSortGetNumberOfPointsPerCell(dm,e,&sum[e]);
573: }
574: } else {
575: DM celldm;
576: PetscBool isda,isplex,isshell;
577: PetscInt p,npoints;
578: PetscInt *swarm_cellid;
580: /* get the number of cells */
581: DMSwarmGetCellDM(dm,&celldm);
582: PetscObjectTypeCompare((PetscObject)celldm,DMDA,&isda);
583: PetscObjectTypeCompare((PetscObject)celldm,DMPLEX,&isplex);
584: PetscObjectTypeCompare((PetscObject)celldm,DMSHELL,&isshell);
585: if (isda) {
586: PetscInt _nel,_npe;
587: const PetscInt *_element;
589: DMDAGetElements(celldm,&_nel,&_npe,&_element);
590: nel = _nel;
591: DMDARestoreElements(celldm,&_nel,&_npe,&_element);
592: } else if (isplex) {
593: PetscInt ps,pe;
595: DMPlexGetHeightStratum(celldm,0,&ps,&pe);
596: nel = pe - ps;
597: } else if (isshell) {
598: PetscErrorCode (*method_DMShellGetNumberOfCells)(DM,PetscInt*);
600: PetscObjectQueryFunction((PetscObject)celldm,"DMGetNumberOfCells_C",&method_DMShellGetNumberOfCells);
601: if (method_DMShellGetNumberOfCells) {
602: method_DMShellGetNumberOfCells(celldm,&nel);
603: } 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);");
604: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Cannot determine the number of cells for a DM not of type DA, PLEX or SHELL");
606: PetscMalloc1(nel,&sum);
607: PetscArrayzero(sum,nel);
608: DMSwarmGetLocalSize(dm,&npoints);
609: DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
610: for (p=0; p<npoints; p++) {
611: if (swarm_cellid[p] != DMLOCATEPOINT_POINT_NOT_FOUND) {
612: sum[ swarm_cellid[p] ]++;
613: }
614: }
615: DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
616: }
617: if (ncells) { *ncells = nel; }
618: *count = sum;
619: return(0);
620: }