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