Actual source code: swarm_migrate.c

  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;
 14:   DMSwarmDataEx  de;
 15:   PetscInt       p,npoints,*rankval,n_points_recv;
 16:   PetscMPIInt    rank,nrank;
 17:   void           *point_buffer,*recv_points;
 18:   size_t         sizeof_dmswarm_point;
 19:   PetscBool      debug = PETSC_FALSE;

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

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

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

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

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

 71:         DMSwarmDataBucketRemovePointAtIndex(swarm->db,p);

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

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

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

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

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

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

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

193:   DMSwarmGetCellDM(dm,&dmcell);

196:   MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);
197:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);

199: #if 1
200:   {
201:     PetscInt *p_cellid;
202:     PetscInt npoints_curr,range = 0;
203:     PetscSFNode *sf_cells;

205:     DMSwarmDataBucketGetSizes(swarm->db,&npoints_curr,NULL,NULL);
206:     PetscMalloc1(npoints_curr, &sf_cells);

208:     DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
209:     DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&p_cellid);
210:     for (p=0; p<npoints_curr; p++) {

212:       sf_cells[p].rank  = 0;
213:       sf_cells[p].index = p_cellid[p];
214:       if (p_cellid[p] > range) {
215:         range = p_cellid[p];
216:       }
217:     }
218:     DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&p_cellid);
219:     DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);

221:     /* PetscSFCreate(PetscObjectComm((PetscObject)dm),&sfcell); */
222:     PetscSFCreate(PETSC_COMM_SELF,&sfcell);
223:     PetscSFSetGraph(sfcell, range, npoints_curr, NULL, PETSC_OWN_POINTER, sf_cells, PETSC_OWN_POINTER);
224:   }
225: #endif

227:   DMSwarmCreateLocalVectorFromField(dm, DMSwarmPICField_coor, &pos);
228:   DMLocatePoints(dmcell, pos, DM_POINTLOCATION_NONE, &sfcell);
229:   DMSwarmDestroyLocalVectorFromField(dm, DMSwarmPICField_coor, &pos);

231:   if (error_check) {
232:     DMSwarmGetSize(dm,&npointsg);
233:   }
234:   DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
235:   DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
236:   PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell);
237:   for (p=0; p<npoints; p++) {
238:     rankval[p] = LA_sfcell[p].index;
239:   }
240:   DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
241:   PetscSFDestroy(&sfcell);

243:   if (size > 1) {
244:     DMSwarmMigrate_DMNeighborScatter(dm,dmcell,remove_sent_points,&npoints_prior_migration);
245:   } else {
246:     DMSwarmDataField PField;
247:     PetscInt npoints_curr;

249:     /* remove points which the domain */
250:     DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,DMSwarmField_rank,&PField);
251:     DMSwarmDataFieldGetEntries(PField,(void**)&rankval);

253:     DMSwarmDataBucketGetSizes(swarm->db,&npoints_curr,NULL,NULL);
254:     for (p=0; p<npoints_curr; p++) {
255:       if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
256:         /* kill point */
257:         DMSwarmDataBucketRemovePointAtIndex(swarm->db,p);
258:         DMSwarmDataBucketGetSizes(swarm->db,&npoints_curr,NULL,NULL); /* you need to update npoints as the list size decreases! */
259:         DMSwarmDataFieldGetEntries(PField,(void**)&rankval); /* update date point increase realloc performed */
260:         p--; /* check replacement point */
261:       }
262:     }
263:     DMSwarmGetSize(dm,&npoints_prior_migration);

265:   }

267:   /* locate points newly received */
268:   DMSwarmDataBucketGetSizes(swarm->db,&npoints2,NULL,NULL);

270: #if 0
271:   { /* safe alternative - however this performs two point locations on: (i) the initial points set and; (ii) the (initial + received) point set */
272:     PetscScalar *LA_coor;
273:     PetscInt bs;
274:     DMSwarmDataField PField;

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

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

283:     PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell);
284:     DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
285:     for (p=0; p<npoints2; p++) {
286:       rankval[p] = LA_sfcell[p].index;
287:     }
288:     PetscSFDestroy(&sfcell);
289:     DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);

291:     /* remove points which left processor */
292:     DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,DMSwarmField_rank,&PField);
293:     DMSwarmDataFieldGetEntries(PField,(void**)&rankval);

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

308:   { /* perform two point locations: (i) on the initial points set prior to communication; and (ii) on the new (received) points */
309:     PetscScalar      *LA_coor;
310:     PetscInt         npoints_from_neighbours,bs;
311:     DMSwarmDataField PField;

313:     npoints_from_neighbours = npoints2 - npoints_prior_migration;

315:     DMSwarmGetField(dm,DMSwarmPICField_coor,&bs,NULL,(void**)&LA_coor);
316:     VecCreateSeqWithArray(PETSC_COMM_SELF,bs,bs*npoints_from_neighbours,(const PetscScalar*)&LA_coor[bs*npoints_prior_migration],&pos);

318:     DMLocatePoints(dmcell,pos,DM_POINTLOCATION_NONE,&sfcell);

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

323:     PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell);
324:     DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
325:     for (p=0; p<npoints_from_neighbours; p++) {
326:       rankval[npoints_prior_migration + p] = LA_sfcell[p].index;
327:     }
328:     DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
329:     PetscSFDestroy(&sfcell);

331:     /* remove points which left processor */
332:     DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,DMSwarmField_rank,&PField);
333:     DMSwarmDataFieldGetEntries(PField,(void**)&rankval);

335:     DMSwarmDataBucketGetSizes(swarm->db,&npoints2,NULL,NULL);
336:     for (p=npoints_prior_migration; p<npoints2; p++) {
337:       if (rankval[p] == DMLOCATEPOINT_POINT_NOT_FOUND) {
338:         /* kill point */
339:         DMSwarmDataBucketRemovePointAtIndex(swarm->db,p);
340:         DMSwarmDataBucketGetSizes(swarm->db,&npoints2,NULL,NULL); /* you need to update npoints as the list size decreases! */
341:         DMSwarmDataFieldGetEntries(PField,(void**)&rankval); /* update date point increase realloc performed */
342:         p--; /* check replacement point */
343:       }
344:     }
345:   }

347:   {
348:     PetscInt *p_cellid;

350:     DMSwarmDataBucketGetSizes(swarm->db,&npoints2,NULL,NULL);
351:     DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
352:     DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&p_cellid);
353:     for (p=0; p<npoints2; p++) {
354:       p_cellid[p] = rankval[p];
355:     }
356:     DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&p_cellid);
357:     DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
358:   }

360:   /* check for error on removed points */
361:   if (error_check) {
362:     DMSwarmGetSize(dm,&npoints2g);
364:   }
365:   return 0;
366: }

368: PetscErrorCode DMSwarmMigrate_CellDMExact(DM dm,PetscBool remove_sent_points)
369: {
370:   return 0;
371: }

373: /*
374:  Redundant as this assumes points can only be sent to a single rank
375: */
376: PetscErrorCode DMSwarmMigrate_GlobalToLocal_Basic(DM dm,PetscInt *globalsize)
377: {
378:   DM_Swarm       *swarm = (DM_Swarm*)dm->data;
379:   DMSwarmDataEx  de;
380:   PetscInt       p,npoints,*rankval,n_points_recv;
381:   PetscMPIInt    rank,nrank,negrank;
382:   void           *point_buffer,*recv_points;
383:   size_t         sizeof_dmswarm_point;

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

432:     DMSwarmDataBucketInsertPackedArray(swarm->db,npoints+p,data_p);
433:   }
434:   DMSwarmDataExView(de);
435:   DMSwarmDataBucketDestroyPackedArray(swarm->db,&point_buffer);
436:   DMSwarmDataExDestroy(de);
437:   return 0;
438: }

440: typedef struct {
441:   PetscMPIInt owner_rank;
442:   PetscReal   min[3],max[3];
443: } CollectBBox;

445: PETSC_EXTERN PetscErrorCode DMSwarmCollect_DMDABoundingBox(DM dm,PetscInt *globalsize)
446: {
447:   DM_Swarm *        swarm = (DM_Swarm*)dm->data;
448:   DMSwarmDataEx     de;
449:   PetscInt          p,pk,npoints,*rankval,n_points_recv,n_bbox_recv,dim,neighbour_cells;
450:   PetscMPIInt       rank,nrank;
451:   void              *point_buffer,*recv_points;
452:   size_t            sizeof_dmswarm_point,sizeof_bbox_ctx;
453:   PetscBool         isdmda;
454:   CollectBBox       *bbox,*recv_bbox;
455:   const PetscMPIInt *dmneighborranks;
456:   DM                dmcell;

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

460:   DMSwarmGetCellDM(dm,&dmcell);
462:   isdmda = PETSC_FALSE;
463:   PetscObjectTypeCompare((PetscObject)dmcell,DMDA,&isdmda);

466:   DMGetDimension(dm,&dim);
467:   sizeof_bbox_ctx = sizeof(CollectBBox);
468:   PetscMalloc1(1,&bbox);
469:   bbox->owner_rank = rank;

471:   /* compute the bounding box based on the overlapping / stenctil size */
472:   {
473:     Vec lcoor;

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

541:     DMSwarmGetField(dm,"coorx",NULL,NULL,(void**)&array_x);
542:     DMSwarmGetField(dm,"coory",NULL,NULL,(void**)&array_y);
543:     for (p=0; p<npoints; p++) {
544:       if ((array_x[p] >= recv_bbox[pk].min[0]) && (array_x[p] <= recv_bbox[pk].max[0])) {
545:         if ((array_y[p] >= recv_bbox[pk].min[1]) && (array_y[p] <= recv_bbox[pk].max[1])) {
546:           DMSwarmDataExAddToSendCount(de,recv_bbox[pk].owner_rank,1);
547:         }
548:       }
549:     }
550:     DMSwarmRestoreField(dm,"coory",NULL,NULL,(void**)&array_y);
551:     DMSwarmRestoreField(dm,"coorx",NULL,NULL,(void**)&array_x);
552:   }
553:   DMSwarmDataExFinalizeSendCount(de);
554:   DMSwarmDataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);
555:   DMSwarmDataExPackInitialize(de,sizeof_dmswarm_point);
556:   for (pk=0; pk<n_bbox_recv; pk++) {
557:     PetscReal *array_x,*array_y;

559:     DMSwarmGetField(dm,"coorx",NULL,NULL,(void**)&array_x);
560:     DMSwarmGetField(dm,"coory",NULL,NULL,(void**)&array_y);
561:     for (p=0; p<npoints; p++) {
562:       if ((array_x[p] >= recv_bbox[pk].min[0]) && (array_x[p] <= recv_bbox[pk].max[0])) {
563:         if ((array_y[p] >= recv_bbox[pk].min[1]) && (array_y[p] <= recv_bbox[pk].max[1])) {
564:           /* copy point into buffer */
565:           DMSwarmDataBucketFillPackedArray(swarm->db,p,point_buffer);
566:           /* insert point buffer into DMSwarmDataExchanger */
567:           DMSwarmDataExPackData(de,recv_bbox[pk].owner_rank,1,point_buffer);
568:         }
569:       }
570:     }
571:     DMSwarmRestoreField(dm,"coory",NULL,NULL,(void**)&array_y);
572:     DMSwarmRestoreField(dm,"coorx",NULL,NULL,(void**)&array_x);
573:   }
574:   DMSwarmDataExPackFinalize(de);
575:   DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
576:   DMSwarmDataExBegin(de);
577:   DMSwarmDataExEnd(de);
578:   DMSwarmDataExGetRecvData(de,&n_points_recv,(void**)&recv_points);
579:   DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
580:   DMSwarmDataBucketSetSizes(swarm->db,npoints + n_points_recv,DMSWARM_DATA_BUCKET_BUFFER_DEFAULT);
581:   for (p=0; p<n_points_recv; p++) {
582:     void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point);

584:     DMSwarmDataBucketInsertPackedArray(swarm->db,npoints+p,data_p);
585:   }
586:   DMSwarmDataBucketDestroyPackedArray(swarm->db,&point_buffer);
587:   PetscFree(bbox);
588:   DMSwarmDataExView(de);
589:   DMSwarmDataExDestroy(de);
590:   return 0;
591: }

593: /* General collection when no order, or neighbour information is provided */
594: /*
595:  User provides context and collect() method
596:  Broadcast user context

598:  for each context / rank {
599:    collect(swarm,context,n,list)
600:  }
601: */
602: PETSC_EXTERN PetscErrorCode DMSwarmCollect_General(DM dm,PetscErrorCode (*collect)(DM,void*,PetscInt*,PetscInt**),size_t ctx_size,void *ctx,PetscInt *globalsize)
603: {
604:   DM_Swarm       *swarm = (DM_Swarm*)dm->data;
605:   DMSwarmDataEx  de;
606:   PetscInt       p,r,npoints,n_points_recv;
607:   PetscMPIInt    size,rank;
608:   void           *point_buffer,*recv_points;
609:   void           *ctxlist;
610:   PetscInt       *n2collect,**collectlist;
611:   size_t         sizeof_dmswarm_point;

613:   MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);
614:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
615:   DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
616:   *globalsize = npoints;
617:   /* Broadcast user context */
618:   PetscMalloc(ctx_size*size,&ctxlist);
619:   MPI_Allgather(ctx,ctx_size,MPI_CHAR,ctxlist,ctx_size,MPI_CHAR,PetscObjectComm((PetscObject)dm));
620:   PetscMalloc1(size,&n2collect);
621:   PetscMalloc1(size,&collectlist);
622:   for (r=0; r<size; r++) {
623:     PetscInt _n2collect;
624:     PetscInt *_collectlist;
625:     void     *_ctx_r;

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

674:     DMSwarmDataBucketInsertPackedArray(swarm->db,npoints+p,data_p);
675:   }
676:   /* Release memory */
677:   for (r=0; r<size; r++) {
678:     if (collectlist[r]) PetscFree(collectlist[r]);
679:   }
680:   PetscFree(collectlist);
681:   PetscFree(n2collect);
682:   PetscFree(ctxlist);
683:   DMSwarmDataBucketDestroyPackedArray(swarm->db,&point_buffer);
684:   DMSwarmDataExView(de);
685:   DMSwarmDataExDestroy(de);
686:   return 0;
687: }