1df21e3a8SDave May 2df21e3a8SDave May #include <petsc/private/dmswarmimpl.h> /*I "petscdmswarm.h" I*/ 3df21e3a8SDave May #include "data_bucket.h" 4df21e3a8SDave May #include "data_ex.h" 5df21e3a8SDave May 6df21e3a8SDave May 7*480eef7bSDave May /* 8*480eef7bSDave May User loads desired location (MPI rank) into field DMSwarm_rank 9*480eef7bSDave May */ 10df21e3a8SDave May #undef __FUNCT__ 11df21e3a8SDave May #define __FUNCT__ "DMSwarmMigrate_Push_Basic" 12df21e3a8SDave May PetscErrorCode DMSwarmMigrate_Push_Basic(DM dm,PetscBool remove_sent_points) 13df21e3a8SDave May { 14df21e3a8SDave May DM_Swarm *swarm = (DM_Swarm*)dm->data; 15df21e3a8SDave May PetscErrorCode ierr; 16df21e3a8SDave May DataEx de; 17df21e3a8SDave May PetscInt p,npoints,*rankval,n_points_recv; 18df21e3a8SDave May PetscMPIInt rank,nrank; 19df21e3a8SDave May void *point_buffer,*recv_points; 20df21e3a8SDave May size_t sizeof_dmswarm_point; 21df21e3a8SDave May 22df21e3a8SDave May ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr); 23df21e3a8SDave May 24df21e3a8SDave May ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 25df21e3a8SDave May ierr = DMSwarmGetField(dm,"DMSwarm_rank",NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 26df21e3a8SDave May 27df21e3a8SDave May de = DataExCreate(PetscObjectComm((PetscObject)dm),0);CHKERRQ(ierr); 28df21e3a8SDave May 29df21e3a8SDave May ierr = DataExTopologyInitialize(de);CHKERRQ(ierr); 30df21e3a8SDave May for (p=0; p<npoints; p++) { 31df21e3a8SDave May nrank = rankval[p]; 32df21e3a8SDave May if (nrank != rank) { 33df21e3a8SDave May ierr = DataExTopologyAddNeighbour(de,nrank);CHKERRQ(ierr); 34df21e3a8SDave May } 35df21e3a8SDave May } 36df21e3a8SDave May ierr = DataExTopologyFinalize(de);CHKERRQ(ierr); 37df21e3a8SDave May 38df21e3a8SDave May ierr = DataExInitializeSendCount(de);CHKERRQ(ierr); 39df21e3a8SDave May for (p=0; p<npoints; p++) { 40df21e3a8SDave May nrank = rankval[p]; 41df21e3a8SDave May if (nrank != rank) { 42df21e3a8SDave May ierr = DataExAddToSendCount(de,nrank,1);CHKERRQ(ierr); 43df21e3a8SDave May } 44df21e3a8SDave May } 45df21e3a8SDave May ierr = DataExFinalizeSendCount(de);CHKERRQ(ierr); 46df21e3a8SDave May 47df21e3a8SDave May ierr = DataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);CHKERRQ(ierr); 48df21e3a8SDave May ierr = DataExPackInitialize(de,sizeof_dmswarm_point);CHKERRQ(ierr); 49df21e3a8SDave May for (p=0; p<npoints; p++) { 50df21e3a8SDave May nrank = rankval[p]; 51df21e3a8SDave May if (nrank != rank) { 52df21e3a8SDave May /* copy point into buffer */ 53df21e3a8SDave May ierr = DataBucketFillPackedArray(swarm->db,p,point_buffer);CHKERRQ(ierr); 54df21e3a8SDave May /* insert point buffer into DataExchanger */ 55df21e3a8SDave May ierr = DataExPackData(de,nrank,1,point_buffer);CHKERRQ(ierr); 56df21e3a8SDave May } 57df21e3a8SDave May } 58df21e3a8SDave May ierr = DataExPackFinalize(de);CHKERRQ(ierr); 59df21e3a8SDave May 60df21e3a8SDave May 61df21e3a8SDave May if (remove_sent_points) { 62df21e3a8SDave May /* remove points which left processor */ 63df21e3a8SDave May ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 64df21e3a8SDave May for (p=0; p<npoints; p++) { 65df21e3a8SDave May nrank = rankval[p]; 66df21e3a8SDave May if (nrank != rank) { 67df21e3a8SDave May /* kill point */ 68df21e3a8SDave May ierr = DataBucketRemovePointAtIndex(swarm->db,p);CHKERRQ(ierr); 69df21e3a8SDave May DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); /* you need to update npoints as the list size decreases! */ 70df21e3a8SDave May p--; /* check replacement point */ 71df21e3a8SDave May } 72df21e3a8SDave May } 73df21e3a8SDave May } 74df21e3a8SDave May ierr = DMSwarmRestoreField(dm,"DMSwarm_rank",NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 75df21e3a8SDave May 76df21e3a8SDave May ierr = DataExBegin(de);CHKERRQ(ierr); 77df21e3a8SDave May ierr = DataExEnd(de);CHKERRQ(ierr); 78df21e3a8SDave May 79df21e3a8SDave May ierr = DataExGetRecvData(de,&n_points_recv,(void**)&recv_points);CHKERRQ(ierr); 80df21e3a8SDave May 81df21e3a8SDave May ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 82df21e3a8SDave May ierr = DataBucketSetSizes(swarm->db,npoints + n_points_recv,-1);CHKERRQ(ierr); 83df21e3a8SDave May for (p=0; p<n_points_recv; p++) { 84df21e3a8SDave May void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point ); 85df21e3a8SDave May 86df21e3a8SDave May ierr = DataBucketInsertPackedArray(swarm->db,npoints+p,data_p);CHKERRQ(ierr); 87df21e3a8SDave May } 88df21e3a8SDave May ierr = DataExView(de);CHKERRQ(ierr); 89df21e3a8SDave May ierr = DataBucketDestroyPackedArray(swarm->db,&point_buffer);CHKERRQ(ierr); 90df21e3a8SDave May ierr = DataExDestroy(de);CHKERRQ(ierr); 91df21e3a8SDave May 92df21e3a8SDave May PetscFunctionReturn(0); 93df21e3a8SDave May } 942712d1f2SDave May 95*480eef7bSDave May //DMLocatePoints 96*480eef7bSDave May 97*480eef7bSDave May #undef __FUNCT__ 98*480eef7bSDave May #define __FUNCT__ "DMSwarmMigrate_CellDM" 99*480eef7bSDave May PetscErrorCode DMSwarmMigrate_CellDM(DM dm,PetscBool remove_sent_points) 100*480eef7bSDave May { 101*480eef7bSDave May DM_Swarm *swarm = (DM_Swarm*)dm->data; 102*480eef7bSDave May PetscErrorCode ierr; 103*480eef7bSDave May PetscInt p,npoints,*rankval; 104*480eef7bSDave May const PetscInt *LA_iscell; 105*480eef7bSDave May DM dmcell; 106*480eef7bSDave May IS iscell; 107*480eef7bSDave May Vec pos; 108*480eef7bSDave May 109*480eef7bSDave May ierr = DMSwarmGetCellDM(dm,&dmcell);CHKERRQ(ierr); 110*480eef7bSDave May if (!dmcell) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only valid if cell DM provided"); 111*480eef7bSDave May 112*480eef7bSDave May ierr = DMSwarmCreateGlobalVectorFromField(dm,"coor",&pos);CHKERRQ(ierr); 113*480eef7bSDave May ierr = DMLocatePoints(dmcell,pos,&iscell);CHKERRQ(ierr); 114*480eef7bSDave May ierr = DMSwarmDestroyGlobalVectorFromField(dm,"coor",&pos);CHKERRQ(ierr); 115*480eef7bSDave May 116*480eef7bSDave May ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 117*480eef7bSDave May ierr = DMSwarmGetField(dm,"DMSwarm_rank",NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 118*480eef7bSDave May ierr = ISGetIndices(iscell,&LA_iscell);CHKERRQ(ierr); 119*480eef7bSDave May for (p=0; p<npoints; p++) { 120*480eef7bSDave May if (LA_iscell[p] == -1) { 121*480eef7bSDave May rankval[p] = -1; 122*480eef7bSDave May } 123*480eef7bSDave May } 124*480eef7bSDave May ierr = ISRestoreIndices(iscell,&LA_iscell);CHKERRQ(ierr); 125*480eef7bSDave May ierr = DMSwarmRestoreField(dm,"DMSwarm_rank",NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 126*480eef7bSDave May ierr = ISDestroy(&iscell);CHKERRQ(ierr); 127*480eef7bSDave May 128*480eef7bSDave May ierr = DMSwarmMigrate_Push_Basic(dm,remove_sent_points);CHKERRQ(ierr); 129*480eef7bSDave May 130*480eef7bSDave May PetscFunctionReturn(0); 131*480eef7bSDave May } 132*480eef7bSDave May 133*480eef7bSDave May /* 134*480eef7bSDave May Redundant as this assumes points can only be sent to a single rank 135*480eef7bSDave May */ 1362712d1f2SDave May #undef __FUNCT__ 1372712d1f2SDave May #define __FUNCT__ "DMSwarmMigrate_GlobalToLocal_Basic" 1382712d1f2SDave May PetscErrorCode DMSwarmMigrate_GlobalToLocal_Basic(DM dm,PetscInt *globalsize) 1392712d1f2SDave May { 1402712d1f2SDave May DM_Swarm *swarm = (DM_Swarm*)dm->data; 1412712d1f2SDave May PetscErrorCode ierr; 1422712d1f2SDave May DataEx de; 1432712d1f2SDave May PetscInt p,npoints,*rankval,n_points_recv; 1442712d1f2SDave May PetscMPIInt rank,nrank,negrank; 1452712d1f2SDave May void *point_buffer,*recv_points; 1462712d1f2SDave May size_t sizeof_dmswarm_point; 1472712d1f2SDave May 1482712d1f2SDave May ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr); 1492712d1f2SDave May 1502712d1f2SDave May ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 1512712d1f2SDave May *globalsize = npoints; 1522712d1f2SDave May ierr = DMSwarmGetField(dm,"DMSwarm_rank",NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 1532712d1f2SDave May 1542712d1f2SDave May de = DataExCreate(PetscObjectComm((PetscObject)dm),0);CHKERRQ(ierr); 1552712d1f2SDave May 1562712d1f2SDave May ierr = DataExTopologyInitialize(de);CHKERRQ(ierr); 1572712d1f2SDave May for (p=0; p<npoints; p++) { 1582712d1f2SDave May negrank = rankval[p]; 1592712d1f2SDave May if (negrank < 0) { 1602712d1f2SDave May nrank = -negrank - 1; 1612712d1f2SDave May ierr = DataExTopologyAddNeighbour(de,nrank);CHKERRQ(ierr); 1622712d1f2SDave May } 1632712d1f2SDave May } 1642712d1f2SDave May ierr = DataExTopologyFinalize(de);CHKERRQ(ierr); 1652712d1f2SDave May 1662712d1f2SDave May ierr = DataExInitializeSendCount(de);CHKERRQ(ierr); 1672712d1f2SDave May for (p=0; p<npoints; p++) { 1682712d1f2SDave May negrank = rankval[p]; 1692712d1f2SDave May if (negrank < 0) { 1702712d1f2SDave May nrank = -negrank - 1; 1712712d1f2SDave May ierr = DataExAddToSendCount(de,nrank,1);CHKERRQ(ierr); 1722712d1f2SDave May } 1732712d1f2SDave May } 1742712d1f2SDave May ierr = DataExFinalizeSendCount(de);CHKERRQ(ierr); 1752712d1f2SDave May 1762712d1f2SDave May ierr = DataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);CHKERRQ(ierr); 1772712d1f2SDave May ierr = DataExPackInitialize(de,sizeof_dmswarm_point);CHKERRQ(ierr); 1782712d1f2SDave May for (p=0; p<npoints; p++) { 1792712d1f2SDave May negrank = rankval[p]; 1802712d1f2SDave May if (negrank < 0) { 1812712d1f2SDave May nrank = -negrank - 1; 1822712d1f2SDave May rankval[p] = nrank; 1832712d1f2SDave May /* copy point into buffer */ 1842712d1f2SDave May ierr = DataBucketFillPackedArray(swarm->db,p,point_buffer);CHKERRQ(ierr); 1852712d1f2SDave May /* insert point buffer into DataExchanger */ 1862712d1f2SDave May ierr = DataExPackData(de,nrank,1,point_buffer);CHKERRQ(ierr); 1872712d1f2SDave May rankval[p] = negrank; 1882712d1f2SDave May } 1892712d1f2SDave May } 1902712d1f2SDave May ierr = DataExPackFinalize(de);CHKERRQ(ierr); 1912712d1f2SDave May 1922712d1f2SDave May ierr = DMSwarmRestoreField(dm,"DMSwarm_rank",NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 1932712d1f2SDave May 1942712d1f2SDave May ierr = DataExBegin(de);CHKERRQ(ierr); 1952712d1f2SDave May ierr = DataExEnd(de);CHKERRQ(ierr); 1962712d1f2SDave May 1972712d1f2SDave May ierr = DataExGetRecvData(de,&n_points_recv,(void**)&recv_points);CHKERRQ(ierr); 1982712d1f2SDave May 1992712d1f2SDave May ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 2002712d1f2SDave May ierr = DataBucketSetSizes(swarm->db,npoints + n_points_recv,-1);CHKERRQ(ierr); 2012712d1f2SDave May for (p=0; p<n_points_recv; p++) { 2022712d1f2SDave May void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point ); 2032712d1f2SDave May 2042712d1f2SDave May ierr = DataBucketInsertPackedArray(swarm->db,npoints+p,data_p);CHKERRQ(ierr); 2052712d1f2SDave May } 2062712d1f2SDave May ierr = DataExView(de);CHKERRQ(ierr); 2072712d1f2SDave May ierr = DataBucketDestroyPackedArray(swarm->db,&point_buffer);CHKERRQ(ierr); 2082712d1f2SDave May ierr = DataExDestroy(de);CHKERRQ(ierr); 2092712d1f2SDave May 2102712d1f2SDave May PetscFunctionReturn(0); 2112712d1f2SDave May } 212b16650c8SDave May 213b16650c8SDave May typedef struct { 214b16650c8SDave May PetscMPIInt owner_rank; 215b16650c8SDave May PetscReal min[3],max[3]; 216b16650c8SDave May } CollectBBox; 217b16650c8SDave May 218b16650c8SDave May #undef __FUNCT__ 219fe39f135SDave May #define __FUNCT__ "DMSwarmCollect_DMDABoundingBox" 220fe39f135SDave May PETSC_EXTERN PetscErrorCode DMSwarmCollect_DMDABoundingBox(DM dm,PetscInt *globalsize) 221b16650c8SDave May { 222b16650c8SDave May DM_Swarm *swarm = (DM_Swarm*)dm->data; 223b16650c8SDave May PetscErrorCode ierr; 224b16650c8SDave May DataEx de; 225b16650c8SDave May PetscInt p,pk,npoints,*rankval,n_points_recv,n_bbox_recv,dim,neighbour_cells; 226b16650c8SDave May PetscMPIInt rank,nrank; 227b16650c8SDave May void *point_buffer,*recv_points; 228b16650c8SDave May size_t sizeof_dmswarm_point,sizeof_bbox_ctx; 229b16650c8SDave May PetscBool isdmda; 230b16650c8SDave May CollectBBox *bbox,*recv_bbox; 231b16650c8SDave May const PetscMPIInt *dmneighborranks; 232b16650c8SDave May DM dmcell; 233b16650c8SDave May 234b16650c8SDave May ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr); 235b16650c8SDave May 236fe39f135SDave May ierr = DMSwarmGetCellDM(dm,&dmcell);CHKERRQ(ierr); 237fe39f135SDave May if (!dmcell) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only valid if cell DM provided"); 238b16650c8SDave May 239b16650c8SDave May isdmda = PETSC_FALSE; 240b16650c8SDave May PetscObjectTypeCompare((PetscObject)dmcell,DMDA,&isdmda); 241b16650c8SDave May if (!isdmda) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only DMDA support for CollectBoundingBox"); 242b16650c8SDave May 243fe39f135SDave May ierr = DMDAGetInfo(dm,&dim,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);CHKERRQ(ierr); 244fe39f135SDave May if (dim == 1) { 245fe39f135SDave May neighbour_cells = 3; 246fe39f135SDave May } else if (dim == 2) { 247fe39f135SDave May neighbour_cells = 9; 248fe39f135SDave May } else { 249fe39f135SDave May neighbour_cells = 27; 250fe39f135SDave May } 251fe39f135SDave May 252b16650c8SDave May sizeof_bbox_ctx = sizeof(CollectBBox); 253b16650c8SDave May PetscMalloc1(1,&bbox); 254b16650c8SDave May bbox->owner_rank = rank; 255b16650c8SDave May 256b16650c8SDave May /* compute the bounding box based on the overlapping / stenctil size */ 257b16650c8SDave May //ierr = DMDAGetLocalBoundingBox(dmcell,bbox->min,bbox->max);CHKERRQ(ierr); 258b16650c8SDave May { 259b16650c8SDave May Vec lcoor; 260b16650c8SDave May 261b16650c8SDave May ierr = DMGetCoordinatesLocal(dmcell,&lcoor);CHKERRQ(ierr); 262b16650c8SDave May 263fe39f135SDave May if (dim >= 1) { 264b16650c8SDave May ierr = VecStrideMin(lcoor,0,NULL,&bbox->min[0]);CHKERRQ(ierr); 265b16650c8SDave May ierr = VecStrideMax(lcoor,0,NULL,&bbox->max[0]);CHKERRQ(ierr); 266fe39f135SDave May } 267fe39f135SDave May 268fe39f135SDave May if (dim >= 2) { 269fe39f135SDave May ierr = VecStrideMin(lcoor,1,NULL,&bbox->min[1]);CHKERRQ(ierr); 270b16650c8SDave May ierr = VecStrideMax(lcoor,1,NULL,&bbox->max[1]);CHKERRQ(ierr); 271b16650c8SDave May } 272b16650c8SDave May 273fe39f135SDave May if (dim == 3) { 274fe39f135SDave May ierr = VecStrideMin(lcoor,2,NULL,&bbox->min[2]);CHKERRQ(ierr); 275fe39f135SDave May ierr = VecStrideMax(lcoor,2,NULL,&bbox->max[2]);CHKERRQ(ierr); 276fe39f135SDave May } 277fe39f135SDave May } 278fe39f135SDave May 279b16650c8SDave May ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 280b16650c8SDave May *globalsize = npoints; 281b16650c8SDave May ierr = DMSwarmGetField(dm,"DMSwarm_rank",NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 282b16650c8SDave May 283b16650c8SDave May de = DataExCreate(PetscObjectComm((PetscObject)dm),0);CHKERRQ(ierr); 284b16650c8SDave May 285b16650c8SDave May /* use DMDA neighbours */ 286b16650c8SDave May ierr = DMDAGetNeighbors(dmcell,&dmneighborranks);CHKERRQ(ierr); 287b16650c8SDave May 288b16650c8SDave May ierr = DataExTopologyInitialize(de);CHKERRQ(ierr); 289b16650c8SDave May for (p=0; p<neighbour_cells; p++) { 290b16650c8SDave May if ( (dmneighborranks[p] >= 0) && (dmneighborranks[p] != rank) ) { 291b16650c8SDave May ierr = DataExTopologyAddNeighbour(de,dmneighborranks[p]);CHKERRQ(ierr); 292b16650c8SDave May } 293b16650c8SDave May } 294b16650c8SDave May ierr = DataExTopologyFinalize(de);CHKERRQ(ierr); 295b16650c8SDave May 296b16650c8SDave May ierr = DataExInitializeSendCount(de);CHKERRQ(ierr); 297b16650c8SDave May for (p=0; p<neighbour_cells; p++) { 298b16650c8SDave May if ( (dmneighborranks[p] >= 0) && (dmneighborranks[p] != rank) ) { 299b16650c8SDave May ierr = DataExAddToSendCount(de,dmneighborranks[p],1);CHKERRQ(ierr); 300b16650c8SDave May } 301b16650c8SDave May } 302b16650c8SDave May ierr = DataExFinalizeSendCount(de);CHKERRQ(ierr); 303b16650c8SDave May 304b16650c8SDave May 305b16650c8SDave May /* send bounding boxes */ 306b16650c8SDave May ierr = DataExPackInitialize(de,sizeof_bbox_ctx);CHKERRQ(ierr); 307b16650c8SDave May for (p=0; p<neighbour_cells; p++) { 308b16650c8SDave May nrank = dmneighborranks[p]; 309b16650c8SDave May if ( (nrank >= 0) && (nrank != rank) ) { 310b16650c8SDave May /* insert bbox buffer into DataExchanger */ 311b16650c8SDave May ierr = DataExPackData(de,nrank,1,bbox);CHKERRQ(ierr); 312b16650c8SDave May } 313b16650c8SDave May } 314b16650c8SDave May ierr = DataExPackFinalize(de);CHKERRQ(ierr); 315b16650c8SDave May 316b16650c8SDave May /* recv bounding boxes */ 317b16650c8SDave May ierr = DataExBegin(de);CHKERRQ(ierr); 318b16650c8SDave May ierr = DataExEnd(de);CHKERRQ(ierr); 319b16650c8SDave May 320b16650c8SDave May ierr = DataExGetRecvData(de,&n_bbox_recv,(void**)&recv_bbox);CHKERRQ(ierr); 321b16650c8SDave May 322b16650c8SDave May 323b16650c8SDave May for (p=0; p<n_bbox_recv; p++) { 324b16650c8SDave May printf("[rank %d]: box from %d : range[%+1.4e,%+1.4e]x[%+1.4e,%+1.4e]\n",rank,recv_bbox[p].owner_rank, 325b16650c8SDave May recv_bbox[p].min[0],recv_bbox[p].max[0],recv_bbox[p].min[1],recv_bbox[p].max[1]); 326b16650c8SDave May } 327b16650c8SDave May 328b16650c8SDave May /* of course this is stupid as this "generic" function should have a better way to know what the coordinates are called */ 329b16650c8SDave May ierr = DataExInitializeSendCount(de);CHKERRQ(ierr); 330b16650c8SDave May for (pk=0; pk<n_bbox_recv; pk++) { 331b16650c8SDave May PetscReal *array_x,*array_y; 332b16650c8SDave May 333b16650c8SDave May ierr = DMSwarmGetField(dm,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr); 334b16650c8SDave May ierr = DMSwarmGetField(dm,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr); 335b16650c8SDave May 336b16650c8SDave May for (p=0; p<npoints; p++) { 337b16650c8SDave May if ((array_x[p] >= recv_bbox[pk].min[0]) && (array_x[p] <= recv_bbox[pk].max[0]) ) { 338b16650c8SDave May if ((array_y[p] >= recv_bbox[pk].min[1]) && (array_y[p] <= recv_bbox[pk].max[1]) ) { 339b16650c8SDave May 340b16650c8SDave May ierr = DataExAddToSendCount(de,recv_bbox[pk].owner_rank,1);CHKERRQ(ierr); 341b16650c8SDave May 342b16650c8SDave May } 343b16650c8SDave May } 344b16650c8SDave May } 345b16650c8SDave May ierr = DMSwarmRestoreField(dm,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr); 346b16650c8SDave May ierr = DMSwarmRestoreField(dm,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr); 347b16650c8SDave May } 348b16650c8SDave May ierr = DataExFinalizeSendCount(de);CHKERRQ(ierr); 349b16650c8SDave May 350b16650c8SDave May 351b16650c8SDave May 352b16650c8SDave May ierr = DataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);CHKERRQ(ierr); 353b16650c8SDave May ierr = DataExPackInitialize(de,sizeof_dmswarm_point);CHKERRQ(ierr); 354b16650c8SDave May for (pk=0; pk<n_bbox_recv; pk++) { 355b16650c8SDave May PetscReal *array_x,*array_y; 356b16650c8SDave May 357b16650c8SDave May ierr = DMSwarmGetField(dm,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr); 358b16650c8SDave May ierr = DMSwarmGetField(dm,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr); 359b16650c8SDave May 360b16650c8SDave May for (p=0; p<npoints; p++) { 361b16650c8SDave May if ((array_x[p] >= recv_bbox[pk].min[0]) && (array_x[p] <= recv_bbox[pk].max[0]) ) { 362b16650c8SDave May if ((array_y[p] >= recv_bbox[pk].min[1]) && (array_y[p] <= recv_bbox[pk].max[1]) ) { 363b16650c8SDave May // copy point into buffer // 364b16650c8SDave May ierr = DataBucketFillPackedArray(swarm->db,p,point_buffer);CHKERRQ(ierr); 365b16650c8SDave May // insert point buffer into DataExchanger // 366b16650c8SDave May ierr = DataExPackData(de,recv_bbox[pk].owner_rank,1,point_buffer);CHKERRQ(ierr); 367b16650c8SDave May } 368b16650c8SDave May } 369b16650c8SDave May } 370b16650c8SDave May 371b16650c8SDave May ierr = DMSwarmRestoreField(dm,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr); 372b16650c8SDave May ierr = DMSwarmRestoreField(dm,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr); 373b16650c8SDave May } 374b16650c8SDave May 375b16650c8SDave May ierr = DataExPackFinalize(de);CHKERRQ(ierr); 376b16650c8SDave May 377b16650c8SDave May ierr = DMSwarmRestoreField(dm,"DMSwarm_rank",NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 378b16650c8SDave May 379b16650c8SDave May ierr = DataExBegin(de);CHKERRQ(ierr); 380b16650c8SDave May ierr = DataExEnd(de);CHKERRQ(ierr); 381b16650c8SDave May 382b16650c8SDave May ierr = DataExGetRecvData(de,&n_points_recv,(void**)&recv_points);CHKERRQ(ierr); 383b16650c8SDave May 384b16650c8SDave May ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 385b16650c8SDave May ierr = DataBucketSetSizes(swarm->db,npoints + n_points_recv,-1);CHKERRQ(ierr); 386b16650c8SDave May for (p=0; p<n_points_recv; p++) { 387b16650c8SDave May void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point ); 388b16650c8SDave May 389b16650c8SDave May ierr = DataBucketInsertPackedArray(swarm->db,npoints+p,data_p);CHKERRQ(ierr); 390b16650c8SDave May } 391b16650c8SDave May 392b16650c8SDave May ierr = DataBucketDestroyPackedArray(swarm->db,&point_buffer);CHKERRQ(ierr); 393b16650c8SDave May PetscFree(bbox); 394b16650c8SDave May ierr = DataExView(de);CHKERRQ(ierr); 395b16650c8SDave May ierr = DataExDestroy(de);CHKERRQ(ierr); 396b16650c8SDave May 397b16650c8SDave May PetscFunctionReturn(0); 398b16650c8SDave May } 399a9fd7477SDave May 400a9fd7477SDave May 401a9fd7477SDave May /* General collection when no order, or neighbour information is provided */ 402a9fd7477SDave May /* 403a9fd7477SDave May User provides context and collect() method 404a9fd7477SDave May Broadcast user context 405a9fd7477SDave May 406a9fd7477SDave May for each context / rank { 407a9fd7477SDave May collect(swarm,context,n,list) 408a9fd7477SDave May } 409a9fd7477SDave May 410a9fd7477SDave May */ 411a9fd7477SDave May #undef __FUNCT__ 412a9fd7477SDave May #define __FUNCT__ "DMSwarmCollect_General" 413a9fd7477SDave May PETSC_EXTERN PetscErrorCode DMSwarmCollect_General(DM dm,PetscErrorCode (*collect)(DM,void*,PetscInt*,PetscInt**),size_t ctx_size,void *ctx,PetscInt *globalsize) 414a9fd7477SDave May { 415a9fd7477SDave May DM_Swarm *swarm = (DM_Swarm*)dm->data; 416a9fd7477SDave May PetscErrorCode ierr; 417a9fd7477SDave May DataEx de; 418a9fd7477SDave May PetscInt p,r,npoints,n_points_recv; 419a9fd7477SDave May PetscMPIInt commsize,rank; 420a9fd7477SDave May void *point_buffer,*recv_points; 421a9fd7477SDave May void *ctxlist; 422a9fd7477SDave May PetscInt *n2collect,**collectlist; 423a9fd7477SDave May size_t sizeof_dmswarm_point; 424a9fd7477SDave May 425a9fd7477SDave May ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm),&commsize);CHKERRQ(ierr); 426a9fd7477SDave May ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr); 427a9fd7477SDave May 428a9fd7477SDave May ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 429a9fd7477SDave May *globalsize = npoints; 430a9fd7477SDave May 431a9fd7477SDave May /* Broadcast user context */ 432a9fd7477SDave May PetscMalloc(ctx_size*commsize,&ctxlist); 433a9fd7477SDave May ierr = MPI_Allgather(ctx,ctx_size,MPI_CHAR,ctxlist,ctx_size,MPI_CHAR,PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 434a9fd7477SDave May 435a9fd7477SDave May PetscMalloc1(commsize,&n2collect); 436a9fd7477SDave May PetscMalloc1(commsize,&collectlist); 437a9fd7477SDave May 438a9fd7477SDave May for (r=0; r<commsize; r++) { 439a9fd7477SDave May PetscInt _n2collect; 440a9fd7477SDave May PetscInt *_collectlist; 441a9fd7477SDave May void *_ctx_r; 442a9fd7477SDave May 443a9fd7477SDave May _n2collect = 0; 444a9fd7477SDave May _collectlist = NULL; 445a9fd7477SDave May 446a9fd7477SDave May if (r != rank) { /* don't collect data from yourself */ 447a9fd7477SDave May _ctx_r = (void*)( (char*)ctxlist + r * ctx_size ); 448a9fd7477SDave May ierr = collect(dm,_ctx_r,&_n2collect,&_collectlist);CHKERRQ(ierr); 449a9fd7477SDave May } 450a9fd7477SDave May 451a9fd7477SDave May n2collect[r] = _n2collect; 452a9fd7477SDave May collectlist[r] = _collectlist; 453a9fd7477SDave May } 454a9fd7477SDave May 455a9fd7477SDave May de = DataExCreate(PetscObjectComm((PetscObject)dm),0);CHKERRQ(ierr); 456a9fd7477SDave May 457a9fd7477SDave May /* Define topology */ 458a9fd7477SDave May ierr = DataExTopologyInitialize(de);CHKERRQ(ierr); 459a9fd7477SDave May for (r=0; r<commsize; r++) { 460a9fd7477SDave May if (n2collect[r] > 0) { 461a9fd7477SDave May ierr = DataExTopologyAddNeighbour(de,(PetscMPIInt)r);CHKERRQ(ierr); 462a9fd7477SDave May } 463a9fd7477SDave May } 464a9fd7477SDave May ierr = DataExTopologyFinalize(de);CHKERRQ(ierr); 465a9fd7477SDave May 466a9fd7477SDave May /* Define send counts */ 467a9fd7477SDave May ierr = DataExInitializeSendCount(de);CHKERRQ(ierr); 468a9fd7477SDave May for (r=0; r<commsize; r++) { 469a9fd7477SDave May if (n2collect[r] > 0) { 470a9fd7477SDave May ierr = DataExAddToSendCount(de,r,n2collect[r]);CHKERRQ(ierr); 471a9fd7477SDave May } 472a9fd7477SDave May } 473a9fd7477SDave May ierr = DataExFinalizeSendCount(de);CHKERRQ(ierr); 474a9fd7477SDave May 475a9fd7477SDave May /* Pack data */ 476a9fd7477SDave May ierr = DataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);CHKERRQ(ierr); 477a9fd7477SDave May 478a9fd7477SDave May ierr = DataExPackInitialize(de,sizeof_dmswarm_point);CHKERRQ(ierr); 479a9fd7477SDave May for (r=0; r<commsize; r++) { 480a9fd7477SDave May 481a9fd7477SDave May for (p=0; p<n2collect[r]; p++) { 482a9fd7477SDave May ierr = DataBucketFillPackedArray(swarm->db,collectlist[r][p],point_buffer);CHKERRQ(ierr); 483a9fd7477SDave May /* insert point buffer into the data exchanger */ 484a9fd7477SDave May ierr = DataExPackData(de,r,1,point_buffer);CHKERRQ(ierr); 485a9fd7477SDave May } 486a9fd7477SDave May } 487a9fd7477SDave May 488a9fd7477SDave May ierr = DataExPackFinalize(de);CHKERRQ(ierr); 489a9fd7477SDave May 490a9fd7477SDave May /* Scatter */ 491a9fd7477SDave May ierr = DataExBegin(de);CHKERRQ(ierr); 492a9fd7477SDave May ierr = DataExEnd(de);CHKERRQ(ierr); 493a9fd7477SDave May 494a9fd7477SDave May /* Collect data in DMSwarm container */ 495a9fd7477SDave May ierr = DataExGetRecvData(de,&n_points_recv,(void**)&recv_points);CHKERRQ(ierr); 496a9fd7477SDave May 497a9fd7477SDave May ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 498a9fd7477SDave May ierr = DataBucketSetSizes(swarm->db,npoints + n_points_recv,-1);CHKERRQ(ierr); 499a9fd7477SDave May for (p=0; p<n_points_recv; p++) { 500a9fd7477SDave May void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point ); 501a9fd7477SDave May 502a9fd7477SDave May ierr = DataBucketInsertPackedArray(swarm->db,npoints+p,data_p);CHKERRQ(ierr); 503a9fd7477SDave May } 504a9fd7477SDave May 505a9fd7477SDave May /* Release memory */ 506a9fd7477SDave May for (r=0; r<commsize; r++) { 507a9fd7477SDave May if (collectlist[r]) PetscFree(collectlist[r]); 508a9fd7477SDave May } 509a9fd7477SDave May PetscFree(collectlist); 510a9fd7477SDave May PetscFree(n2collect); 511a9fd7477SDave May PetscFree(ctxlist); 512a9fd7477SDave May ierr = DataBucketDestroyPackedArray(swarm->db,&point_buffer);CHKERRQ(ierr); 513a9fd7477SDave May ierr = DataExView(de);CHKERRQ(ierr); 514a9fd7477SDave May ierr = DataExDestroy(de);CHKERRQ(ierr); 515a9fd7477SDave May 516a9fd7477SDave May PetscFunctionReturn(0); 517a9fd7477SDave May } 518a9fd7477SDave May 519