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