xref: /petsc/src/dm/impls/swarm/swarm_migrate.c (revision 480eef7b45975a89e98bc4cecdc59060a0a0bd97)
1 
2 #include <petsc/private/dmswarmimpl.h>    /*I   "petscdmswarm.h"   I*/
3 #include "data_bucket.h"
4 #include "data_ex.h"
5 
6 
7 /*
8  User loads desired location (MPI rank) into field DMSwarm_rank
9 */
10 #undef __FUNCT__
11 #define __FUNCT__ "DMSwarmMigrate_Push_Basic"
12 PetscErrorCode DMSwarmMigrate_Push_Basic(DM dm,PetscBool remove_sent_points)
13 {
14   DM_Swarm *swarm = (DM_Swarm*)dm->data;
15   PetscErrorCode ierr;
16   DataEx de;
17   PetscInt p,npoints,*rankval,n_points_recv;
18   PetscMPIInt rank,nrank;
19   void *point_buffer,*recv_points;
20   size_t sizeof_dmswarm_point;
21 
22   ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr);
23 
24   ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr);
25   ierr = DMSwarmGetField(dm,"DMSwarm_rank",NULL,NULL,(void**)&rankval);CHKERRQ(ierr);
26 
27   de = DataExCreate(PetscObjectComm((PetscObject)dm),0);CHKERRQ(ierr);
28 
29   ierr = DataExTopologyInitialize(de);CHKERRQ(ierr);
30   for (p=0; p<npoints; p++) {
31     nrank = rankval[p];
32     if (nrank != rank) {
33       ierr = DataExTopologyAddNeighbour(de,nrank);CHKERRQ(ierr);
34     }
35   }
36   ierr = DataExTopologyFinalize(de);CHKERRQ(ierr);
37 
38   ierr = DataExInitializeSendCount(de);CHKERRQ(ierr);
39   for (p=0; p<npoints; p++) {
40     nrank = rankval[p];
41     if (nrank != rank) {
42       ierr = DataExAddToSendCount(de,nrank,1);CHKERRQ(ierr);
43     }
44   }
45   ierr = DataExFinalizeSendCount(de);CHKERRQ(ierr);
46 
47   ierr = DataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);CHKERRQ(ierr);
48   ierr = DataExPackInitialize(de,sizeof_dmswarm_point);CHKERRQ(ierr);
49   for (p=0; p<npoints; p++) {
50     nrank = rankval[p];
51     if (nrank != rank) {
52       /* copy point into buffer */
53       ierr = DataBucketFillPackedArray(swarm->db,p,point_buffer);CHKERRQ(ierr);
54       /* insert point buffer into DataExchanger */
55       ierr = DataExPackData(de,nrank,1,point_buffer);CHKERRQ(ierr);
56     }
57   }
58   ierr = DataExPackFinalize(de);CHKERRQ(ierr);
59 
60 
61   if (remove_sent_points) {
62     /* remove points which left processor */
63     ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr);
64     for (p=0; p<npoints; p++) {
65       nrank = rankval[p];
66       if (nrank != rank) {
67         /* kill point */
68         ierr = DataBucketRemovePointAtIndex(swarm->db,p);CHKERRQ(ierr);
69         DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); /* you need to update npoints as the list size decreases! */
70         p--; /* check replacement point */
71       }
72     }
73   }
74   ierr = DMSwarmRestoreField(dm,"DMSwarm_rank",NULL,NULL,(void**)&rankval);CHKERRQ(ierr);
75 
76   ierr = DataExBegin(de);CHKERRQ(ierr);
77   ierr = DataExEnd(de);CHKERRQ(ierr);
78 
79   ierr = DataExGetRecvData(de,&n_points_recv,(void**)&recv_points);CHKERRQ(ierr);
80 
81 	ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr);
82 	ierr = DataBucketSetSizes(swarm->db,npoints + n_points_recv,-1);CHKERRQ(ierr);
83 	for (p=0; p<n_points_recv; p++) {
84     void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point );
85 
86     ierr = DataBucketInsertPackedArray(swarm->db,npoints+p,data_p);CHKERRQ(ierr);
87   }
88   ierr = DataExView(de);CHKERRQ(ierr);
89   ierr = DataBucketDestroyPackedArray(swarm->db,&point_buffer);CHKERRQ(ierr);
90   ierr = DataExDestroy(de);CHKERRQ(ierr);
91 
92   PetscFunctionReturn(0);
93 }
94 
95 //DMLocatePoints
96 
97 #undef __FUNCT__
98 #define __FUNCT__ "DMSwarmMigrate_CellDM"
99 PetscErrorCode DMSwarmMigrate_CellDM(DM dm,PetscBool remove_sent_points)
100 {
101   DM_Swarm *swarm = (DM_Swarm*)dm->data;
102   PetscErrorCode ierr;
103   PetscInt p,npoints,*rankval;
104   const PetscInt *LA_iscell;
105   DM dmcell;
106   IS iscell;
107   Vec pos;
108 
109   ierr = DMSwarmGetCellDM(dm,&dmcell);CHKERRQ(ierr);
110   if (!dmcell) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only valid if cell DM provided");
111 
112   ierr = DMSwarmCreateGlobalVectorFromField(dm,"coor",&pos);CHKERRQ(ierr);
113   ierr = DMLocatePoints(dmcell,pos,&iscell);CHKERRQ(ierr);
114   ierr = DMSwarmDestroyGlobalVectorFromField(dm,"coor",&pos);CHKERRQ(ierr);
115 
116   ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr);
117   ierr = DMSwarmGetField(dm,"DMSwarm_rank",NULL,NULL,(void**)&rankval);CHKERRQ(ierr);
118   ierr = ISGetIndices(iscell,&LA_iscell);CHKERRQ(ierr);
119   for (p=0; p<npoints; p++) {
120     if (LA_iscell[p] == -1) {
121       rankval[p] = -1;
122     }
123   }
124   ierr = ISRestoreIndices(iscell,&LA_iscell);CHKERRQ(ierr);
125   ierr = DMSwarmRestoreField(dm,"DMSwarm_rank",NULL,NULL,(void**)&rankval);CHKERRQ(ierr);
126   ierr = ISDestroy(&iscell);CHKERRQ(ierr);
127 
128   ierr = DMSwarmMigrate_Push_Basic(dm,remove_sent_points);CHKERRQ(ierr);
129 
130   PetscFunctionReturn(0);
131 }
132 
133 /*
134  Redundant as this assumes points can only be sent to a single rank
135 */
136 #undef __FUNCT__
137 #define __FUNCT__ "DMSwarmMigrate_GlobalToLocal_Basic"
138 PetscErrorCode DMSwarmMigrate_GlobalToLocal_Basic(DM dm,PetscInt *globalsize)
139 {
140   DM_Swarm *swarm = (DM_Swarm*)dm->data;
141   PetscErrorCode ierr;
142   DataEx de;
143   PetscInt p,npoints,*rankval,n_points_recv;
144   PetscMPIInt rank,nrank,negrank;
145   void *point_buffer,*recv_points;
146   size_t sizeof_dmswarm_point;
147 
148   ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr);
149 
150   ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr);
151   *globalsize = npoints;
152   ierr = DMSwarmGetField(dm,"DMSwarm_rank",NULL,NULL,(void**)&rankval);CHKERRQ(ierr);
153 
154   de = DataExCreate(PetscObjectComm((PetscObject)dm),0);CHKERRQ(ierr);
155 
156   ierr = DataExTopologyInitialize(de);CHKERRQ(ierr);
157   for (p=0; p<npoints; p++) {
158     negrank = rankval[p];
159     if (negrank < 0) {
160       nrank = -negrank - 1;
161       ierr = DataExTopologyAddNeighbour(de,nrank);CHKERRQ(ierr);
162     }
163   }
164   ierr = DataExTopologyFinalize(de);CHKERRQ(ierr);
165 
166   ierr = DataExInitializeSendCount(de);CHKERRQ(ierr);
167   for (p=0; p<npoints; p++) {
168     negrank = rankval[p];
169     if (negrank < 0) {
170       nrank = -negrank - 1;
171       ierr = DataExAddToSendCount(de,nrank,1);CHKERRQ(ierr);
172     }
173   }
174   ierr = DataExFinalizeSendCount(de);CHKERRQ(ierr);
175 
176   ierr = DataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);CHKERRQ(ierr);
177   ierr = DataExPackInitialize(de,sizeof_dmswarm_point);CHKERRQ(ierr);
178   for (p=0; p<npoints; p++) {
179     negrank = rankval[p];
180     if (negrank < 0) {
181       nrank = -negrank - 1;
182       rankval[p] = nrank;
183       /* copy point into buffer */
184       ierr = DataBucketFillPackedArray(swarm->db,p,point_buffer);CHKERRQ(ierr);
185       /* insert point buffer into DataExchanger */
186       ierr = DataExPackData(de,nrank,1,point_buffer);CHKERRQ(ierr);
187       rankval[p] = negrank;
188     }
189   }
190   ierr = DataExPackFinalize(de);CHKERRQ(ierr);
191 
192   ierr = DMSwarmRestoreField(dm,"DMSwarm_rank",NULL,NULL,(void**)&rankval);CHKERRQ(ierr);
193 
194   ierr = DataExBegin(de);CHKERRQ(ierr);
195   ierr = DataExEnd(de);CHKERRQ(ierr);
196 
197   ierr = DataExGetRecvData(de,&n_points_recv,(void**)&recv_points);CHKERRQ(ierr);
198 
199 	ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr);
200 	ierr = DataBucketSetSizes(swarm->db,npoints + n_points_recv,-1);CHKERRQ(ierr);
201 	for (p=0; p<n_points_recv; p++) {
202     void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point );
203 
204     ierr = DataBucketInsertPackedArray(swarm->db,npoints+p,data_p);CHKERRQ(ierr);
205   }
206   ierr = DataExView(de);CHKERRQ(ierr);
207   ierr = DataBucketDestroyPackedArray(swarm->db,&point_buffer);CHKERRQ(ierr);
208   ierr = DataExDestroy(de);CHKERRQ(ierr);
209 
210   PetscFunctionReturn(0);
211 }
212 
213 typedef struct {
214   PetscMPIInt owner_rank;
215   PetscReal min[3],max[3];
216 } CollectBBox;
217 
218 #undef __FUNCT__
219 #define __FUNCT__ "DMSwarmCollect_DMDABoundingBox"
220 PETSC_EXTERN PetscErrorCode DMSwarmCollect_DMDABoundingBox(DM dm,PetscInt *globalsize)
221 {
222   DM_Swarm *swarm = (DM_Swarm*)dm->data;
223   PetscErrorCode ierr;
224   DataEx de;
225   PetscInt p,pk,npoints,*rankval,n_points_recv,n_bbox_recv,dim,neighbour_cells;
226   PetscMPIInt rank,nrank;
227   void *point_buffer,*recv_points;
228   size_t sizeof_dmswarm_point,sizeof_bbox_ctx;
229   PetscBool isdmda;
230   CollectBBox *bbox,*recv_bbox;
231   const PetscMPIInt *dmneighborranks;
232   DM dmcell;
233 
234   ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr);
235 
236   ierr = DMSwarmGetCellDM(dm,&dmcell);CHKERRQ(ierr);
237   if (!dmcell) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only valid if cell DM provided");
238 
239   isdmda = PETSC_FALSE;
240   PetscObjectTypeCompare((PetscObject)dmcell,DMDA,&isdmda);
241   if (!isdmda) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only DMDA support for CollectBoundingBox");
242 
243   ierr = DMDAGetInfo(dm,&dim,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);CHKERRQ(ierr);
244   if (dim == 1) {
245     neighbour_cells = 3;
246   } else if (dim == 2) {
247     neighbour_cells = 9;
248   } else {
249     neighbour_cells = 27;
250   }
251 
252   sizeof_bbox_ctx = sizeof(CollectBBox);
253   PetscMalloc1(1,&bbox);
254   bbox->owner_rank = rank;
255 
256   /* compute the bounding box based on the overlapping / stenctil size */
257   //ierr = DMDAGetLocalBoundingBox(dmcell,bbox->min,bbox->max);CHKERRQ(ierr);
258   {
259     Vec lcoor;
260 
261     ierr = DMGetCoordinatesLocal(dmcell,&lcoor);CHKERRQ(ierr);
262 
263     if (dim >= 1) {
264       ierr = VecStrideMin(lcoor,0,NULL,&bbox->min[0]);CHKERRQ(ierr);
265       ierr = VecStrideMax(lcoor,0,NULL,&bbox->max[0]);CHKERRQ(ierr);
266     }
267 
268     if (dim >= 2) {
269       ierr = VecStrideMin(lcoor,1,NULL,&bbox->min[1]);CHKERRQ(ierr);
270       ierr = VecStrideMax(lcoor,1,NULL,&bbox->max[1]);CHKERRQ(ierr);
271     }
272 
273     if (dim == 3) {
274       ierr = VecStrideMin(lcoor,2,NULL,&bbox->min[2]);CHKERRQ(ierr);
275       ierr = VecStrideMax(lcoor,2,NULL,&bbox->max[2]);CHKERRQ(ierr);
276     }
277 }
278 
279   ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr);
280   *globalsize = npoints;
281   ierr = DMSwarmGetField(dm,"DMSwarm_rank",NULL,NULL,(void**)&rankval);CHKERRQ(ierr);
282 
283   de = DataExCreate(PetscObjectComm((PetscObject)dm),0);CHKERRQ(ierr);
284 
285   /* use DMDA neighbours */
286 	ierr = DMDAGetNeighbors(dmcell,&dmneighborranks);CHKERRQ(ierr);
287 
288   ierr = DataExTopologyInitialize(de);CHKERRQ(ierr);
289   for (p=0; p<neighbour_cells; p++) {
290 		if ( (dmneighborranks[p] >= 0) && (dmneighborranks[p] != rank) ) {
291       ierr = DataExTopologyAddNeighbour(de,dmneighborranks[p]);CHKERRQ(ierr);
292     }
293   }
294   ierr = DataExTopologyFinalize(de);CHKERRQ(ierr);
295 
296   ierr = DataExInitializeSendCount(de);CHKERRQ(ierr);
297   for (p=0; p<neighbour_cells; p++) {
298 		if ( (dmneighborranks[p] >= 0) && (dmneighborranks[p] != rank) ) {
299       ierr = DataExAddToSendCount(de,dmneighborranks[p],1);CHKERRQ(ierr);
300     }
301   }
302   ierr = DataExFinalizeSendCount(de);CHKERRQ(ierr);
303 
304 
305   /* send bounding boxes */
306   ierr = DataExPackInitialize(de,sizeof_bbox_ctx);CHKERRQ(ierr);
307   for (p=0; p<neighbour_cells; p++) {
308     nrank = dmneighborranks[p];
309 		if ( (nrank >= 0) && (nrank != rank) ) {
310       /* insert bbox buffer into DataExchanger */
311       ierr = DataExPackData(de,nrank,1,bbox);CHKERRQ(ierr);
312     }
313   }
314   ierr = DataExPackFinalize(de);CHKERRQ(ierr);
315 
316   /* recv bounding boxes */
317   ierr = DataExBegin(de);CHKERRQ(ierr);
318   ierr = DataExEnd(de);CHKERRQ(ierr);
319 
320   ierr = DataExGetRecvData(de,&n_bbox_recv,(void**)&recv_bbox);CHKERRQ(ierr);
321 
322 
323   for (p=0; p<n_bbox_recv; p++) {
324     printf("[rank %d]: box from %d : range[%+1.4e,%+1.4e]x[%+1.4e,%+1.4e]\n",rank,recv_bbox[p].owner_rank,
325            recv_bbox[p].min[0],recv_bbox[p].max[0],recv_bbox[p].min[1],recv_bbox[p].max[1]);
326   }
327 
328   /* of course this is stupid as this "generic" function should have a better way to know what the coordinates are called */
329   ierr = DataExInitializeSendCount(de);CHKERRQ(ierr);
330   for (pk=0; pk<n_bbox_recv; pk++) {
331     PetscReal *array_x,*array_y;
332 
333     ierr = DMSwarmGetField(dm,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr);
334     ierr = DMSwarmGetField(dm,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr);
335 
336     for (p=0; p<npoints; p++) {
337       if ((array_x[p] >= recv_bbox[pk].min[0]) && (array_x[p] <= recv_bbox[pk].max[0]) ) {
338         if ((array_y[p] >= recv_bbox[pk].min[1]) && (array_y[p] <= recv_bbox[pk].max[1]) ) {
339 
340           ierr = DataExAddToSendCount(de,recv_bbox[pk].owner_rank,1);CHKERRQ(ierr);
341 
342         }
343       }
344     }
345     ierr = DMSwarmRestoreField(dm,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr);
346     ierr = DMSwarmRestoreField(dm,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr);
347   }
348   ierr = DataExFinalizeSendCount(de);CHKERRQ(ierr);
349 
350 
351 
352   ierr = DataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);CHKERRQ(ierr);
353   ierr = DataExPackInitialize(de,sizeof_dmswarm_point);CHKERRQ(ierr);
354   for (pk=0; pk<n_bbox_recv; pk++) {
355     PetscReal *array_x,*array_y;
356 
357     ierr = DMSwarmGetField(dm,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr);
358     ierr = DMSwarmGetField(dm,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr);
359 
360     for (p=0; p<npoints; p++) {
361       if ((array_x[p] >= recv_bbox[pk].min[0]) && (array_x[p] <= recv_bbox[pk].max[0]) ) {
362         if ((array_y[p] >= recv_bbox[pk].min[1]) && (array_y[p] <= recv_bbox[pk].max[1]) ) {
363           // copy point into buffer //
364           ierr = DataBucketFillPackedArray(swarm->db,p,point_buffer);CHKERRQ(ierr);
365           // insert point buffer into DataExchanger //
366           ierr = DataExPackData(de,recv_bbox[pk].owner_rank,1,point_buffer);CHKERRQ(ierr);
367         }
368       }
369     }
370 
371     ierr = DMSwarmRestoreField(dm,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr);
372     ierr = DMSwarmRestoreField(dm,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr);
373   }
374 
375   ierr = DataExPackFinalize(de);CHKERRQ(ierr);
376 
377   ierr = DMSwarmRestoreField(dm,"DMSwarm_rank",NULL,NULL,(void**)&rankval);CHKERRQ(ierr);
378 
379   ierr = DataExBegin(de);CHKERRQ(ierr);
380   ierr = DataExEnd(de);CHKERRQ(ierr);
381 
382   ierr = DataExGetRecvData(de,&n_points_recv,(void**)&recv_points);CHKERRQ(ierr);
383 
384 	ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr);
385 	ierr = DataBucketSetSizes(swarm->db,npoints + n_points_recv,-1);CHKERRQ(ierr);
386 	for (p=0; p<n_points_recv; p++) {
387     void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point );
388 
389     ierr = DataBucketInsertPackedArray(swarm->db,npoints+p,data_p);CHKERRQ(ierr);
390   }
391 
392   ierr = DataBucketDestroyPackedArray(swarm->db,&point_buffer);CHKERRQ(ierr);
393   PetscFree(bbox);
394   ierr = DataExView(de);CHKERRQ(ierr);
395   ierr = DataExDestroy(de);CHKERRQ(ierr);
396 
397   PetscFunctionReturn(0);
398 }
399 
400 
401 /* General collection when no order, or neighbour information is provided */
402 /*
403  User provides context and collect() method
404  Broadcast user context
405 
406  for each context / rank {
407    collect(swarm,context,n,list)
408  }
409 
410 */
411 #undef __FUNCT__
412 #define __FUNCT__ "DMSwarmCollect_General"
413 PETSC_EXTERN PetscErrorCode DMSwarmCollect_General(DM dm,PetscErrorCode (*collect)(DM,void*,PetscInt*,PetscInt**),size_t ctx_size,void *ctx,PetscInt *globalsize)
414 {
415   DM_Swarm       *swarm = (DM_Swarm*)dm->data;
416   PetscErrorCode ierr;
417   DataEx         de;
418   PetscInt       p,r,npoints,n_points_recv;
419   PetscMPIInt    commsize,rank;
420   void           *point_buffer,*recv_points;
421   void           *ctxlist;
422   PetscInt       *n2collect,**collectlist;
423   size_t         sizeof_dmswarm_point;
424 
425   ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm),&commsize);CHKERRQ(ierr);
426   ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr);
427 
428   ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr);
429   *globalsize = npoints;
430 
431   /* Broadcast user context */
432   PetscMalloc(ctx_size*commsize,&ctxlist);
433   ierr = MPI_Allgather(ctx,ctx_size,MPI_CHAR,ctxlist,ctx_size,MPI_CHAR,PetscObjectComm((PetscObject)dm));CHKERRQ(ierr);
434 
435   PetscMalloc1(commsize,&n2collect);
436   PetscMalloc1(commsize,&collectlist);
437 
438   for (r=0; r<commsize; r++) {
439     PetscInt _n2collect;
440     PetscInt *_collectlist;
441     void     *_ctx_r;
442 
443     _n2collect   = 0;
444     _collectlist = NULL;
445 
446     if (r != rank) { /* don't collect data from yourself */
447       _ctx_r = (void*)( (char*)ctxlist + r * ctx_size );
448       ierr = collect(dm,_ctx_r,&_n2collect,&_collectlist);CHKERRQ(ierr);
449     }
450 
451     n2collect[r]   = _n2collect;
452     collectlist[r] = _collectlist;
453   }
454 
455   de = DataExCreate(PetscObjectComm((PetscObject)dm),0);CHKERRQ(ierr);
456 
457   /* Define topology */
458   ierr = DataExTopologyInitialize(de);CHKERRQ(ierr);
459   for (r=0; r<commsize; r++) {
460     if (n2collect[r] > 0) {
461       ierr = DataExTopologyAddNeighbour(de,(PetscMPIInt)r);CHKERRQ(ierr);
462     }
463   }
464   ierr = DataExTopologyFinalize(de);CHKERRQ(ierr);
465 
466   /* Define send counts */
467   ierr = DataExInitializeSendCount(de);CHKERRQ(ierr);
468   for (r=0; r<commsize; r++) {
469     if (n2collect[r] > 0) {
470       ierr = DataExAddToSendCount(de,r,n2collect[r]);CHKERRQ(ierr);
471     }
472   }
473   ierr = DataExFinalizeSendCount(de);CHKERRQ(ierr);
474 
475   /* Pack data */
476   ierr = DataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);CHKERRQ(ierr);
477 
478   ierr = DataExPackInitialize(de,sizeof_dmswarm_point);CHKERRQ(ierr);
479   for (r=0; r<commsize; r++) {
480 
481     for (p=0; p<n2collect[r]; p++) {
482       ierr = DataBucketFillPackedArray(swarm->db,collectlist[r][p],point_buffer);CHKERRQ(ierr);
483       /* insert point buffer into the data exchanger */
484       ierr = DataExPackData(de,r,1,point_buffer);CHKERRQ(ierr);
485     }
486   }
487 
488   ierr = DataExPackFinalize(de);CHKERRQ(ierr);
489 
490   /* Scatter */
491   ierr = DataExBegin(de);CHKERRQ(ierr);
492   ierr = DataExEnd(de);CHKERRQ(ierr);
493 
494   /* Collect data in DMSwarm container */
495   ierr = DataExGetRecvData(de,&n_points_recv,(void**)&recv_points);CHKERRQ(ierr);
496 
497 	ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr);
498 	ierr = DataBucketSetSizes(swarm->db,npoints + n_points_recv,-1);CHKERRQ(ierr);
499 	for (p=0; p<n_points_recv; p++) {
500     void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point );
501 
502     ierr = DataBucketInsertPackedArray(swarm->db,npoints+p,data_p);CHKERRQ(ierr);
503   }
504 
505   /* Release memory */
506   for (r=0; r<commsize; r++) {
507     if (collectlist[r]) PetscFree(collectlist[r]);
508   }
509   PetscFree(collectlist);
510   PetscFree(n2collect);
511   PetscFree(ctxlist);
512   ierr = DataBucketDestroyPackedArray(swarm->db,&point_buffer);CHKERRQ(ierr);
513   ierr = DataExView(de);CHKERRQ(ierr);
514   ierr = DataExDestroy(de);CHKERRQ(ierr);
515 
516   PetscFunctionReturn(0);
517 }
518 
519