xref: /petsc/src/dm/impls/swarm/swarmpic.c (revision 71844bb8d9b6b640b934b4f3526f4bb2bc397e8c)
1 
2 #define PETSCDM_DLL
3 #include <petsc/private/dmswarmimpl.h>    /*I   "petscdmswarm.h"   I*/
4 #include <petscsf.h>
5 #include <petscdmda.h>
6 #include <petscdmplex.h>
7 
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 }
18 
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 {
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   Vec pos;
52   PetscScalar *_pos;
53   PetscReal *swarm_coor;
54   PetscInt *swarm_cellid;
55   PetscSF sfcell = NULL;
56   const PetscSFNode *LA_sfcell;
57 
58   PetscFunctionBegin;
59   DMSWARMPICVALID(dm);
60   ierr = DMSwarmGetCellDM(dm,&celldm);CHKERRQ(ierr);
61   ierr = DMGetCoordinatesLocal(celldm,&coorlocal);CHKERRQ(ierr);
62   ierr = VecGetSize(coorlocal,&N);CHKERRQ(ierr);
63   ierr = VecGetBlockSize(coorlocal,&bs);CHKERRQ(ierr);
64   N = N / bs;
65   ierr = VecGetArrayRead(coorlocal,&_coor);CHKERRQ(ierr);
66   for (i=0; i<N; i++) {
67     for (b=0; b<bs; b++) {
68       gmin[b] = PetscMin(gmin[b],_coor[bs*i+b]);
69       gmax[b] = PetscMax(gmax[b],_coor[bs*i+b]);
70     }
71   }
72   ierr = VecRestoreArrayRead(coorlocal,&_coor);CHKERRQ(ierr);
73 
74   for (b=0; b<bs; b++) {
75     if (npoints[b] > 1) {
76       dx[b] = (max[b] - min[b])/((PetscReal)(npoints[b]-1));
77     } else {
78       dx[b] = 0.0;
79     }
80   }
81 
82   /* determine number of points living in the bounding box */
83   n_estimate = 0;
84   if (bs == 2) { npoints[2] = 1; }
85   for (k=0; k<npoints[2]; k++) {
86     for (j=0; j<npoints[1]; j++) {
87       for (i=0; i<npoints[0]; i++) {
88         PetscReal xp[] = {0.0,0.0,0.0};
89         PetscInt ijk[3];
90         PetscBool point_inside = PETSC_TRUE;
91 
92         ijk[0] = i;
93         ijk[1] = j;
94         ijk[2] = k;
95         for (b=0; b<bs; b++) {
96           xp[b] = min[b] + ijk[b] * dx[b];
97         }
98         for (b=0; b<bs; b++) {
99           if (xp[b] < gmin[b]) { point_inside = PETSC_FALSE; }
100           if (xp[b] > gmax[b]) { point_inside = PETSC_FALSE; }
101         }
102         if (point_inside) { n_estimate++; }
103       }
104     }
105   }
106 
107   /* create candidate list */
108   ierr = VecCreate(PetscObjectComm((PetscObject)dm),&pos);CHKERRQ(ierr);
109   ierr = VecSetSizes(pos,bs*n_estimate,PETSC_DECIDE);CHKERRQ(ierr);
110   ierr = VecSetBlockSize(pos,bs);CHKERRQ(ierr);
111   ierr = VecSetFromOptions(pos);CHKERRQ(ierr);
112   ierr = VecGetArray(pos,&_pos);CHKERRQ(ierr);
113 
114   n_estimate = 0;
115   for (k=0; k<npoints[2]; k++) {
116     for (j=0; j<npoints[1]; j++) {
117       for (i=0; i<npoints[0]; i++) {
118         PetscReal xp[] = {0.0,0.0,0.0};
119         PetscInt ijk[3];
120         PetscBool point_inside = PETSC_TRUE;
121 
122         ijk[0] = i;
123         ijk[1] = j;
124         ijk[2] = k;
125         for (b=0; b<bs; b++) {
126           xp[b] = min[b] + ijk[b] * dx[b];
127         }
128         for (b=0; b<bs; b++) {
129           if (xp[b] < gmin[b]) { point_inside = PETSC_FALSE; }
130           if (xp[b] > gmax[b]) { point_inside = PETSC_FALSE; }
131         }
132         if (point_inside) {
133           for (b=0; b<bs; b++) {
134             _pos[bs*n_estimate+b] = xp[b];
135           }
136           n_estimate++;
137         }
138       }
139     }
140   }
141   ierr = VecRestoreArray(pos,&_pos);CHKERRQ(ierr);
142 
143   /* locate points */
144   ierr = DMLocatePoints(celldm,pos,DM_POINTLOCATION_NONE,&sfcell);CHKERRQ(ierr);
145 
146   ierr = PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell);CHKERRQ(ierr);
147   n_found = 0;
148   for (p=0; p<n_estimate; p++) {
149     if (LA_sfcell[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) {
150       n_found++;
151     }
152   }
153 
154   /* adjust size */
155   if (mode == ADD_VALUES) {
156     ierr = DMSwarmGetLocalSize(dm,&n_curr);CHKERRQ(ierr);
157     n_new_est = n_curr + n_found;
158     ierr = DMSwarmSetLocalSizes(dm,n_new_est,-1);CHKERRQ(ierr);
159   }
160   if (mode == INSERT_VALUES) {
161     n_curr = 0;
162     n_new_est = n_found;
163     ierr = DMSwarmSetLocalSizes(dm,n_new_est,-1);CHKERRQ(ierr);
164   }
165 
166   /* initialize new coords, cell owners, pid */
167   ierr = VecGetArrayRead(pos,&_coor);CHKERRQ(ierr);
168   ierr = DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr);
169   ierr = DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr);
170   n_found = 0;
171   for (p=0; p<n_estimate; p++) {
172     if (LA_sfcell[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) {
173       for (b=0; b<bs; b++) {
174         swarm_coor[bs*(n_curr + n_found) + b] = _coor[bs*p+b];
175       }
176       swarm_cellid[n_curr + n_found] = LA_sfcell[p].index;
177       n_found++;
178     }
179   }
180   ierr = DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr);
181   ierr = DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr);
182   ierr = VecRestoreArrayRead(pos,&_coor);CHKERRQ(ierr);
183 
184   ierr = PetscSFDestroy(&sfcell);CHKERRQ(ierr);
185   ierr = VecDestroy(&pos);CHKERRQ(ierr);
186 
187   PetscFunctionReturn(0);
188 }
189 
190 /*@C
191    DMSwarmSetPointCoordinates - Set point coordinates in a DMSwarm from a user defined list
192 
193    Collective on DM
194 
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)
201 
202    Level: beginner
203 
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.
208 
209 .seealso: DMSwarmSetType(), DMSwarmSetCellDM(), DMSwarmType, DMSwarmSetPointsUniformCoordinates()
210 @*/
211 PETSC_EXTERN PetscErrorCode DMSwarmSetPointCoordinates(DM dm,PetscInt npoints,PetscReal coor[],PetscBool redundant,InsertMode mode)
212 {
213   PetscErrorCode ierr;
214   PetscReal gmin[] = {PETSC_MAX_REAL ,PETSC_MAX_REAL, PETSC_MAX_REAL};
215   PetscReal gmax[] = {PETSC_MIN_REAL, PETSC_MIN_REAL, PETSC_MIN_REAL};
216   PetscInt i,N,bs,b,n_estimate,n_curr,n_new_est,p,n_found;
217   Vec coorlocal;
218   const PetscScalar *_coor;
219   DM celldm;
220   Vec pos;
221   PetscScalar *_pos;
222   PetscReal *swarm_coor;
223   PetscInt *swarm_cellid;
224   PetscSF sfcell = NULL;
225   const PetscSFNode *LA_sfcell;
226   PetscReal *my_coor;
227   PetscInt my_npoints;
228   PetscMPIInt rank;
229   MPI_Comm comm;
230 
231   PetscFunctionBegin;
232   DMSWARMPICVALID(dm);
233   ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
234   ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
235 
236   ierr = DMSwarmGetCellDM(dm,&celldm);CHKERRQ(ierr);
237   ierr = DMGetCoordinatesLocal(celldm,&coorlocal);CHKERRQ(ierr);
238   ierr = VecGetSize(coorlocal,&N);CHKERRQ(ierr);
239   ierr = VecGetBlockSize(coorlocal,&bs);CHKERRQ(ierr);
240   N = N / bs;
241   ierr = VecGetArrayRead(coorlocal,&_coor);CHKERRQ(ierr);
242   for (i=0; i<N; i++) {
243     for (b=0; b<bs; b++) {
244       gmin[b] = PetscMin(gmin[b],_coor[bs*i+b]);
245       gmax[b] = PetscMax(gmax[b],_coor[bs*i+b]);
246     }
247   }
248   ierr = VecRestoreArrayRead(coorlocal,&_coor);CHKERRQ(ierr);
249 
250   /* broadcast points from rank 0 if requested */
251   if (redundant) {
252     my_npoints = npoints;
253     ierr = MPI_Bcast(&my_npoints,1,MPIU_INT,0,comm);CHKERRQ(ierr);
254 
255     if (rank > 0) { /* allocate space */
256       ierr = PetscMalloc1(bs*my_npoints,&my_coor);CHKERRQ(ierr);
257     } else {
258       my_coor = coor;
259     }
260     ierr = MPI_Bcast(my_coor,bs*my_npoints,MPIU_REAL,0,comm);CHKERRQ(ierr);
261   } else {
262     my_npoints = npoints;
263     my_coor = coor;
264   }
265 
266   /* determine the number of points living in the bounding box */
267   n_estimate = 0;
268   for (i=0; i<my_npoints; i++) {
269     PetscBool point_inside = PETSC_TRUE;
270 
271     for (b=0; b<bs; b++) {
272       if (my_coor[bs*i+b] < gmin[b]) { point_inside = PETSC_FALSE; }
273       if (my_coor[bs*i+b] > gmax[b]) { point_inside = PETSC_FALSE; }
274     }
275     if (point_inside) { n_estimate++; }
276   }
277 
278   /* create candidate list */
279   ierr = VecCreate(PetscObjectComm((PetscObject)dm),&pos);CHKERRQ(ierr);
280   ierr = VecSetSizes(pos,bs*n_estimate,PETSC_DECIDE);CHKERRQ(ierr);
281   ierr = VecSetBlockSize(pos,bs);CHKERRQ(ierr);
282   ierr = VecSetFromOptions(pos);CHKERRQ(ierr);
283   ierr = VecGetArray(pos,&_pos);CHKERRQ(ierr);
284 
285   n_estimate = 0;
286   for (i=0; i<my_npoints; i++) {
287     PetscBool point_inside = PETSC_TRUE;
288 
289     for (b=0; b<bs; b++) {
290       if (my_coor[bs*i+b] < gmin[b]) { point_inside = PETSC_FALSE; }
291       if (my_coor[bs*i+b] > gmax[b]) { point_inside = PETSC_FALSE; }
292     }
293     if (point_inside) {
294       for (b=0; b<bs; b++) {
295         _pos[bs*n_estimate+b] = my_coor[bs*i+b];
296       }
297       n_estimate++;
298     }
299   }
300   ierr = VecRestoreArray(pos,&_pos);CHKERRQ(ierr);
301 
302   /* locate points */
303   ierr = DMLocatePoints(celldm,pos,DM_POINTLOCATION_NONE,&sfcell);CHKERRQ(ierr);
304 
305   ierr = PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell);CHKERRQ(ierr);
306   n_found = 0;
307   for (p=0; p<n_estimate; p++) {
308     if (LA_sfcell[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) {
309       n_found++;
310     }
311   }
312 
313   /* adjust size */
314   if (mode == ADD_VALUES) {
315     ierr = DMSwarmGetLocalSize(dm,&n_curr);CHKERRQ(ierr);
316     n_new_est = n_curr + n_found;
317     ierr = DMSwarmSetLocalSizes(dm,n_new_est,-1);CHKERRQ(ierr);
318   }
319   if (mode == INSERT_VALUES) {
320     n_curr = 0;
321     n_new_est = n_found;
322     ierr = DMSwarmSetLocalSizes(dm,n_new_est,-1);CHKERRQ(ierr);
323   }
324 
325   /* initialize new coords, cell owners, pid */
326   ierr = VecGetArrayRead(pos,&_coor);CHKERRQ(ierr);
327   ierr = DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr);
328   ierr = DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr);
329   n_found = 0;
330   for (p=0; p<n_estimate; p++) {
331     if (LA_sfcell[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) {
332       for (b=0; b<bs; b++) {
333         swarm_coor[bs*(n_curr + n_found) + b] = _coor[bs*p+b];
334       }
335       swarm_cellid[n_curr + n_found] = LA_sfcell[p].index;
336       n_found++;
337     }
338   }
339   ierr = DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr);
340   ierr = DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr);
341   ierr = VecRestoreArrayRead(pos,&_coor);CHKERRQ(ierr);
342 
343   if (redundant) {
344     if (rank > 0) {
345       ierr = PetscFree(my_coor);CHKERRQ(ierr);
346     }
347   }
348   ierr = PetscSFDestroy(&sfcell);CHKERRQ(ierr);
349   ierr = VecDestroy(&pos);CHKERRQ(ierr);
350 
351   PetscFunctionReturn(0);
352 }
353 
354 extern PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_DA(DM,DM,DMSwarmPICLayoutType,PetscInt);
355 extern PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX(DM,DM,DMSwarmPICLayoutType,PetscInt);
356 
357 /*@C
358    DMSwarmInsertPointsUsingCellDM - Insert point coordinates within each cell
359 
360    Not collective
361 
362    Input parameters:
363 +  dm - the DMSwarm
364 .  layout_type - method used to fill each cell with the cell DM
365 -  fill_param - parameter controlling how many points per cell are added (the meaning of this parameter is dependent on the layout type)
366 
367  Level: beginner
368 
369  Notes:
370  The insert method will reset any previous defined points within the DMSwarm
371 
372 .seealso: DMSwarmPICLayoutType, DMSwarmSetType(), DMSwarmSetCellDM(), DMSwarmType
373 @*/
374 PETSC_EXTERN PetscErrorCode DMSwarmInsertPointsUsingCellDM(DM dm,DMSwarmPICLayoutType layout_type,PetscInt fill_param)
375 {
376   PetscErrorCode ierr;
377   DM celldm;
378   PetscBool isDA,isPLEX;
379 
380   PetscFunctionBegin;
381   DMSWARMPICVALID(dm);
382   ierr = DMSwarmGetCellDM(dm,&celldm);CHKERRQ(ierr);
383   ierr = PetscObjectTypeCompare((PetscObject)celldm,DMDA,&isDA);CHKERRQ(ierr);
384   ierr = PetscObjectTypeCompare((PetscObject)celldm,DMPLEX,&isPLEX);CHKERRQ(ierr);
385   if (isDA) {
386     ierr = private_DMSwarmInsertPointsUsingCellDM_DA(dm,celldm,layout_type,fill_param);CHKERRQ(ierr);
387   } else if (isPLEX) {
388     ierr = private_DMSwarmInsertPointsUsingCellDM_PLEX(dm,celldm,layout_type,fill_param);CHKERRQ(ierr);
389   } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only supported for cell DMs of type DMDA and DMPLEX");
390 
391   PetscFunctionReturn(0);
392 }
393 
394 /*
395 PETSC_EXTERN PetscErrorCode DMSwarmAddPointCoordinatesCellWise(DM dm,PetscInt cell,PetscInt npoints,PetscReal xi[],PetscBool proximity_initialization)
396 {
397   PetscFunctionBegin;
398   PetscFunctionReturn(0);
399 }
400 */
401 
402 /* Field projection API */
403 /*
404 PETSC_EXTERN PetscErrorCode DMSwarmProjectFields(DM dm,PetscInt project_type,PetscInt nfields,const char *fieldnames[],Vec *fields)
405 {
406   PetscFunctionBegin;
407   PetscFunctionReturn(0);
408 }
409 */
410 
411 /*@C
412    DMSwarmCreatePointPerCellCount - Count the number of points within all cells in the cell DM
413 
414    Not collective
415 
416    Input parameter:
417 .  dm - the DMSwarm
418 
419    Output parameters:
420 +  ncells - the number of cells in the cell DM (optional argument, pass NULL to ignore)
421 -  count - array of length ncells containing the number of points per cell
422 
423    Level: beginner
424 
425    Notes:
426    The array count is allocated internally and must be free'd by the user.
427 
428 .seealso: DMSwarmSetType(), DMSwarmSetCellDM(), DMSwarmType
429 @*/
430 PETSC_EXTERN PetscErrorCode DMSwarmCreatePointPerCellCount(DM dm,PetscInt *ncells,PetscInt **count)
431 {
432   PetscErrorCode ierr;
433   PetscBool      isvalid;
434   PetscInt       nel;
435   PetscInt       *sum;
436 
437   PetscFunctionBegin;
438   ierr = DMSwarmSortGetIsValid(dm,&isvalid);CHKERRQ(ierr);
439   nel = 0;
440   if (isvalid) {
441     PetscInt e;
442 
443     ierr = DMSwarmSortGetSizes(dm,&nel,NULL);CHKERRQ(ierr);
444 
445     ierr = PetscMalloc1(nel,&sum);CHKERRQ(ierr);
446     for (e=0; e<nel; e++) {
447       ierr = DMSwarmSortGetNumberOfPointsPerCell(dm,e,&sum[e]);CHKERRQ(ierr);
448     }
449   } else {
450     DM        celldm;
451     PetscBool isda,isplex,isshell;
452     PetscInt  p,npoints;
453     PetscInt *swarm_cellid;
454 
455     /* get the number of cells */
456     ierr = DMSwarmGetCellDM(dm,&celldm);CHKERRQ(ierr);
457     ierr = PetscObjectTypeCompare((PetscObject)celldm,DMDA,&isda);CHKERRQ(ierr);
458     ierr = PetscObjectTypeCompare((PetscObject)celldm,DMPLEX,&isplex);CHKERRQ(ierr);
459     ierr = PetscObjectTypeCompare((PetscObject)celldm,DMSHELL,&isshell);CHKERRQ(ierr);
460     if (isda) {
461       PetscInt _nel,_npe;
462       const PetscInt *_element;
463 
464       ierr = DMDAGetElements(celldm,&_nel,&_npe,&_element);CHKERRQ(ierr);
465       nel = _nel;
466       ierr = DMDARestoreElements(celldm,&_nel,&_npe,&_element);CHKERRQ(ierr);
467     } else if (isplex) {
468       PetscInt ps,pe;
469 
470       ierr = DMPlexGetHeightStratum(celldm,0,&ps,&pe);CHKERRQ(ierr);
471       nel = pe - ps;
472     } else if (isshell) {
473       PetscErrorCode (*method_DMShellGetNumberOfCells)(DM,PetscInt*);
474 
475       ierr = PetscObjectQueryFunction((PetscObject)celldm,"DMGetNumberOfCells_C",&method_DMShellGetNumberOfCells);CHKERRQ(ierr);
476       if (method_DMShellGetNumberOfCells) {
477         ierr = method_DMShellGetNumberOfCells(celldm,&nel);CHKERRQ(ierr);
478       } 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 );");
479     } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Cannot determine the number of cells for a DM not of type DA, PLEX or SHELL");
480 
481     ierr = PetscMalloc1(nel,&sum);CHKERRQ(ierr);
482     ierr = PetscMemzero(sum,sizeof(PetscInt)*nel);CHKERRQ(ierr);
483     ierr = DMSwarmGetLocalSize(dm,&npoints);CHKERRQ(ierr);
484     ierr = DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr);
485     for (p=0; p<npoints; p++) {
486       if (swarm_cellid[p] != DMLOCATEPOINT_POINT_NOT_FOUND) {
487         sum[ swarm_cellid[p] ]++;
488       }
489     }
490     ierr = DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr);
491   }
492   if (ncells) { *ncells = nel; }
493   *count  = sum;
494   PetscFunctionReturn(0);
495 }
496