xref: /petsc/src/dm/impls/swarm/swarm_migrate.c (revision 009b43ef5b40b325ad97cf697c109976f0bd0406) 
1 
2 #include <petscsf.h>
3 #include <petscdmswarm.h>
4 #include <petsc/private/dmswarmimpl.h>    /*I   "petscdmswarm.h"   I*/
5 #include "data_bucket.h"
6 #include "data_ex.h"
7 
8 
9 /*
10  User loads desired location (MPI rank) into field DMSwarm_rank
11 */
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   PetscFunctionBegin;
23   ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr);
24 
25   ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr);
26   ierr = DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr);
27   ierr = DataExCreate(PetscObjectComm((PetscObject)dm),0, &de);CHKERRQ(ierr);
28   ierr = DataExTopologyInitialize(de);CHKERRQ(ierr);
29   for (p = 0; p < npoints; ++p) {
30     nrank = rankval[p];
31     if (nrank != rank) {
32       ierr = DataExTopologyAddNeighbour(de,nrank);CHKERRQ(ierr);
33     }
34   }
35   ierr = DataExTopologyFinalize(de);CHKERRQ(ierr);
36   ierr = DataExInitializeSendCount(de);CHKERRQ(ierr);
37   for (p=0; p<npoints; p++) {
38     nrank = rankval[p];
39     if (nrank != rank) {
40       ierr = DataExAddToSendCount(de,nrank,1);CHKERRQ(ierr);
41     }
42   }
43   ierr = DataExFinalizeSendCount(de);CHKERRQ(ierr);
44   ierr = DataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);CHKERRQ(ierr);
45   ierr = DataExPackInitialize(de,sizeof_dmswarm_point);CHKERRQ(ierr);
46   for (p=0; p<npoints; p++) {
47     nrank = rankval[p];
48     if (nrank != rank) {
49       /* copy point into buffer */
50       ierr = DataBucketFillPackedArray(swarm->db,p,point_buffer);CHKERRQ(ierr);
51       /* insert point buffer into DataExchanger */
52       ierr = DataExPackData(de,nrank,1,point_buffer);CHKERRQ(ierr);
53     }
54   }
55   ierr = DataExPackFinalize(de);CHKERRQ(ierr);
56   ierr = DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr);
57 
58   if (remove_sent_points) {
59     DataField gfield;
60 
61     ierr = DataBucketGetDataFieldByName(swarm->db,DMSwarmField_rank,&gfield);CHKERRQ(ierr);
62     ierr = DataFieldGetAccess(gfield);CHKERRQ(ierr);
63     ierr = DataFieldGetEntries(gfield,(void**)&rankval);CHKERRQ(ierr);
64 
65     /* remove points which left processor */
66     ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr);
67     for (p=0; p<npoints; p++) {
68       nrank = rankval[p];
69       if (nrank != rank) {
70         /* kill point */
71         ierr = DataFieldRestoreAccess(gfield);CHKERRQ(ierr);
72 
73         ierr = DataBucketRemovePointAtIndex(swarm->db,p);CHKERRQ(ierr);
74 
75         ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); /* you need to update npoints as the list size decreases! */
76         ierr = DataFieldGetAccess(gfield);CHKERRQ(ierr);
77         ierr = DataFieldGetEntries(gfield,(void**)&rankval);CHKERRQ(ierr);
78         p--; /* check replacement point */
79       }
80     }
81     ierr = DataFieldRestoreEntries(gfield,(void**)&rankval);CHKERRQ(ierr);
82     ierr = DataFieldRestoreAccess(gfield);CHKERRQ(ierr);
83   }
84   ierr = DataExBegin(de);CHKERRQ(ierr);
85   ierr = DataExEnd(de);CHKERRQ(ierr);
86   ierr = DataExGetRecvData(de,&n_points_recv,(void**)&recv_points);CHKERRQ(ierr);
87   ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr);
88   ierr = DataBucketSetSizes(swarm->db,npoints + n_points_recv,DATA_BUCKET_BUFFER_DEFAULT);CHKERRQ(ierr);
89   for (p=0; p<n_points_recv; p++) {
90     void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point );
91 
92     ierr = DataBucketInsertPackedArray(swarm->db,npoints+p,data_p);CHKERRQ(ierr);
93   }
94   ierr = DataExView(de);CHKERRQ(ierr);
95   ierr = DataBucketDestroyPackedArray(swarm->db,&point_buffer);CHKERRQ(ierr);
96   ierr = DataExDestroy(de);CHKERRQ(ierr);
97   PetscFunctionReturn(0);
98 }
99 
100 PetscErrorCode DMSwarmMigrate_DMNeighborScatter(DM dm,DM dmcell,PetscBool remove_sent_points,PetscInt *npoints_prior_migration)
101 {
102   DM_Swarm *swarm = (DM_Swarm*)dm->data;
103   PetscErrorCode ierr;
104   DataEx de;
105   PetscInt r,p,npoints,*rankval,n_points_recv;
106   PetscMPIInt rank,_rank;
107   const PetscMPIInt *neighbourranks;
108   void *point_buffer,*recv_points;
109   size_t sizeof_dmswarm_point;
110   PetscInt nneighbors;
111   PetscMPIInt mynneigh,*myneigh;
112 
113   PetscFunctionBegin;
114   ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr);
115   ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr);
116   ierr = DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr);
117   ierr = DataExCreate(PetscObjectComm((PetscObject)dm),0,&de);CHKERRQ(ierr);
118   ierr = DMGetNeighbors(dmcell,&nneighbors,&neighbourranks);CHKERRQ(ierr);
119   ierr = DataExTopologyInitialize(de);CHKERRQ(ierr);
120   for (r=0; r<nneighbors; r++) {
121     _rank = neighbourranks[r];
122     if ((_rank != rank) && (_rank > 0)) {
123       ierr = DataExTopologyAddNeighbour(de,_rank);CHKERRQ(ierr);
124     }
125   }
126   ierr = DataExTopologyFinalize(de);CHKERRQ(ierr);
127   ierr = DataExTopologyGetNeighbours(de,&mynneigh,&myneigh);CHKERRQ(ierr);
128   ierr = DataExInitializeSendCount(de);CHKERRQ(ierr);
129   for (p=0; p<npoints; p++) {
130     if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
131       for (r=0; r<mynneigh; r++) {
132         _rank = myneigh[r];
133         ierr = DataExAddToSendCount(de,_rank,1);CHKERRQ(ierr);
134       }
135     }
136   }
137   ierr = DataExFinalizeSendCount(de);CHKERRQ(ierr);
138   ierr = DataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);CHKERRQ(ierr);
139   ierr = DataExPackInitialize(de,sizeof_dmswarm_point);CHKERRQ(ierr);
140   for (p=0; p<npoints; p++) {
141     if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
142       for (r=0; r<mynneigh; r++) {
143         _rank = myneigh[r];
144         /* copy point into buffer */
145         ierr = DataBucketFillPackedArray(swarm->db,p,point_buffer);CHKERRQ(ierr);
146         /* insert point buffer into DataExchanger */
147         ierr = DataExPackData(de,_rank,1,point_buffer);CHKERRQ(ierr);
148       }
149     }
150   }
151   ierr = DataExPackFinalize(de);CHKERRQ(ierr);
152   ierr = DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr);
153   if (remove_sent_points) {
154     DataField PField;
155 
156     ierr = DataBucketGetDataFieldByName(swarm->db,DMSwarmField_rank,&PField);CHKERRQ(ierr);
157     ierr = DataFieldGetEntries(PField,(void**)&rankval);CHKERRQ(ierr);
158     /* remove points which left processor */
159     ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr);
160     for (p=0; p<npoints; p++) {
161       if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
162         /* kill point */
163         ierr = DataBucketRemovePointAtIndex(swarm->db,p);CHKERRQ(ierr);
164         ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); /* you need to update npoints as the list size decreases! */
165         ierr = DataFieldGetEntries(PField,(void**)&rankval);CHKERRQ(ierr); /* update date point increase realloc performed */
166         p--; /* check replacement point */
167       }
168     }
169   }
170   ierr = DataBucketGetSizes(swarm->db,npoints_prior_migration,NULL,NULL);CHKERRQ(ierr);
171   ierr = DataExBegin(de);CHKERRQ(ierr);
172   ierr = DataExEnd(de);CHKERRQ(ierr);
173   ierr = DataExGetRecvData(de,&n_points_recv,(void**)&recv_points);CHKERRQ(ierr);
174   ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr);
175   ierr = DataBucketSetSizes(swarm->db,npoints + n_points_recv,DATA_BUCKET_BUFFER_DEFAULT);CHKERRQ(ierr);
176   for (p=0; p<n_points_recv; p++) {
177     void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point );
178 
179     ierr = DataBucketInsertPackedArray(swarm->db,npoints+p,data_p);CHKERRQ(ierr);
180   }
181   ierr = DataBucketDestroyPackedArray(swarm->db,&point_buffer);CHKERRQ(ierr);
182   ierr = DataExDestroy(de);CHKERRQ(ierr);
183   PetscFunctionReturn(0);
184 }
185 
186 PetscErrorCode DMSwarmMigrate_CellDMScatter(DM dm,PetscBool remove_sent_points)
187 {
188   DM_Swarm *swarm = (DM_Swarm*)dm->data;
189   PetscErrorCode ierr;
190   PetscInt p,npoints,npointsg=0,npoints2,npoints2g,*rankval,npoints_prior_migration;
191   PetscSF sfcell = NULL;
192   const PetscSFNode *LA_sfcell;
193   DM dmcell;
194   Vec pos;
195   PetscBool error_check = swarm->migrate_error_on_missing_point;
196   PetscMPIInt commsize,rank;
197 
198   PetscFunctionBegin;
199   ierr = DMSwarmGetCellDM(dm,&dmcell);CHKERRQ(ierr);
200   if (!dmcell) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only valid if cell DM provided");
201 
202   ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm),&commsize);CHKERRQ(ierr);
203   ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr);
204 
205   ierr = DMSwarmCreateLocalVectorFromField(dm, DMSwarmPICField_coor, &pos);CHKERRQ(ierr);
206   ierr = DMLocatePoints(dmcell, pos, DM_POINTLOCATION_NONE, &sfcell);CHKERRQ(ierr);
207   ierr = DMSwarmDestroyLocalVectorFromField(dm, DMSwarmPICField_coor, &pos);CHKERRQ(ierr);
208 
209   if (error_check) {
210     ierr = DMSwarmGetSize(dm,&npointsg);CHKERRQ(ierr);
211   }
212   ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr);
213   ierr = DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr);
214   ierr = PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell);CHKERRQ(ierr);
215   for (p=0; p<npoints; p++) {
216     rankval[p] = LA_sfcell[p].index;
217   }
218   ierr = DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr);
219   ierr = PetscSFDestroy(&sfcell);CHKERRQ(ierr);
220 
221   if (commsize > 1) {
222     ierr = DMSwarmMigrate_DMNeighborScatter(dm,dmcell,remove_sent_points,&npoints_prior_migration);CHKERRQ(ierr);
223   } else {
224     ierr = DMSwarmGetSize(dm,&npoints_prior_migration);CHKERRQ(ierr);
225   }
226 
227   /* locate points newly recevied */
228   ierr = DataBucketGetSizes(swarm->db,&npoints2,NULL,NULL);CHKERRQ(ierr);
229 
230 #if 0
231   { // safe alternative - however this performs two point locations on: (i) the intial points set and; (ii) the (intial + recieved) point set
232     PetscScalar *LA_coor;
233     PetscInt bs;
234     DataField PField;
235 
236     ierr = DMSwarmGetField(dm,DMSwarmPICField_coor,&bs,NULL,(void**)&LA_coor);CHKERRQ(ierr);
237     ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,bs,bs*npoints2,(const PetscScalar*)LA_coor,&pos);CHKERRQ(ierr);
238     ierr = DMLocatePoints(dmcell,pos,DM_POINTLOCATION_NONE,&sfcell);CHKERRQ(ierr);
239 
240     ierr = VecDestroy(&pos);CHKERRQ(ierr);
241     ierr = DMSwarmRestoreField(dm,DMSwarmPICField_coor,&bs,NULL,(void**)&LA_coor);CHKERRQ(ierr);
242 
243     ierr = PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell);CHKERRQ(ierr);
244     ierr = DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr);
245     for (p=0; p<npoints2; p++) {
246       rankval[p] = LA_sfcell[p].index;
247     }
248     ierr = PetscSFDestroy(&sfcell);CHKERRQ(ierr);
249     ierr = DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr);
250 
251     /* remove points which left processor */
252     ierr = DataBucketGetDataFieldByName(swarm->db,DMSwarmField_rank,&PField);CHKERRQ(ierr);
253     ierr = DataFieldGetEntries(PField,(void**)&rankval);CHKERRQ(ierr);
254 
255     ierr = DataBucketGetSizes(swarm->db,&npoints2,NULL,NULL);CHKERRQ(ierr);
256     for (p=0; p<npoints2; p++) {
257       if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
258         /* kill point */
259         ierr = DataBucketRemovePointAtIndex(swarm->db,p);CHKERRQ(ierr);
260         ierr = DataBucketGetSizes(swarm->db,&npoints2,NULL,NULL);CHKERRQ(ierr); /* you need to update npoints as the list size decreases! */
261         ierr = DataFieldGetEntries(PField,(void**)&rankval);CHKERRQ(ierr); /* update date point increase realloc performed */
262         p--; /* check replacement point */
263       }
264     }
265   }
266 #endif
267 
268   { // this performs two point locations: (i) on the intial points set prior to communication; and (ii) on the new (recieved) points
269     PetscScalar *LA_coor;
270     PetscInt npoints_from_neighbours,bs;
271     DataField PField;
272 
273     npoints_from_neighbours = npoints2 - npoints_prior_migration;
274 
275     ierr = DMSwarmGetField(dm,DMSwarmPICField_coor,&bs,NULL,(void**)&LA_coor);CHKERRQ(ierr);
276     ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,bs,bs*npoints_from_neighbours,(const PetscScalar*)&LA_coor[bs*npoints_prior_migration],&pos);CHKERRQ(ierr);
277 
278     ierr = DMLocatePoints(dmcell,pos,DM_POINTLOCATION_NONE,&sfcell);CHKERRQ(ierr);
279 
280     ierr = VecDestroy(&pos);CHKERRQ(ierr);
281     ierr = DMSwarmRestoreField(dm,DMSwarmPICField_coor,&bs,NULL,(void**)&LA_coor);CHKERRQ(ierr);
282 
283     ierr = PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell);CHKERRQ(ierr);
284     ierr = DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr);
285     for (p=0; p<npoints_from_neighbours; p++) {
286       rankval[npoints_prior_migration + p] = LA_sfcell[p].index;
287     }
288     ierr = DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr);
289     ierr = PetscSFDestroy(&sfcell);CHKERRQ(ierr);
290 
291     /* remove points which left processor */
292     ierr = DataBucketGetDataFieldByName(swarm->db,DMSwarmField_rank,&PField);CHKERRQ(ierr);
293     ierr = DataFieldGetEntries(PField,(void**)&rankval);CHKERRQ(ierr);
294 
295     ierr = DataBucketGetSizes(swarm->db,&npoints2,NULL,NULL);CHKERRQ(ierr);
296     for (p=npoints_prior_migration; p<npoints2; p++) {
297       if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
298         /* kill point */
299         ierr = DataBucketRemovePointAtIndex(swarm->db,p);CHKERRQ(ierr);
300         ierr = DataBucketGetSizes(swarm->db,&npoints2,NULL,NULL);CHKERRQ(ierr); /* you need to update npoints as the list size decreases! */
301         ierr = DataFieldGetEntries(PField,(void**)&rankval);CHKERRQ(ierr); /* update date point increase realloc performed */
302         p--; /* check replacement point */
303       }
304     }
305   }
306 
307   /* check for error on removed points */
308   if (error_check) {
309     ierr = DMSwarmGetSize(dm,&npoints2g);CHKERRQ(ierr);
310     if (npointsg != npoints2g) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Points from the DMSwarm must remain constant during migration (initial %D - final %D)",npointsg,npoints2g);
311   }
312   PetscFunctionReturn(0);
313 }
314 
315 PetscErrorCode DMSwarmMigrate_CellDMExact(DM dm,PetscBool remove_sent_points)
316 {
317   PetscFunctionBegin;
318   PetscFunctionReturn(0);
319 }
320 
321 /*
322  Redundant as this assumes points can only be sent to a single rank
323 */
324 PetscErrorCode DMSwarmMigrate_GlobalToLocal_Basic(DM dm,PetscInt *globalsize)
325 {
326   DM_Swarm *swarm = (DM_Swarm*)dm->data;
327   PetscErrorCode ierr;
328   DataEx de;
329   PetscInt p,npoints,*rankval,n_points_recv;
330   PetscMPIInt rank,nrank,negrank;
331   void *point_buffer,*recv_points;
332   size_t sizeof_dmswarm_point;
333 
334   PetscFunctionBegin;
335   ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr);
336   ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr);
337   *globalsize = npoints;
338   ierr = DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr);
339   ierr = DataExCreate(PetscObjectComm((PetscObject)dm),0,&de);CHKERRQ(ierr);
340   ierr = DataExTopologyInitialize(de);CHKERRQ(ierr);
341   for (p=0; p<npoints; p++) {
342     negrank = rankval[p];
343     if (negrank < 0) {
344       nrank = -negrank - 1;
345       ierr = DataExTopologyAddNeighbour(de,nrank);CHKERRQ(ierr);
346     }
347   }
348   ierr = DataExTopologyFinalize(de);CHKERRQ(ierr);
349   ierr = DataExInitializeSendCount(de);CHKERRQ(ierr);
350   for (p=0; p<npoints; p++) {
351     negrank = rankval[p];
352     if (negrank < 0) {
353       nrank = -negrank - 1;
354       ierr = DataExAddToSendCount(de,nrank,1);CHKERRQ(ierr);
355     }
356   }
357   ierr = DataExFinalizeSendCount(de);CHKERRQ(ierr);
358   ierr = DataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);CHKERRQ(ierr);
359   ierr = DataExPackInitialize(de,sizeof_dmswarm_point);CHKERRQ(ierr);
360   for (p=0; p<npoints; p++) {
361     negrank = rankval[p];
362     if (negrank < 0) {
363       nrank = -negrank - 1;
364       rankval[p] = nrank;
365       /* copy point into buffer */
366       ierr = DataBucketFillPackedArray(swarm->db,p,point_buffer);CHKERRQ(ierr);
367       /* insert point buffer into DataExchanger */
368       ierr = DataExPackData(de,nrank,1,point_buffer);CHKERRQ(ierr);
369       rankval[p] = negrank;
370     }
371   }
372   ierr = DataExPackFinalize(de);CHKERRQ(ierr);
373   ierr = DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr);
374   ierr = DataExBegin(de);CHKERRQ(ierr);
375   ierr = DataExEnd(de);CHKERRQ(ierr);
376   ierr = DataExGetRecvData(de,&n_points_recv,(void**)&recv_points);CHKERRQ(ierr);
377   ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr);
378   ierr = DataBucketSetSizes(swarm->db,npoints + n_points_recv,DATA_BUCKET_BUFFER_DEFAULT);CHKERRQ(ierr);
379   for (p=0; p<n_points_recv; p++) {
380     void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point );
381 
382     ierr = DataBucketInsertPackedArray(swarm->db,npoints+p,data_p);CHKERRQ(ierr);
383   }
384   ierr = DataExView(de);CHKERRQ(ierr);
385   ierr = DataBucketDestroyPackedArray(swarm->db,&point_buffer);CHKERRQ(ierr);
386   ierr = DataExDestroy(de);CHKERRQ(ierr);
387   PetscFunctionReturn(0);
388 }
389 
390 typedef struct {
391   PetscMPIInt owner_rank;
392   PetscReal min[3],max[3];
393 } CollectBBox;
394 
395 PETSC_EXTERN PetscErrorCode DMSwarmCollect_DMDABoundingBox(DM dm,PetscInt *globalsize)
396 {
397   DM_Swarm *swarm = (DM_Swarm*)dm->data;
398   PetscErrorCode ierr;
399   DataEx de;
400   PetscInt p,pk,npoints,*rankval,n_points_recv,n_bbox_recv,dim,neighbour_cells;
401   PetscMPIInt rank,nrank;
402   void *point_buffer,*recv_points;
403   size_t sizeof_dmswarm_point,sizeof_bbox_ctx;
404   PetscBool isdmda;
405   CollectBBox *bbox,*recv_bbox;
406   const PetscMPIInt *dmneighborranks;
407   DM dmcell;
408 
409   PetscFunctionBegin;
410   ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr);
411 
412   ierr = DMSwarmGetCellDM(dm,&dmcell);CHKERRQ(ierr);
413   if (!dmcell) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only valid if cell DM provided");
414   isdmda = PETSC_FALSE;
415   PetscObjectTypeCompare((PetscObject)dmcell,DMDA,&isdmda);
416   if (!isdmda) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only DMDA support for CollectBoundingBox");
417 
418   ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr);
419   sizeof_bbox_ctx = sizeof(CollectBBox);
420   PetscMalloc1(1,&bbox);
421   bbox->owner_rank = rank;
422 
423   /* compute the bounding box based on the overlapping / stenctil size */
424   {
425     Vec lcoor;
426 
427     ierr = DMGetCoordinatesLocal(dmcell,&lcoor);CHKERRQ(ierr);
428     if (dim >= 1) {
429       ierr = VecStrideMin(lcoor,0,NULL,&bbox->min[0]);CHKERRQ(ierr);
430       ierr = VecStrideMax(lcoor,0,NULL,&bbox->max[0]);CHKERRQ(ierr);
431     }
432     if (dim >= 2) {
433       ierr = VecStrideMin(lcoor,1,NULL,&bbox->min[1]);CHKERRQ(ierr);
434       ierr = VecStrideMax(lcoor,1,NULL,&bbox->max[1]);CHKERRQ(ierr);
435     }
436     if (dim == 3) {
437       ierr = VecStrideMin(lcoor,2,NULL,&bbox->min[2]);CHKERRQ(ierr);
438       ierr = VecStrideMax(lcoor,2,NULL,&bbox->max[2]);CHKERRQ(ierr);
439     }
440   }
441   ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr);
442   *globalsize = npoints;
443   ierr = DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr);
444   ierr = DataExCreate(PetscObjectComm((PetscObject)dm),0,&de);CHKERRQ(ierr);
445   /* use DMDA neighbours */
446   ierr = DMDAGetNeighbors(dmcell,&dmneighborranks);CHKERRQ(ierr);
447   if (dim == 1) {
448     neighbour_cells = 3;
449   } else if (dim == 2) {
450     neighbour_cells = 9;
451   } else {
452     neighbour_cells = 27;
453   }
454   ierr = DataExTopologyInitialize(de);CHKERRQ(ierr);
455   for (p=0; p<neighbour_cells; p++) {
456     if ( (dmneighborranks[p] >= 0) && (dmneighborranks[p] != rank) ) {
457       ierr = DataExTopologyAddNeighbour(de,dmneighborranks[p]);CHKERRQ(ierr);
458     }
459   }
460   ierr = DataExTopologyFinalize(de);CHKERRQ(ierr);
461   ierr = DataExInitializeSendCount(de);CHKERRQ(ierr);
462   for (p=0; p<neighbour_cells; p++) {
463     if ( (dmneighborranks[p] >= 0) && (dmneighborranks[p] != rank) ) {
464       ierr = DataExAddToSendCount(de,dmneighborranks[p],1);CHKERRQ(ierr);
465     }
466   }
467   ierr = DataExFinalizeSendCount(de);CHKERRQ(ierr);
468   /* send bounding boxes */
469   ierr = DataExPackInitialize(de,sizeof_bbox_ctx);CHKERRQ(ierr);
470   for (p=0; p<neighbour_cells; p++) {
471     nrank = dmneighborranks[p];
472     if ( (nrank >= 0) && (nrank != rank) ) {
473       /* insert bbox buffer into DataExchanger */
474       ierr = DataExPackData(de,nrank,1,bbox);CHKERRQ(ierr);
475     }
476   }
477   ierr = DataExPackFinalize(de);CHKERRQ(ierr);
478   /* recv bounding boxes */
479   ierr = DataExBegin(de);CHKERRQ(ierr);
480   ierr = DataExEnd(de);CHKERRQ(ierr);
481   ierr = DataExGetRecvData(de,&n_bbox_recv,(void**)&recv_bbox);CHKERRQ(ierr);
482   for (p=0; p<n_bbox_recv; p++) {
483     PetscPrintf(PETSC_COMM_SELF,"[rank %d]: box from %d : range[%+1.4e,%+1.4e]x[%+1.4e,%+1.4e]\n",rank,recv_bbox[p].owner_rank,
484            (double)recv_bbox[p].min[0],(double)recv_bbox[p].max[0],(double)recv_bbox[p].min[1],(double)recv_bbox[p].max[1]);
485   }
486   /* of course this is stupid as this "generic" function should have a better way to know what the coordinates are called */
487   ierr = DataExInitializeSendCount(de);CHKERRQ(ierr);
488   for (pk=0; pk<n_bbox_recv; pk++) {
489     PetscReal *array_x,*array_y;
490 
491     ierr = DMSwarmGetField(dm,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr);
492     ierr = DMSwarmGetField(dm,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr);
493     for (p=0; p<npoints; p++) {
494       if ((array_x[p] >= recv_bbox[pk].min[0]) && (array_x[p] <= recv_bbox[pk].max[0]) ) {
495         if ((array_y[p] >= recv_bbox[pk].min[1]) && (array_y[p] <= recv_bbox[pk].max[1]) ) {
496           ierr = DataExAddToSendCount(de,recv_bbox[pk].owner_rank,1);CHKERRQ(ierr);
497         }
498       }
499     }
500     ierr = DMSwarmRestoreField(dm,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr);
501     ierr = DMSwarmRestoreField(dm,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr);
502   }
503   ierr = DataExFinalizeSendCount(de);CHKERRQ(ierr);
504   ierr = DataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);CHKERRQ(ierr);
505   ierr = DataExPackInitialize(de,sizeof_dmswarm_point);CHKERRQ(ierr);
506   for (pk=0; pk<n_bbox_recv; pk++) {
507     PetscReal *array_x,*array_y;
508 
509     ierr = DMSwarmGetField(dm,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr);
510     ierr = DMSwarmGetField(dm,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr);
511     for (p=0; p<npoints; p++) {
512       if ((array_x[p] >= recv_bbox[pk].min[0]) && (array_x[p] <= recv_bbox[pk].max[0]) ) {
513         if ((array_y[p] >= recv_bbox[pk].min[1]) && (array_y[p] <= recv_bbox[pk].max[1]) ) {
514           /* copy point into buffer */
515           ierr = DataBucketFillPackedArray(swarm->db,p,point_buffer);CHKERRQ(ierr);
516           /* insert point buffer into DataExchanger */
517           ierr = DataExPackData(de,recv_bbox[pk].owner_rank,1,point_buffer);CHKERRQ(ierr);
518         }
519       }
520     }
521     ierr = DMSwarmRestoreField(dm,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr);
522     ierr = DMSwarmRestoreField(dm,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr);
523   }
524   ierr = DataExPackFinalize(de);CHKERRQ(ierr);
525   ierr = DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr);
526   ierr = DataExBegin(de);CHKERRQ(ierr);
527   ierr = DataExEnd(de);CHKERRQ(ierr);
528   ierr = DataExGetRecvData(de,&n_points_recv,(void**)&recv_points);CHKERRQ(ierr);
529   ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr);
530   ierr = DataBucketSetSizes(swarm->db,npoints + n_points_recv,DATA_BUCKET_BUFFER_DEFAULT);CHKERRQ(ierr);
531   for (p=0; p<n_points_recv; p++) {
532     void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point );
533 
534     ierr = DataBucketInsertPackedArray(swarm->db,npoints+p,data_p);CHKERRQ(ierr);
535   }
536   ierr = DataBucketDestroyPackedArray(swarm->db,&point_buffer);CHKERRQ(ierr);
537   PetscFree(bbox);
538   ierr = DataExView(de);CHKERRQ(ierr);
539   ierr = DataExDestroy(de);CHKERRQ(ierr);
540   PetscFunctionReturn(0);
541 }
542 
543 
544 /* General collection when no order, or neighbour information is provided */
545 /*
546  User provides context and collect() method
547  Broadcast user context
548 
549  for each context / rank {
550    collect(swarm,context,n,list)
551  }
552 */
553 PETSC_EXTERN PetscErrorCode DMSwarmCollect_General(DM dm,PetscErrorCode (*collect)(DM,void*,PetscInt*,PetscInt**),size_t ctx_size,void *ctx,PetscInt *globalsize)
554 {
555   DM_Swarm       *swarm = (DM_Swarm*)dm->data;
556   PetscErrorCode ierr;
557   DataEx         de;
558   PetscInt       p,r,npoints,n_points_recv;
559   PetscMPIInt    commsize,rank;
560   void           *point_buffer,*recv_points;
561   void           *ctxlist;
562   PetscInt       *n2collect,**collectlist;
563   size_t         sizeof_dmswarm_point;
564 
565   PetscFunctionBegin;
566   ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm),&commsize);CHKERRQ(ierr);
567   ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr);
568   ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr);
569   *globalsize = npoints;
570   /* Broadcast user context */
571   PetscMalloc(ctx_size*commsize,&ctxlist);
572   ierr = MPI_Allgather(ctx,ctx_size,MPI_CHAR,ctxlist,ctx_size,MPI_CHAR,PetscObjectComm((PetscObject)dm));CHKERRQ(ierr);
573   ierr = PetscMalloc1(commsize,&n2collect);CHKERRQ(ierr);
574   ierr = PetscMalloc1(commsize,&collectlist);CHKERRQ(ierr);
575   for (r=0; r<commsize; r++) {
576     PetscInt _n2collect;
577     PetscInt *_collectlist;
578     void     *_ctx_r;
579 
580     _n2collect   = 0;
581     _collectlist = NULL;
582     if (r != rank) { /* don't collect data from yourself */
583       _ctx_r = (void*)( (char*)ctxlist + r * ctx_size );
584       ierr = collect(dm,_ctx_r,&_n2collect,&_collectlist);CHKERRQ(ierr);
585     }
586     n2collect[r]   = _n2collect;
587     collectlist[r] = _collectlist;
588   }
589   ierr = DataExCreate(PetscObjectComm((PetscObject)dm),0,&de);CHKERRQ(ierr);
590   /* Define topology */
591   ierr = DataExTopologyInitialize(de);CHKERRQ(ierr);
592   for (r=0; r<commsize; r++) {
593     if (n2collect[r] > 0) {
594       ierr = DataExTopologyAddNeighbour(de,(PetscMPIInt)r);CHKERRQ(ierr);
595     }
596   }
597   ierr = DataExTopologyFinalize(de);CHKERRQ(ierr);
598   /* Define send counts */
599   ierr = DataExInitializeSendCount(de);CHKERRQ(ierr);
600   for (r=0; r<commsize; r++) {
601     if (n2collect[r] > 0) {
602       ierr = DataExAddToSendCount(de,r,n2collect[r]);CHKERRQ(ierr);
603     }
604   }
605   ierr = DataExFinalizeSendCount(de);CHKERRQ(ierr);
606   /* Pack data */
607   ierr = DataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);CHKERRQ(ierr);
608   ierr = DataExPackInitialize(de,sizeof_dmswarm_point);CHKERRQ(ierr);
609   for (r=0; r<commsize; r++) {
610     for (p=0; p<n2collect[r]; p++) {
611       ierr = DataBucketFillPackedArray(swarm->db,collectlist[r][p],point_buffer);CHKERRQ(ierr);
612       /* insert point buffer into the data exchanger */
613       ierr = DataExPackData(de,r,1,point_buffer);CHKERRQ(ierr);
614     }
615   }
616   ierr = DataExPackFinalize(de);CHKERRQ(ierr);
617   /* Scatter */
618   ierr = DataExBegin(de);CHKERRQ(ierr);
619   ierr = DataExEnd(de);CHKERRQ(ierr);
620   /* Collect data in DMSwarm container */
621   ierr = DataExGetRecvData(de,&n_points_recv,(void**)&recv_points);CHKERRQ(ierr);
622   ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr);
623   ierr = DataBucketSetSizes(swarm->db,npoints + n_points_recv,DATA_BUCKET_BUFFER_DEFAULT);CHKERRQ(ierr);
624   for (p=0; p<n_points_recv; p++) {
625     void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point );
626 
627     ierr = DataBucketInsertPackedArray(swarm->db,npoints+p,data_p);CHKERRQ(ierr);
628   }
629   /* Release memory */
630   for (r=0; r<commsize; r++) {
631     if (collectlist[r]) PetscFree(collectlist[r]);
632   }
633   ierr = PetscFree(collectlist);CHKERRQ(ierr);
634   ierr = PetscFree(n2collect);CHKERRQ(ierr);
635   ierr = PetscFree(ctxlist);CHKERRQ(ierr);
636   ierr = DataBucketDestroyPackedArray(swarm->db,&point_buffer);CHKERRQ(ierr);
637   ierr = DataExView(de);CHKERRQ(ierr);
638   ierr = DataExDestroy(de);CHKERRQ(ierr);
639   PetscFunctionReturn(0);
640 }
641 
642