Actual source code: swarm_migrate.c

petsc-3.13.6 2020-09-29
Report Typos and Errors
  1:  #include <petscsf.h>
  2:  #include <petscdmswarm.h>
  3:  #include <petscdmda.h>
  4:  #include <petsc/private/dmswarmimpl.h>
  5: #include "../src/dm/impls/swarm/data_bucket.h"
  6: #include "../src/dm/impls/swarm/data_ex.h"

  8: /*
  9:  User loads desired location (MPI rank) into field DMSwarm_rank
 10: */
 11: PetscErrorCode DMSwarmMigrate_Push_Basic(DM dm,PetscBool remove_sent_points)
 12: {
 13:   DM_Swarm       *swarm = (DM_Swarm*)dm->data;
 15:   DMSwarmDataEx  de;
 16:   PetscInt       p,npoints,*rankval,n_points_recv;
 17:   PetscMPIInt    rank,nrank;
 18:   void           *point_buffer,*recv_points;
 19:   size_t         sizeof_dmswarm_point;

 22:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);

 24:   DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
 25:   DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
 26:   DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm),0, &de);
 27:   DMSwarmDataExTopologyInitialize(de);
 28:   for (p = 0; p < npoints; ++p) {
 29:     nrank = rankval[p];
 30:     if (nrank != rank) {
 31:       DMSwarmDataExTopologyAddNeighbour(de,nrank);
 32:     }
 33:   }
 34:   DMSwarmDataExTopologyFinalize(de);
 35:   DMSwarmDataExInitializeSendCount(de);
 36:   for (p=0; p<npoints; p++) {
 37:     nrank = rankval[p];
 38:     if (nrank != rank) {
 39:       DMSwarmDataExAddToSendCount(de,nrank,1);
 40:     }
 41:   }
 42:   DMSwarmDataExFinalizeSendCount(de);
 43:   DMSwarmDataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);
 44:   DMSwarmDataExPackInitialize(de,sizeof_dmswarm_point);
 45:   for (p=0; p<npoints; p++) {
 46:     nrank = rankval[p];
 47:     if (nrank != rank) {
 48:       /* copy point into buffer */
 49:       DMSwarmDataBucketFillPackedArray(swarm->db,p,point_buffer);
 50:       /* insert point buffer into DMSwarmDataExchanger */
 51:       DMSwarmDataExPackData(de,nrank,1,point_buffer);
 52:     }
 53:   }
 54:   DMSwarmDataExPackFinalize(de);
 55:   DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);

 57:   if (remove_sent_points) {
 58:     DMSwarmDataField gfield;

 60:     DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,DMSwarmField_rank,&gfield);
 61:     DMSwarmDataFieldGetAccess(gfield);
 62:     DMSwarmDataFieldGetEntries(gfield,(void**)&rankval);

 64:     /* remove points which left processor */
 65:     DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
 66:     for (p=0; p<npoints; p++) {
 67:       nrank = rankval[p];
 68:       if (nrank != rank) {
 69:         /* kill point */
 70:         DMSwarmDataFieldRestoreAccess(gfield);
 71:         
 72:         DMSwarmDataBucketRemovePointAtIndex(swarm->db,p);

 74:         DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL); /* you need to update npoints as the list size decreases! */
 75:         DMSwarmDataFieldGetAccess(gfield);
 76:         DMSwarmDataFieldGetEntries(gfield,(void**)&rankval);
 77:         p--; /* check replacement point */
 78:       }
 79:     }
 80:     DMSwarmDataFieldRestoreEntries(gfield,(void**)&rankval);
 81:     DMSwarmDataFieldRestoreAccess(gfield);
 82:   }
 83:   DMSwarmDataExBegin(de);
 84:   DMSwarmDataExEnd(de);
 85:   DMSwarmDataExGetRecvData(de,&n_points_recv,(void**)&recv_points);
 86:   DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
 87:   DMSwarmDataBucketSetSizes(swarm->db,npoints + n_points_recv,DMSWARM_DATA_BUCKET_BUFFER_DEFAULT);
 88:   for (p=0; p<n_points_recv; p++) {
 89:     void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point );

 91:     DMSwarmDataBucketInsertPackedArray(swarm->db,npoints+p,data_p);
 92:   }
 93:   DMSwarmDataExView(de);
 94:   DMSwarmDataBucketDestroyPackedArray(swarm->db,&point_buffer);
 95:   DMSwarmDataExDestroy(de);
 96:   return(0);
 97: }

 99: PetscErrorCode DMSwarmMigrate_DMNeighborScatter(DM dm,DM dmcell,PetscBool remove_sent_points,PetscInt *npoints_prior_migration)
100: {
101:   DM_Swarm          *swarm = (DM_Swarm*)dm->data;
102:   PetscErrorCode    ierr;
103:   DMSwarmDataEx     de;
104:   PetscInt          r,p,npoints,*rankval,n_points_recv;
105:   PetscMPIInt       rank,_rank;
106:   const PetscMPIInt *neighbourranks;
107:   void              *point_buffer,*recv_points;
108:   size_t            sizeof_dmswarm_point;
109:   PetscInt          nneighbors;
110:   PetscMPIInt       mynneigh,*myneigh;

113:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
114:   DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
115:   DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
116:   DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm),0,&de);
117:   DMGetNeighbors(dmcell,&nneighbors,&neighbourranks);
118:   DMSwarmDataExTopologyInitialize(de);
119:   for (r=0; r<nneighbors; r++) {
120:     _rank = neighbourranks[r];
121:     if ((_rank != rank) && (_rank > 0)) {
122:       DMSwarmDataExTopologyAddNeighbour(de,_rank);
123:     }
124:   }
125:   DMSwarmDataExTopologyFinalize(de);
126:   DMSwarmDataExTopologyGetNeighbours(de,&mynneigh,&myneigh);
127:   DMSwarmDataExInitializeSendCount(de);
128:   for (p=0; p<npoints; p++) {
129:     if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
130:       for (r=0; r<mynneigh; r++) {
131:         _rank = myneigh[r];
132:         DMSwarmDataExAddToSendCount(de,_rank,1);
133:       }
134:     }
135:   }
136:   DMSwarmDataExFinalizeSendCount(de);
137:   DMSwarmDataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);
138:   DMSwarmDataExPackInitialize(de,sizeof_dmswarm_point);
139:   for (p=0; p<npoints; p++) {
140:     if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
141:       for (r=0; r<mynneigh; r++) {
142:         _rank = myneigh[r];
143:         /* copy point into buffer */
144:         DMSwarmDataBucketFillPackedArray(swarm->db,p,point_buffer);
145:         /* insert point buffer into DMSwarmDataExchanger */
146:         DMSwarmDataExPackData(de,_rank,1,point_buffer);
147:       }
148:     }
149:   }
150:   DMSwarmDataExPackFinalize(de);
151:   DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
152:   if (remove_sent_points) {
153:     DMSwarmDataField PField;

155:     DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,DMSwarmField_rank,&PField);
156:     DMSwarmDataFieldGetEntries(PField,(void**)&rankval);
157:     /* remove points which left processor */
158:     DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
159:     for (p=0; p<npoints; p++) {
160:       if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
161:         /* kill point */
162:         DMSwarmDataBucketRemovePointAtIndex(swarm->db,p);
163:         DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL); /* you need to update npoints as the list size decreases! */
164:         DMSwarmDataFieldGetEntries(PField,(void**)&rankval); /* update date point increase realloc performed */
165:         p--; /* check replacement point */
166:       }
167:     }
168:   }
169:   DMSwarmDataBucketGetSizes(swarm->db,npoints_prior_migration,NULL,NULL);
170:   DMSwarmDataExBegin(de);
171:   DMSwarmDataExEnd(de);
172:   DMSwarmDataExGetRecvData(de,&n_points_recv,(void**)&recv_points);
173:   DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
174:   DMSwarmDataBucketSetSizes(swarm->db,npoints + n_points_recv,DMSWARM_DATA_BUCKET_BUFFER_DEFAULT);
175:   for (p=0; p<n_points_recv; p++) {
176:     void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point );

178:     DMSwarmDataBucketInsertPackedArray(swarm->db,npoints+p,data_p);
179:   }
180:   DMSwarmDataBucketDestroyPackedArray(swarm->db,&point_buffer);
181:   DMSwarmDataExDestroy(de);
182:   return(0);
183: }

185: PetscErrorCode DMSwarmMigrate_CellDMScatter(DM dm,PetscBool remove_sent_points)
186: {
187:   DM_Swarm          *swarm = (DM_Swarm*)dm->data;
188:   PetscErrorCode    ierr;
189:   PetscInt          p,npoints,npointsg=0,npoints2,npoints2g,*rankval,npoints_prior_migration;
190:   PetscSF           sfcell = NULL;
191:   const PetscSFNode *LA_sfcell;
192:   DM                dmcell;
193:   Vec               pos;
194:   PetscBool         error_check = swarm->migrate_error_on_missing_point;
195:   PetscMPIInt       size,rank;

198:   DMSwarmGetCellDM(dm,&dmcell);
199:   if (!dmcell) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only valid if cell DM provided");

201:   MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);
202:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);

204: #if 1
205:   {
206:     PetscInt *p_cellid;
207:     PetscInt npoints_curr,range = 0;
208:     PetscSFNode *sf_cells;

210:     
211:     DMSwarmDataBucketGetSizes(swarm->db,&npoints_curr,NULL,NULL);
212:     PetscMalloc1(npoints_curr, &sf_cells);

214:     DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
215:     DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&p_cellid);
216:     for (p=0; p<npoints_curr; p++) {

218:       sf_cells[p].rank  = 0;
219:       sf_cells[p].index = p_cellid[p];
220:       if (p_cellid[p] > range) {
221:         range = p_cellid[p];
222:       }
223:     }
224:     DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&p_cellid);
225:     DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);

227:     /*PetscSFCreate(PetscObjectComm((PetscObject)dm),&sfcell);*/
228:     PetscSFCreate(PETSC_COMM_SELF,&sfcell);
229:     PetscSFSetGraph(sfcell, range, npoints_curr, NULL, PETSC_OWN_POINTER, sf_cells, PETSC_OWN_POINTER);
230:   }
231: #endif
232:   
233:   DMSwarmCreateLocalVectorFromField(dm, DMSwarmPICField_coor, &pos);
234:   DMLocatePoints(dmcell, pos, DM_POINTLOCATION_NONE, &sfcell);
235:   DMSwarmDestroyLocalVectorFromField(dm, DMSwarmPICField_coor, &pos);

237:   if (error_check) {
238:     DMSwarmGetSize(dm,&npointsg);
239:   }
240:   DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
241:   DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
242:   PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell);
243:   for (p=0; p<npoints; p++) {
244:     rankval[p] = LA_sfcell[p].index;
245:   }
246:   DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
247:   PetscSFDestroy(&sfcell);

249:   if (size > 1) {
250:     DMSwarmMigrate_DMNeighborScatter(dm,dmcell,remove_sent_points,&npoints_prior_migration);
251:   } else {
252:     DMSwarmDataField PField;
253:     PetscInt npoints_curr;
254:     
255:     /* remove points which the domain */
256:     DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,DMSwarmField_rank,&PField);
257:     DMSwarmDataFieldGetEntries(PField,(void**)&rankval);
258:     
259:     DMSwarmDataBucketGetSizes(swarm->db,&npoints_curr,NULL,NULL);
260:     for (p=0; p<npoints_curr; p++) {
261:       if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
262:         /* kill point */
263:         DMSwarmDataBucketRemovePointAtIndex(swarm->db,p);
264:         DMSwarmDataBucketGetSizes(swarm->db,&npoints_curr,NULL,NULL); /* you need to update npoints as the list size decreases! */
265:         DMSwarmDataFieldGetEntries(PField,(void**)&rankval); /* update date point increase realloc performed */
266:         p--; /* check replacement point */
267:       }
268:     }
269:     DMSwarmGetSize(dm,&npoints_prior_migration);
270:     
271:   }

273:   /* locate points newly received */
274:   DMSwarmDataBucketGetSizes(swarm->db,&npoints2,NULL,NULL);
275:   
276: #if 0
277:   { /* safe alternative - however this performs two point locations on: (i) the initial points set and; (ii) the (initial + received) point set */
278:     PetscScalar *LA_coor;
279:     PetscInt bs;
280:     DMSwarmDataField PField;

282:     DMSwarmGetField(dm,DMSwarmPICField_coor,&bs,NULL,(void**)&LA_coor);
283:     VecCreateSeqWithArray(PETSC_COMM_SELF,bs,bs*npoints2,(const PetscScalar*)LA_coor,&pos);
284:     DMLocatePoints(dmcell,pos,DM_POINTLOCATION_NONE,&sfcell);

286:     VecDestroy(&pos);
287:     DMSwarmRestoreField(dm,DMSwarmPICField_coor,&bs,NULL,(void**)&LA_coor);

289:     PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell);
290:     DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
291:     for (p=0; p<npoints2; p++) {
292:       rankval[p] = LA_sfcell[p].index;
293:     }
294:     PetscSFDestroy(&sfcell);
295:     DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);

297:     /* remove points which left processor */
298:     DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,DMSwarmField_rank,&PField);
299:     DMSwarmDataFieldGetEntries(PField,(void**)&rankval);

301:     DMSwarmDataBucketGetSizes(swarm->db,&npoints2,NULL,NULL);
302:     for (p=0; p<npoints2; p++) {
303:       if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
304:         /* kill point */
305:         DMSwarmDataBucketRemovePointAtIndex(swarm->db,p);
306:         DMSwarmDataBucketGetSizes(swarm->db,&npoints2,NULL,NULL); /* you need to update npoints as the list size decreases! */
307:         DMSwarmDataFieldGetEntries(PField,(void**)&rankval); /* update date point increase realloc performed */
308:         p--; /* check replacement point */
309:       }
310:     }
311:   }
312: #endif

314:   { /* this performs two point locations: (i) on the initial points set prior to communication; and (ii) on the new (received) points */
315:     PetscScalar      *LA_coor;
316:     PetscInt         npoints_from_neighbours,bs;
317:     DMSwarmDataField PField;
318:     
319:     npoints_from_neighbours = npoints2 - npoints_prior_migration;
320:     
321:     DMSwarmGetField(dm,DMSwarmPICField_coor,&bs,NULL,(void**)&LA_coor);
322:     VecCreateSeqWithArray(PETSC_COMM_SELF,bs,bs*npoints_from_neighbours,(const PetscScalar*)&LA_coor[bs*npoints_prior_migration],&pos);

324:     DMLocatePoints(dmcell,pos,DM_POINTLOCATION_NONE,&sfcell);
325:     
326:     VecDestroy(&pos);
327:     DMSwarmRestoreField(dm,DMSwarmPICField_coor,&bs,NULL,(void**)&LA_coor);
328:     
329:     PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell);
330:     DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
331:     for (p=0; p<npoints_from_neighbours; p++) {
332:       rankval[npoints_prior_migration + p] = LA_sfcell[p].index;
333:     }
334:     DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
335:     PetscSFDestroy(&sfcell);
336:     
337:     /* remove points which left processor */
338:     DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,DMSwarmField_rank,&PField);
339:     DMSwarmDataFieldGetEntries(PField,(void**)&rankval);
340:     
341:     DMSwarmDataBucketGetSizes(swarm->db,&npoints2,NULL,NULL);
342:     for (p=npoints_prior_migration; p<npoints2; p++) {
343:       if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
344:         /* kill point */
345:         DMSwarmDataBucketRemovePointAtIndex(swarm->db,p);
346:         DMSwarmDataBucketGetSizes(swarm->db,&npoints2,NULL,NULL); /* you need to update npoints as the list size decreases! */
347:         DMSwarmDataFieldGetEntries(PField,(void**)&rankval); /* update date point increase realloc performed */
348:         p--; /* check replacement point */
349:       }
350:     }
351:   }
352:   
353:   {
354:     PetscInt *p_cellid;
355:     
356:     DMSwarmDataBucketGetSizes(swarm->db,&npoints2,NULL,NULL);
357:     DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
358:     DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&p_cellid);
359:     for (p=0; p<npoints2; p++) {
360:       p_cellid[p] = rankval[p];
361:     }
362:     DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&p_cellid);
363:     DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
364:   }
365:   
366:   /* check for error on removed points */
367:   if (error_check) {
368:     DMSwarmGetSize(dm,&npoints2g);
369:     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);
370:   }
371:   return(0);
372: }

374: PetscErrorCode DMSwarmMigrate_CellDMExact(DM dm,PetscBool remove_sent_points)
375: {
377:   return(0);
378: }

380: /*
381:  Redundant as this assumes points can only be sent to a single rank
382: */
383: PetscErrorCode DMSwarmMigrate_GlobalToLocal_Basic(DM dm,PetscInt *globalsize)
384: {
385:   DM_Swarm       *swarm = (DM_Swarm*)dm->data;
387:   DMSwarmDataEx  de;
388:   PetscInt       p,npoints,*rankval,n_points_recv;
389:   PetscMPIInt    rank,nrank,negrank;
390:   void           *point_buffer,*recv_points;
391:   size_t         sizeof_dmswarm_point;

394:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
395:   DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
396:   *globalsize = npoints;
397:   DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
398:   DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm),0,&de);
399:   DMSwarmDataExTopologyInitialize(de);
400:   for (p=0; p<npoints; p++) {
401:     negrank = rankval[p];
402:     if (negrank < 0) {
403:       nrank = -negrank - 1;
404:       DMSwarmDataExTopologyAddNeighbour(de,nrank);
405:     }
406:   }
407:   DMSwarmDataExTopologyFinalize(de);
408:   DMSwarmDataExInitializeSendCount(de);
409:   for (p=0; p<npoints; p++) {
410:     negrank = rankval[p];
411:     if (negrank < 0) {
412:       nrank = -negrank - 1;
413:       DMSwarmDataExAddToSendCount(de,nrank,1);
414:     }
415:   }
416:   DMSwarmDataExFinalizeSendCount(de);
417:   DMSwarmDataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);
418:   DMSwarmDataExPackInitialize(de,sizeof_dmswarm_point);
419:   for (p=0; p<npoints; p++) {
420:     negrank = rankval[p];
421:     if (negrank < 0) {
422:       nrank = -negrank - 1;
423:       rankval[p] = nrank;
424:       /* copy point into buffer */
425:       DMSwarmDataBucketFillPackedArray(swarm->db,p,point_buffer);
426:       /* insert point buffer into DMSwarmDataExchanger */
427:       DMSwarmDataExPackData(de,nrank,1,point_buffer);
428:       rankval[p] = negrank;
429:     }
430:   }
431:   DMSwarmDataExPackFinalize(de);
432:   DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
433:   DMSwarmDataExBegin(de);
434:   DMSwarmDataExEnd(de);
435:   DMSwarmDataExGetRecvData(de,&n_points_recv,(void**)&recv_points);
436:   DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
437:   DMSwarmDataBucketSetSizes(swarm->db,npoints + n_points_recv,DMSWARM_DATA_BUCKET_BUFFER_DEFAULT);
438:   for (p=0; p<n_points_recv; p++) {
439:     void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point );

441:     DMSwarmDataBucketInsertPackedArray(swarm->db,npoints+p,data_p);
442:   }
443:   DMSwarmDataExView(de);
444:   DMSwarmDataBucketDestroyPackedArray(swarm->db,&point_buffer);
445:   DMSwarmDataExDestroy(de);
446:   return(0);
447: }

449: typedef struct {
450:   PetscMPIInt owner_rank;
451:   PetscReal   min[3],max[3];
452: } CollectBBox;

454: PETSC_EXTERN PetscErrorCode DMSwarmCollect_DMDABoundingBox(DM dm,PetscInt *globalsize)
455: {
456:   DM_Swarm *        swarm = (DM_Swarm*)dm->data;
457:   PetscErrorCode    ierr;
458:   DMSwarmDataEx     de;
459:   PetscInt          p,pk,npoints,*rankval,n_points_recv,n_bbox_recv,dim,neighbour_cells;
460:   PetscMPIInt       rank,nrank;
461:   void              *point_buffer,*recv_points;
462:   size_t            sizeof_dmswarm_point,sizeof_bbox_ctx;
463:   PetscBool         isdmda;
464:   CollectBBox       *bbox,*recv_bbox;
465:   const PetscMPIInt *dmneighborranks;
466:   DM                dmcell;

469:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);

471:   DMSwarmGetCellDM(dm,&dmcell);
472:   if (!dmcell) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only valid if cell DM provided");
473:   isdmda = PETSC_FALSE;
474:   PetscObjectTypeCompare((PetscObject)dmcell,DMDA,&isdmda);
475:   if (!isdmda) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only DMDA support for CollectBoundingBox");

477:   DMGetDimension(dm,&dim);
478:   sizeof_bbox_ctx = sizeof(CollectBBox);
479:   PetscMalloc1(1,&bbox);
480:   bbox->owner_rank = rank;

482:   /* compute the bounding box based on the overlapping / stenctil size */
483:   {
484:     Vec lcoor;

486:     DMGetCoordinatesLocal(dmcell,&lcoor);
487:     if (dim >= 1) {
488:       VecStrideMin(lcoor,0,NULL,&bbox->min[0]);
489:       VecStrideMax(lcoor,0,NULL,&bbox->max[0]);
490:     }
491:     if (dim >= 2) {
492:       VecStrideMin(lcoor,1,NULL,&bbox->min[1]);
493:       VecStrideMax(lcoor,1,NULL,&bbox->max[1]);
494:     }
495:     if (dim == 3) {
496:       VecStrideMin(lcoor,2,NULL,&bbox->min[2]);
497:       VecStrideMax(lcoor,2,NULL,&bbox->max[2]);
498:     }
499:   }
500:   DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
501:   *globalsize = npoints;
502:   DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
503:   DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm),0,&de);
504:   /* use DMDA neighbours */
505:   DMDAGetNeighbors(dmcell,&dmneighborranks);
506:   if (dim == 1) {
507:     neighbour_cells = 3;
508:   } else if (dim == 2) {
509:     neighbour_cells = 9;
510:   } else {
511:     neighbour_cells = 27;
512:   }
513:   DMSwarmDataExTopologyInitialize(de);
514:   for (p=0; p<neighbour_cells; p++) {
515:     if ( (dmneighborranks[p] >= 0) && (dmneighborranks[p] != rank) ) {
516:       DMSwarmDataExTopologyAddNeighbour(de,dmneighborranks[p]);
517:     }
518:   }
519:   DMSwarmDataExTopologyFinalize(de);
520:   DMSwarmDataExInitializeSendCount(de);
521:   for (p=0; p<neighbour_cells; p++) {
522:     if ( (dmneighborranks[p] >= 0) && (dmneighborranks[p] != rank) ) {
523:       DMSwarmDataExAddToSendCount(de,dmneighborranks[p],1);
524:     }
525:   }
526:   DMSwarmDataExFinalizeSendCount(de);
527:   /* send bounding boxes */
528:   DMSwarmDataExPackInitialize(de,sizeof_bbox_ctx);
529:   for (p=0; p<neighbour_cells; p++) {
530:     nrank = dmneighborranks[p];
531:     if ( (nrank >= 0) && (nrank != rank) ) {
532:       /* insert bbox buffer into DMSwarmDataExchanger */
533:       DMSwarmDataExPackData(de,nrank,1,bbox);
534:     }
535:   }
536:   DMSwarmDataExPackFinalize(de);
537:   /* recv bounding boxes */
538:   DMSwarmDataExBegin(de);
539:   DMSwarmDataExEnd(de);
540:   DMSwarmDataExGetRecvData(de,&n_bbox_recv,(void**)&recv_bbox);
541:   /*  Wrong, should not be using PETSC_COMM_WORLD */
542:   for (p=0; p<n_bbox_recv; p++) {
543:     PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[rank %d]: box from %d : range[%+1.4e,%+1.4e]x[%+1.4e,%+1.4e]\n",rank,recv_bbox[p].owner_rank,
544:            (double)recv_bbox[p].min[0],(double)recv_bbox[p].max[0],(double)recv_bbox[p].min[1],(double)recv_bbox[p].max[1]);
545:   }
546:   PetscSynchronizedFlush(PETSC_COMM_WORLD,stdout);
547:   /* of course this is stupid as this "generic" function should have a better way to know what the coordinates are called */
548:   DMSwarmDataExInitializeSendCount(de);
549:   for (pk=0; pk<n_bbox_recv; pk++) {
550:     PetscReal *array_x,*array_y;

552:     DMSwarmGetField(dm,"coorx",NULL,NULL,(void**)&array_x);
553:     DMSwarmGetField(dm,"coory",NULL,NULL,(void**)&array_y);
554:     for (p=0; p<npoints; p++) {
555:       if ((array_x[p] >= recv_bbox[pk].min[0]) && (array_x[p] <= recv_bbox[pk].max[0]) ) {
556:         if ((array_y[p] >= recv_bbox[pk].min[1]) && (array_y[p] <= recv_bbox[pk].max[1]) ) {
557:           DMSwarmDataExAddToSendCount(de,recv_bbox[pk].owner_rank,1);
558:         }
559:       }
560:     }
561:     DMSwarmRestoreField(dm,"coory",NULL,NULL,(void**)&array_y);
562:     DMSwarmRestoreField(dm,"coorx",NULL,NULL,(void**)&array_x);
563:   }
564:   DMSwarmDataExFinalizeSendCount(de);
565:   DMSwarmDataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);
566:   DMSwarmDataExPackInitialize(de,sizeof_dmswarm_point);
567:   for (pk=0; pk<n_bbox_recv; pk++) {
568:     PetscReal *array_x,*array_y;

570:     DMSwarmGetField(dm,"coorx",NULL,NULL,(void**)&array_x);
571:     DMSwarmGetField(dm,"coory",NULL,NULL,(void**)&array_y);
572:     for (p=0; p<npoints; p++) {
573:       if ((array_x[p] >= recv_bbox[pk].min[0]) && (array_x[p] <= recv_bbox[pk].max[0]) ) {
574:         if ((array_y[p] >= recv_bbox[pk].min[1]) && (array_y[p] <= recv_bbox[pk].max[1]) ) {
575:           /* copy point into buffer */
576:           DMSwarmDataBucketFillPackedArray(swarm->db,p,point_buffer);
577:           /* insert point buffer into DMSwarmDataExchanger */
578:           DMSwarmDataExPackData(de,recv_bbox[pk].owner_rank,1,point_buffer);
579:         }
580:       }
581:     }
582:     DMSwarmRestoreField(dm,"coory",NULL,NULL,(void**)&array_y);
583:     DMSwarmRestoreField(dm,"coorx",NULL,NULL,(void**)&array_x);
584:   }
585:   DMSwarmDataExPackFinalize(de);
586:   DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
587:   DMSwarmDataExBegin(de);
588:   DMSwarmDataExEnd(de);
589:   DMSwarmDataExGetRecvData(de,&n_points_recv,(void**)&recv_points);
590:   DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
591:   DMSwarmDataBucketSetSizes(swarm->db,npoints + n_points_recv,DMSWARM_DATA_BUCKET_BUFFER_DEFAULT);
592:   for (p=0; p<n_points_recv; p++) {
593:     void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point );

595:     DMSwarmDataBucketInsertPackedArray(swarm->db,npoints+p,data_p);
596:   }
597:   DMSwarmDataBucketDestroyPackedArray(swarm->db,&point_buffer);
598:   PetscFree(bbox);
599:   DMSwarmDataExView(de);
600:   DMSwarmDataExDestroy(de);
601:   return(0);
602: }


605: /* General collection when no order, or neighbour information is provided */
606: /*
607:  User provides context and collect() method
608:  Broadcast user context

610:  for each context / rank {
611:    collect(swarm,context,n,list)
612:  }
613: */
614: PETSC_EXTERN PetscErrorCode DMSwarmCollect_General(DM dm,PetscErrorCode (*collect)(DM,void*,PetscInt*,PetscInt**),size_t ctx_size,void *ctx,PetscInt *globalsize)
615: {
616:   DM_Swarm       *swarm = (DM_Swarm*)dm->data;
618:   DMSwarmDataEx         de;
619:   PetscInt       p,r,npoints,n_points_recv;
620:   PetscMPIInt    size,rank;
621:   void           *point_buffer,*recv_points;
622:   void           *ctxlist;
623:   PetscInt       *n2collect,**collectlist;
624:   size_t         sizeof_dmswarm_point;

627:   MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);
628:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
629:   DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
630:   *globalsize = npoints;
631:   /* Broadcast user context */
632:   PetscMalloc(ctx_size*size,&ctxlist);
633:   MPI_Allgather(ctx,ctx_size,MPI_CHAR,ctxlist,ctx_size,MPI_CHAR,PetscObjectComm((PetscObject)dm));
634:   PetscMalloc1(size,&n2collect);
635:   PetscMalloc1(size,&collectlist);
636:   for (r=0; r<size; r++) {
637:     PetscInt _n2collect;
638:     PetscInt *_collectlist;
639:     void     *_ctx_r;

641:     _n2collect   = 0;
642:     _collectlist = NULL;
643:     if (r != rank) { /* don't collect data from yourself */
644:       _ctx_r = (void*)( (char*)ctxlist + r * ctx_size );
645:       collect(dm,_ctx_r,&_n2collect,&_collectlist);
646:     }
647:     n2collect[r]   = _n2collect;
648:     collectlist[r] = _collectlist;
649:   }
650:   DMSwarmDataExCreate(PetscObjectComm((PetscObject)dm),0,&de);
651:   /* Define topology */
652:   DMSwarmDataExTopologyInitialize(de);
653:   for (r=0; r<size; r++) {
654:     if (n2collect[r] > 0) {
655:       DMSwarmDataExTopologyAddNeighbour(de,(PetscMPIInt)r);
656:     }
657:   }
658:   DMSwarmDataExTopologyFinalize(de);
659:   /* Define send counts */
660:   DMSwarmDataExInitializeSendCount(de);
661:   for (r=0; r<size; r++) {
662:     if (n2collect[r] > 0) {
663:       DMSwarmDataExAddToSendCount(de,r,n2collect[r]);
664:     }
665:   }
666:   DMSwarmDataExFinalizeSendCount(de);
667:   /* Pack data */
668:   DMSwarmDataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);
669:   DMSwarmDataExPackInitialize(de,sizeof_dmswarm_point);
670:   for (r=0; r<size; r++) {
671:     for (p=0; p<n2collect[r]; p++) {
672:       DMSwarmDataBucketFillPackedArray(swarm->db,collectlist[r][p],point_buffer);
673:       /* insert point buffer into the data exchanger */
674:       DMSwarmDataExPackData(de,r,1,point_buffer);
675:     }
676:   }
677:   DMSwarmDataExPackFinalize(de);
678:   /* Scatter */
679:   DMSwarmDataExBegin(de);
680:   DMSwarmDataExEnd(de);
681:   /* Collect data in DMSwarm container */
682:   DMSwarmDataExGetRecvData(de,&n_points_recv,(void**)&recv_points);
683:   DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
684:   DMSwarmDataBucketSetSizes(swarm->db,npoints + n_points_recv,DMSWARM_DATA_BUCKET_BUFFER_DEFAULT);
685:   for (p=0; p<n_points_recv; p++) {
686:     void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point );

688:     DMSwarmDataBucketInsertPackedArray(swarm->db,npoints+p,data_p);
689:   }
690:   /* Release memory */
691:   for (r=0; r<size; r++) {
692:     if (collectlist[r]) PetscFree(collectlist[r]);
693:   }
694:   PetscFree(collectlist);
695:   PetscFree(n2collect);
696:   PetscFree(ctxlist);
697:   DMSwarmDataBucketDestroyPackedArray(swarm->db,&point_buffer);
698:   DMSwarmDataExView(de);
699:   DMSwarmDataExDestroy(de);
700:   return(0);
701: }