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