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