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