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;

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

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

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

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

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

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

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

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

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

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

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

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

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

192:   DMSwarmGetCellDM(dm,&dmcell);

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

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

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

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

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

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

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

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

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

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

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

264:   }

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

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

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

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

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

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

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

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

312:     npoints_from_neighbours = npoints2 - npoints_prior_migration;

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

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

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

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

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

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

346:   {
347:     PetscInt *p_cellid;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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