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;
     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);

     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;

    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

    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;

    255:     /* remove points which the domain */
    256:     DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,DMSwarmField_rank,&PField);
    257:     DMSwarmDataFieldGetEntries(PField,(void**)&rankval);

    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);

    271:   }

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

    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;

    319:     npoints_from_neighbours = npoints2 - npoints_prior_migration;

    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);

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

    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);

    337:     /* remove points which left processor */
    338:     DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,DMSwarmField_rank,&PField);
    339:     DMSwarmDataFieldGetEntries(PField,(void**)&rankval);

    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:   }

    353:   {
    354:     PetscInt *p_cellid;

    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:   }

    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: }