1 2 #include <petsc/private/dmswarmimpl.h> /*I "petscdmswarm.h" I*/ 3 #include "data_bucket.h" 4 #include "data_ex.h" 5 6 7 /* 8 User loads desired location (MPI rank) into field DMSwarm_rank 9 */ 10 #undef __FUNCT__ 11 #define __FUNCT__ "DMSwarmMigrate_Push_Basic" 12 PetscErrorCode DMSwarmMigrate_Push_Basic(DM dm,PetscBool remove_sent_points) 13 { 14 DM_Swarm *swarm = (DM_Swarm*)dm->data; 15 PetscErrorCode ierr; 16 DataEx de; 17 PetscInt p,npoints,*rankval,n_points_recv; 18 PetscMPIInt rank,nrank; 19 void *point_buffer,*recv_points; 20 size_t sizeof_dmswarm_point; 21 22 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr); 23 24 ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 25 ierr = DMSwarmGetField(dm,"DMSwarm_rank",NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 26 27 de = DataExCreate(PetscObjectComm((PetscObject)dm),0);CHKERRQ(ierr); 28 29 ierr = DataExTopologyInitialize(de);CHKERRQ(ierr); 30 for (p=0; p<npoints; p++) { 31 nrank = rankval[p]; 32 if (nrank != rank) { 33 ierr = DataExTopologyAddNeighbour(de,nrank);CHKERRQ(ierr); 34 } 35 } 36 ierr = DataExTopologyFinalize(de);CHKERRQ(ierr); 37 38 ierr = DataExInitializeSendCount(de);CHKERRQ(ierr); 39 for (p=0; p<npoints; p++) { 40 nrank = rankval[p]; 41 if (nrank != rank) { 42 ierr = DataExAddToSendCount(de,nrank,1);CHKERRQ(ierr); 43 } 44 } 45 ierr = DataExFinalizeSendCount(de);CHKERRQ(ierr); 46 47 ierr = DataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);CHKERRQ(ierr); 48 ierr = DataExPackInitialize(de,sizeof_dmswarm_point);CHKERRQ(ierr); 49 for (p=0; p<npoints; p++) { 50 nrank = rankval[p]; 51 if (nrank != rank) { 52 /* copy point into buffer */ 53 ierr = DataBucketFillPackedArray(swarm->db,p,point_buffer);CHKERRQ(ierr); 54 /* insert point buffer into DataExchanger */ 55 ierr = DataExPackData(de,nrank,1,point_buffer);CHKERRQ(ierr); 56 } 57 } 58 ierr = DataExPackFinalize(de);CHKERRQ(ierr); 59 60 61 if (remove_sent_points) { 62 /* remove points which left processor */ 63 ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 64 for (p=0; p<npoints; p++) { 65 nrank = rankval[p]; 66 if (nrank != rank) { 67 /* kill point */ 68 ierr = DataBucketRemovePointAtIndex(swarm->db,p);CHKERRQ(ierr); 69 DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); /* you need to update npoints as the list size decreases! */ 70 p--; /* check replacement point */ 71 } 72 } 73 } 74 ierr = DMSwarmRestoreField(dm,"DMSwarm_rank",NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 75 76 ierr = DataExBegin(de);CHKERRQ(ierr); 77 ierr = DataExEnd(de);CHKERRQ(ierr); 78 79 ierr = DataExGetRecvData(de,&n_points_recv,(void**)&recv_points);CHKERRQ(ierr); 80 81 ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 82 ierr = DataBucketSetSizes(swarm->db,npoints + n_points_recv,-1);CHKERRQ(ierr); 83 for (p=0; p<n_points_recv; p++) { 84 void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point ); 85 86 ierr = DataBucketInsertPackedArray(swarm->db,npoints+p,data_p);CHKERRQ(ierr); 87 } 88 ierr = DataExView(de);CHKERRQ(ierr); 89 ierr = DataBucketDestroyPackedArray(swarm->db,&point_buffer);CHKERRQ(ierr); 90 ierr = DataExDestroy(de);CHKERRQ(ierr); 91 92 PetscFunctionReturn(0); 93 } 94 95 #undef __FUNCT__ 96 #define __FUNCT__ "DMSwarmMigrate_DMNeighborScatter" 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 PetscErrorCode ierr; 101 DataEx 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 109 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr); 110 111 ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 112 ierr = DMSwarmGetField(dm,"DMSwarm_rank",NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 113 114 de = DataExCreate(PetscObjectComm((PetscObject)dm),0);CHKERRQ(ierr); 115 116 ierr = DMGetNeighbors(dmcell,&nneighbors,&neighbourranks);CHKERRQ(ierr); 117 ierr = DataExTopologyInitialize(de);CHKERRQ(ierr); 118 for (r=0; r<nneighbors; r++) { 119 _rank = neighbourranks[r]; 120 if ((_rank != rank) && (_rank > 0)) { 121 ierr = DataExTopologyAddNeighbour(de,_rank);CHKERRQ(ierr); 122 } 123 } 124 ierr = DataExTopologyFinalize(de);CHKERRQ(ierr); 125 126 ierr = DataExInitializeSendCount(de);CHKERRQ(ierr); 127 for (p=0; p<npoints; p++) { 128 if (rankval[p] == -1) { 129 for (r=0; r<nneighbors; r++) { 130 _rank = neighbourranks[r]; 131 if ((_rank != rank) && (_rank > 0)) { 132 ierr = DataExAddToSendCount(de,_rank,1);CHKERRQ(ierr); 133 } 134 } 135 } 136 } 137 ierr = DataExFinalizeSendCount(de);CHKERRQ(ierr); 138 139 ierr = DataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);CHKERRQ(ierr); 140 ierr = DataExPackInitialize(de,sizeof_dmswarm_point);CHKERRQ(ierr); 141 for (p=0; p<npoints; p++) { 142 if (rankval[p] == -1) { 143 for (r=0; r<nneighbors; r++) { 144 _rank = neighbourranks[r]; 145 if ((_rank != rank) && (_rank > 0)) { 146 /* copy point into buffer */ 147 ierr = DataBucketFillPackedArray(swarm->db,p,point_buffer);CHKERRQ(ierr); 148 /* insert point buffer into DataExchanger */ 149 ierr = DataExPackData(de,_rank,1,point_buffer);CHKERRQ(ierr); 150 } 151 } 152 } 153 } 154 ierr = DataExPackFinalize(de);CHKERRQ(ierr); 155 156 157 if (remove_sent_points) { 158 /* remove points which left processor */ 159 ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 160 for (p=0; p<npoints; p++) { 161 if (rankval[p] == -1) { 162 /* kill point */ 163 ierr = DataBucketRemovePointAtIndex(swarm->db,p);CHKERRQ(ierr); 164 DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); /* you need to update npoints as the list size decreases! */ 165 p--; /* check replacement point */ 166 } 167 } 168 } 169 ierr = DMSwarmRestoreField(dm,"DMSwarm_rank",NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 170 ierr = DataBucketGetSizes(swarm->db,npoints_prior_migration,NULL,NULL);CHKERRQ(ierr); 171 172 ierr = DataExBegin(de);CHKERRQ(ierr); 173 ierr = DataExEnd(de);CHKERRQ(ierr); 174 175 ierr = DataExGetRecvData(de,&n_points_recv,(void**)&recv_points);CHKERRQ(ierr); 176 177 ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 178 ierr = DataBucketSetSizes(swarm->db,npoints + n_points_recv,-1);CHKERRQ(ierr); 179 for (p=0; p<n_points_recv; p++) { 180 void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point ); 181 182 ierr = DataBucketInsertPackedArray(swarm->db,npoints+p,data_p);CHKERRQ(ierr); 183 } 184 ierr = DataExView(de);CHKERRQ(ierr); 185 ierr = DataBucketDestroyPackedArray(swarm->db,&point_buffer);CHKERRQ(ierr); 186 ierr = DataExDestroy(de);CHKERRQ(ierr); 187 188 PetscFunctionReturn(0); 189 } 190 191 #undef __FUNCT__ 192 #define __FUNCT__ "DMSwarmMigrate_CellDM" 193 PetscErrorCode DMSwarmMigrate_CellDM(DM dm,PetscBool remove_sent_points) 194 { 195 DM_Swarm *swarm = (DM_Swarm*)dm->data; 196 PetscErrorCode ierr; 197 PetscInt p,npoints,npointsg,npoints2,npoints2g,len,*rankval,npoints_prior_migration; 198 const PetscInt *LA_iscell; 199 DM dmcell; 200 IS iscell; 201 Vec pos; 202 PetscBool error_check = swarm->migrate_error_on_missing_point; 203 204 ierr = DMSwarmGetCellDM(dm,&dmcell);CHKERRQ(ierr); 205 if (!dmcell) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only valid if cell DM provided"); 206 207 ierr = DMSwarmCreateGlobalVectorFromField(dm,"coor",&pos);CHKERRQ(ierr); 208 ierr = DMLocatePoints(dmcell,pos,&iscell);CHKERRQ(ierr); 209 ierr = DMSwarmDestroyGlobalVectorFromField(dm,"coor",&pos);CHKERRQ(ierr); 210 211 if (error_check) { 212 ierr = DMSwarmGetSize(dm,&npointsg);CHKERRQ(ierr); 213 } 214 ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 215 ierr = DMSwarmGetField(dm,"DMSwarm_rank",NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 216 ierr = ISGetIndices(iscell,&LA_iscell);CHKERRQ(ierr); 217 for (p=0; p<npoints; p++) { 218 if (LA_iscell[p] == -1) { 219 rankval[p] = -1; 220 } 221 } 222 ierr = ISRestoreIndices(iscell,&LA_iscell);CHKERRQ(ierr); 223 ierr = DMSwarmRestoreField(dm,"DMSwarm_rank",NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 224 ierr = ISDestroy(&iscell);CHKERRQ(ierr); 225 226 ierr = DMSwarmMigrate_DMNeighborScatter(dm,dmcell,remove_sent_points,&npoints_prior_migration);CHKERRQ(ierr); 227 228 /* locate points newly recevied */ 229 ierr = DataBucketGetSizes(swarm->db,&npoints2,NULL,NULL);CHKERRQ(ierr); 230 len = npoints2 - npoints_prior_migration; 231 if (len > 0) { 232 PetscScalar *LA_coor; 233 PetscInt bs = 1; 234 235 ierr = DMSwarmGetField(dm,"coor",NULL,NULL,(void**)&LA_coor);CHKERRQ(ierr); 236 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,bs,len,(const PetscScalar*)&LA_coor[npoints_prior_migration],&pos);CHKERRQ(ierr); 237 238 ierr = DMLocatePoints(dmcell,pos,&iscell);CHKERRQ(ierr); 239 240 ierr = VecDestroy(&pos);CHKERRQ(ierr); 241 ierr = DMSwarmRestoreField(dm,"coor",NULL,NULL,(void**)&LA_coor);CHKERRQ(ierr); 242 243 ierr = ISGetIndices(iscell,&LA_iscell);CHKERRQ(ierr); 244 ierr = DMSwarmGetField(dm,"DMSwarm_rank",NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 245 for (p=0; p<len; p++) { 246 if (LA_iscell[p] == -1) { 247 rankval[npoints_prior_migration+p] = -1; 248 } 249 } 250 ierr = ISRestoreIndices(iscell,&LA_iscell);CHKERRQ(ierr); 251 ierr = ISDestroy(&iscell);CHKERRQ(ierr); 252 253 /* remove points which left processor */ 254 ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 255 for (p=npoints_prior_migration; p<npoints2; p++) { 256 if (rankval[p] == -1) { 257 /* kill point */ 258 ierr = DataBucketRemovePointAtIndex(swarm->db,p);CHKERRQ(ierr); 259 DataBucketGetSizes(swarm->db,&npoints2,NULL,NULL);CHKERRQ(ierr); /* you need to update npoints as the list size decreases! */ 260 p--; /* check replacement point */ 261 } 262 } 263 ierr = DMSwarmRestoreField(dm,"DMSwarm_rank",NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 264 265 } 266 267 /* check for error on removed points */ 268 if (error_check) { 269 ierr = DMSwarmGetSize(dm,&npoints2g);CHKERRQ(ierr); 270 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); 271 } 272 PetscFunctionReturn(0); 273 } 274 275 /* 276 Redundant as this assumes points can only be sent to a single rank 277 */ 278 #undef __FUNCT__ 279 #define __FUNCT__ "DMSwarmMigrate_GlobalToLocal_Basic" 280 PetscErrorCode DMSwarmMigrate_GlobalToLocal_Basic(DM dm,PetscInt *globalsize) 281 { 282 DM_Swarm *swarm = (DM_Swarm*)dm->data; 283 PetscErrorCode ierr; 284 DataEx de; 285 PetscInt p,npoints,*rankval,n_points_recv; 286 PetscMPIInt rank,nrank,negrank; 287 void *point_buffer,*recv_points; 288 size_t sizeof_dmswarm_point; 289 290 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr); 291 292 ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 293 *globalsize = npoints; 294 ierr = DMSwarmGetField(dm,"DMSwarm_rank",NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 295 296 de = DataExCreate(PetscObjectComm((PetscObject)dm),0);CHKERRQ(ierr); 297 298 ierr = DataExTopologyInitialize(de);CHKERRQ(ierr); 299 for (p=0; p<npoints; p++) { 300 negrank = rankval[p]; 301 if (negrank < 0) { 302 nrank = -negrank - 1; 303 ierr = DataExTopologyAddNeighbour(de,nrank);CHKERRQ(ierr); 304 } 305 } 306 ierr = DataExTopologyFinalize(de);CHKERRQ(ierr); 307 308 ierr = DataExInitializeSendCount(de);CHKERRQ(ierr); 309 for (p=0; p<npoints; p++) { 310 negrank = rankval[p]; 311 if (negrank < 0) { 312 nrank = -negrank - 1; 313 ierr = DataExAddToSendCount(de,nrank,1);CHKERRQ(ierr); 314 } 315 } 316 ierr = DataExFinalizeSendCount(de);CHKERRQ(ierr); 317 318 ierr = DataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);CHKERRQ(ierr); 319 ierr = DataExPackInitialize(de,sizeof_dmswarm_point);CHKERRQ(ierr); 320 for (p=0; p<npoints; p++) { 321 negrank = rankval[p]; 322 if (negrank < 0) { 323 nrank = -negrank - 1; 324 rankval[p] = nrank; 325 /* copy point into buffer */ 326 ierr = DataBucketFillPackedArray(swarm->db,p,point_buffer);CHKERRQ(ierr); 327 /* insert point buffer into DataExchanger */ 328 ierr = DataExPackData(de,nrank,1,point_buffer);CHKERRQ(ierr); 329 rankval[p] = negrank; 330 } 331 } 332 ierr = DataExPackFinalize(de);CHKERRQ(ierr); 333 334 ierr = DMSwarmRestoreField(dm,"DMSwarm_rank",NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 335 336 ierr = DataExBegin(de);CHKERRQ(ierr); 337 ierr = DataExEnd(de);CHKERRQ(ierr); 338 339 ierr = DataExGetRecvData(de,&n_points_recv,(void**)&recv_points);CHKERRQ(ierr); 340 341 ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 342 ierr = DataBucketSetSizes(swarm->db,npoints + n_points_recv,-1);CHKERRQ(ierr); 343 for (p=0; p<n_points_recv; p++) { 344 void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point ); 345 346 ierr = DataBucketInsertPackedArray(swarm->db,npoints+p,data_p);CHKERRQ(ierr); 347 } 348 ierr = DataExView(de);CHKERRQ(ierr); 349 ierr = DataBucketDestroyPackedArray(swarm->db,&point_buffer);CHKERRQ(ierr); 350 ierr = DataExDestroy(de);CHKERRQ(ierr); 351 352 PetscFunctionReturn(0); 353 } 354 355 typedef struct { 356 PetscMPIInt owner_rank; 357 PetscReal min[3],max[3]; 358 } CollectBBox; 359 360 #undef __FUNCT__ 361 #define __FUNCT__ "DMSwarmCollect_DMDABoundingBox" 362 PETSC_EXTERN PetscErrorCode DMSwarmCollect_DMDABoundingBox(DM dm,PetscInt *globalsize) 363 { 364 DM_Swarm *swarm = (DM_Swarm*)dm->data; 365 PetscErrorCode ierr; 366 DataEx de; 367 PetscInt p,pk,npoints,*rankval,n_points_recv,n_bbox_recv,dim,neighbour_cells; 368 PetscMPIInt rank,nrank; 369 void *point_buffer,*recv_points; 370 size_t sizeof_dmswarm_point,sizeof_bbox_ctx; 371 PetscBool isdmda; 372 CollectBBox *bbox,*recv_bbox; 373 const PetscMPIInt *dmneighborranks; 374 DM dmcell; 375 376 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr); 377 378 ierr = DMSwarmGetCellDM(dm,&dmcell);CHKERRQ(ierr); 379 if (!dmcell) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only valid if cell DM provided"); 380 381 isdmda = PETSC_FALSE; 382 PetscObjectTypeCompare((PetscObject)dmcell,DMDA,&isdmda); 383 if (!isdmda) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only DMDA support for CollectBoundingBox"); 384 385 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr); 386 387 sizeof_bbox_ctx = sizeof(CollectBBox); 388 PetscMalloc1(1,&bbox); 389 bbox->owner_rank = rank; 390 391 /* compute the bounding box based on the overlapping / stenctil size */ 392 //ierr = DMDAGetLocalBoundingBox(dmcell,bbox->min,bbox->max);CHKERRQ(ierr); 393 { 394 Vec lcoor; 395 396 ierr = DMGetCoordinatesLocal(dmcell,&lcoor);CHKERRQ(ierr); 397 398 if (dim >= 1) { 399 ierr = VecStrideMin(lcoor,0,NULL,&bbox->min[0]);CHKERRQ(ierr); 400 ierr = VecStrideMax(lcoor,0,NULL,&bbox->max[0]);CHKERRQ(ierr); 401 } 402 403 if (dim >= 2) { 404 ierr = VecStrideMin(lcoor,1,NULL,&bbox->min[1]);CHKERRQ(ierr); 405 ierr = VecStrideMax(lcoor,1,NULL,&bbox->max[1]);CHKERRQ(ierr); 406 } 407 408 if (dim == 3) { 409 ierr = VecStrideMin(lcoor,2,NULL,&bbox->min[2]);CHKERRQ(ierr); 410 ierr = VecStrideMax(lcoor,2,NULL,&bbox->max[2]);CHKERRQ(ierr); 411 } 412 } 413 414 ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 415 *globalsize = npoints; 416 ierr = DMSwarmGetField(dm,"DMSwarm_rank",NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 417 418 de = DataExCreate(PetscObjectComm((PetscObject)dm),0);CHKERRQ(ierr); 419 420 /* use DMDA neighbours */ 421 ierr = DMDAGetNeighbors(dmcell,&dmneighborranks);CHKERRQ(ierr); 422 if (dim == 1) { 423 neighbour_cells = 3; 424 } else if (dim == 2) { 425 neighbour_cells = 9; 426 } else { 427 neighbour_cells = 27; 428 } 429 430 ierr = DataExTopologyInitialize(de);CHKERRQ(ierr); 431 for (p=0; p<neighbour_cells; p++) { 432 if ( (dmneighborranks[p] >= 0) && (dmneighborranks[p] != rank) ) { 433 ierr = DataExTopologyAddNeighbour(de,dmneighborranks[p]);CHKERRQ(ierr); 434 } 435 } 436 ierr = DataExTopologyFinalize(de);CHKERRQ(ierr); 437 438 ierr = DataExInitializeSendCount(de);CHKERRQ(ierr); 439 for (p=0; p<neighbour_cells; p++) { 440 if ( (dmneighborranks[p] >= 0) && (dmneighborranks[p] != rank) ) { 441 ierr = DataExAddToSendCount(de,dmneighborranks[p],1);CHKERRQ(ierr); 442 } 443 } 444 ierr = DataExFinalizeSendCount(de);CHKERRQ(ierr); 445 446 447 /* send bounding boxes */ 448 ierr = DataExPackInitialize(de,sizeof_bbox_ctx);CHKERRQ(ierr); 449 for (p=0; p<neighbour_cells; p++) { 450 nrank = dmneighborranks[p]; 451 if ( (nrank >= 0) && (nrank != rank) ) { 452 /* insert bbox buffer into DataExchanger */ 453 ierr = DataExPackData(de,nrank,1,bbox);CHKERRQ(ierr); 454 } 455 } 456 ierr = DataExPackFinalize(de);CHKERRQ(ierr); 457 458 /* recv bounding boxes */ 459 ierr = DataExBegin(de);CHKERRQ(ierr); 460 ierr = DataExEnd(de);CHKERRQ(ierr); 461 462 ierr = DataExGetRecvData(de,&n_bbox_recv,(void**)&recv_bbox);CHKERRQ(ierr); 463 464 465 for (p=0; p<n_bbox_recv; p++) { 466 printf("[rank %d]: box from %d : range[%+1.4e,%+1.4e]x[%+1.4e,%+1.4e]\n",rank,recv_bbox[p].owner_rank, 467 recv_bbox[p].min[0],recv_bbox[p].max[0],recv_bbox[p].min[1],recv_bbox[p].max[1]); 468 } 469 470 /* of course this is stupid as this "generic" function should have a better way to know what the coordinates are called */ 471 ierr = DataExInitializeSendCount(de);CHKERRQ(ierr); 472 for (pk=0; pk<n_bbox_recv; pk++) { 473 PetscReal *array_x,*array_y; 474 475 ierr = DMSwarmGetField(dm,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr); 476 ierr = DMSwarmGetField(dm,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr); 477 478 for (p=0; p<npoints; p++) { 479 if ((array_x[p] >= recv_bbox[pk].min[0]) && (array_x[p] <= recv_bbox[pk].max[0]) ) { 480 if ((array_y[p] >= recv_bbox[pk].min[1]) && (array_y[p] <= recv_bbox[pk].max[1]) ) { 481 482 ierr = DataExAddToSendCount(de,recv_bbox[pk].owner_rank,1);CHKERRQ(ierr); 483 484 } 485 } 486 } 487 ierr = DMSwarmRestoreField(dm,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr); 488 ierr = DMSwarmRestoreField(dm,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr); 489 } 490 ierr = DataExFinalizeSendCount(de);CHKERRQ(ierr); 491 492 493 494 ierr = DataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);CHKERRQ(ierr); 495 ierr = DataExPackInitialize(de,sizeof_dmswarm_point);CHKERRQ(ierr); 496 for (pk=0; pk<n_bbox_recv; pk++) { 497 PetscReal *array_x,*array_y; 498 499 ierr = DMSwarmGetField(dm,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr); 500 ierr = DMSwarmGetField(dm,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr); 501 502 for (p=0; p<npoints; p++) { 503 if ((array_x[p] >= recv_bbox[pk].min[0]) && (array_x[p] <= recv_bbox[pk].max[0]) ) { 504 if ((array_y[p] >= recv_bbox[pk].min[1]) && (array_y[p] <= recv_bbox[pk].max[1]) ) { 505 // copy point into buffer // 506 ierr = DataBucketFillPackedArray(swarm->db,p,point_buffer);CHKERRQ(ierr); 507 // insert point buffer into DataExchanger // 508 ierr = DataExPackData(de,recv_bbox[pk].owner_rank,1,point_buffer);CHKERRQ(ierr); 509 } 510 } 511 } 512 513 ierr = DMSwarmRestoreField(dm,"coory",NULL,NULL,(void**)&array_y);CHKERRQ(ierr); 514 ierr = DMSwarmRestoreField(dm,"coorx",NULL,NULL,(void**)&array_x);CHKERRQ(ierr); 515 } 516 517 ierr = DataExPackFinalize(de);CHKERRQ(ierr); 518 519 ierr = DMSwarmRestoreField(dm,"DMSwarm_rank",NULL,NULL,(void**)&rankval);CHKERRQ(ierr); 520 521 ierr = DataExBegin(de);CHKERRQ(ierr); 522 ierr = DataExEnd(de);CHKERRQ(ierr); 523 524 ierr = DataExGetRecvData(de,&n_points_recv,(void**)&recv_points);CHKERRQ(ierr); 525 526 ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 527 ierr = DataBucketSetSizes(swarm->db,npoints + n_points_recv,-1);CHKERRQ(ierr); 528 for (p=0; p<n_points_recv; p++) { 529 void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point ); 530 531 ierr = DataBucketInsertPackedArray(swarm->db,npoints+p,data_p);CHKERRQ(ierr); 532 } 533 534 ierr = DataBucketDestroyPackedArray(swarm->db,&point_buffer);CHKERRQ(ierr); 535 PetscFree(bbox); 536 ierr = DataExView(de);CHKERRQ(ierr); 537 ierr = DataExDestroy(de);CHKERRQ(ierr); 538 539 PetscFunctionReturn(0); 540 } 541 542 543 /* General collection when no order, or neighbour information is provided */ 544 /* 545 User provides context and collect() method 546 Broadcast user context 547 548 for each context / rank { 549 collect(swarm,context,n,list) 550 } 551 552 */ 553 #undef __FUNCT__ 554 #define __FUNCT__ "DMSwarmCollect_General" 555 PETSC_EXTERN PetscErrorCode DMSwarmCollect_General(DM dm,PetscErrorCode (*collect)(DM,void*,PetscInt*,PetscInt**),size_t ctx_size,void *ctx,PetscInt *globalsize) 556 { 557 DM_Swarm *swarm = (DM_Swarm*)dm->data; 558 PetscErrorCode ierr; 559 DataEx de; 560 PetscInt p,r,npoints,n_points_recv; 561 PetscMPIInt commsize,rank; 562 void *point_buffer,*recv_points; 563 void *ctxlist; 564 PetscInt *n2collect,**collectlist; 565 size_t sizeof_dmswarm_point; 566 567 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm),&commsize);CHKERRQ(ierr); 568 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr); 569 570 ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 571 *globalsize = npoints; 572 573 /* Broadcast user context */ 574 PetscMalloc(ctx_size*commsize,&ctxlist); 575 ierr = MPI_Allgather(ctx,ctx_size,MPI_CHAR,ctxlist,ctx_size,MPI_CHAR,PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 576 577 PetscMalloc1(commsize,&n2collect); 578 PetscMalloc1(commsize,&collectlist); 579 580 for (r=0; r<commsize; r++) { 581 PetscInt _n2collect; 582 PetscInt *_collectlist; 583 void *_ctx_r; 584 585 _n2collect = 0; 586 _collectlist = NULL; 587 588 if (r != rank) { /* don't collect data from yourself */ 589 _ctx_r = (void*)( (char*)ctxlist + r * ctx_size ); 590 ierr = collect(dm,_ctx_r,&_n2collect,&_collectlist);CHKERRQ(ierr); 591 } 592 593 n2collect[r] = _n2collect; 594 collectlist[r] = _collectlist; 595 } 596 597 de = DataExCreate(PetscObjectComm((PetscObject)dm),0);CHKERRQ(ierr); 598 599 /* Define topology */ 600 ierr = DataExTopologyInitialize(de);CHKERRQ(ierr); 601 for (r=0; r<commsize; r++) { 602 if (n2collect[r] > 0) { 603 ierr = DataExTopologyAddNeighbour(de,(PetscMPIInt)r);CHKERRQ(ierr); 604 } 605 } 606 ierr = DataExTopologyFinalize(de);CHKERRQ(ierr); 607 608 /* Define send counts */ 609 ierr = DataExInitializeSendCount(de);CHKERRQ(ierr); 610 for (r=0; r<commsize; r++) { 611 if (n2collect[r] > 0) { 612 ierr = DataExAddToSendCount(de,r,n2collect[r]);CHKERRQ(ierr); 613 } 614 } 615 ierr = DataExFinalizeSendCount(de);CHKERRQ(ierr); 616 617 /* Pack data */ 618 ierr = DataBucketCreatePackedArray(swarm->db,&sizeof_dmswarm_point,&point_buffer);CHKERRQ(ierr); 619 620 ierr = DataExPackInitialize(de,sizeof_dmswarm_point);CHKERRQ(ierr); 621 for (r=0; r<commsize; r++) { 622 623 for (p=0; p<n2collect[r]; p++) { 624 ierr = DataBucketFillPackedArray(swarm->db,collectlist[r][p],point_buffer);CHKERRQ(ierr); 625 /* insert point buffer into the data exchanger */ 626 ierr = DataExPackData(de,r,1,point_buffer);CHKERRQ(ierr); 627 } 628 } 629 630 ierr = DataExPackFinalize(de);CHKERRQ(ierr); 631 632 /* Scatter */ 633 ierr = DataExBegin(de);CHKERRQ(ierr); 634 ierr = DataExEnd(de);CHKERRQ(ierr); 635 636 /* Collect data in DMSwarm container */ 637 ierr = DataExGetRecvData(de,&n_points_recv,(void**)&recv_points);CHKERRQ(ierr); 638 639 ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr); 640 ierr = DataBucketSetSizes(swarm->db,npoints + n_points_recv,-1);CHKERRQ(ierr); 641 for (p=0; p<n_points_recv; p++) { 642 void *data_p = (void*)( (char*)recv_points + p*sizeof_dmswarm_point ); 643 644 ierr = DataBucketInsertPackedArray(swarm->db,npoints+p,data_p);CHKERRQ(ierr); 645 } 646 647 /* Release memory */ 648 for (r=0; r<commsize; r++) { 649 if (collectlist[r]) PetscFree(collectlist[r]); 650 } 651 PetscFree(collectlist); 652 PetscFree(n2collect); 653 PetscFree(ctxlist); 654 ierr = DataBucketDestroyPackedArray(swarm->db,&point_buffer);CHKERRQ(ierr); 655 ierr = DataExView(de);CHKERRQ(ierr); 656 ierr = DataExDestroy(de);CHKERRQ(ierr); 657 658 PetscFunctionReturn(0); 659 } 660 661