1 2 #define PETSCDM_DLL 3 #include <petsc/private/dmswarmimpl.h> /*I "petscdmswarm.h" I*/ 4 #include <petscsf.h> 5 6 /* 7 Error chceking macto to ensure the swarm type is correct and that a cell DM has been set 8 */ 9 #define DMSWARMPICVALID(dm) \ 10 { \ 11 DM_Swarm *_swarm = (DM_Swarm*)(dm)->data; \ 12 if (_swarm->swarm_type != DMSWARM_PIC) SETERRQ(PetscObjectComm((PetscObject)(dm)),PETSC_ERR_SUP,"Only valid for DMSwarm-PIC. You must call DMSwarmSetType(dm,DMSWARM_PIC)"); \ 13 else \ 14 if (!_swarm->dmcell) SETERRQ(PetscObjectComm((PetscObject)(dm)),PETSC_ERR_SUP,"Only valid for DMSwarmPIC if the cell DM is set. You must call DMSwarmSetCellDM(dm,celldm)"); \ 15 } 16 17 /* Coordinate insertition/addition API */ 18 /*@C 19 DMSwarmSetPointsUniformCoordinates - Set point coordinates in a DMSwarm on a regular (ijk) grid 20 21 Collective on DM 22 23 Input parameters: 24 + dm - the DMSwarm 25 . min - minimum coordinate values in the x, y, z directions (array of length dim) 26 . max - maximum coordinate values in the x, y, z directions (array of length dim) 27 . npoints - number of points in each spatial direction (array of length dim) 28 - mode - indicates whether to append points to the swarm (ADD_VALUES), or over-ride existing points (INSERT_VALUES) 29 30 Level: beginner 31 32 Notes: 33 When using mode = INSERT_VALUES, this method will reset the number of particles in the DMSwarm 34 to be npoints[0]*npoints[1] (2D) or npoints[0]*npoints[1]*npoints[2] (3D). When using mode = ADD_VALUES, 35 new points will be appended to any already existing in the DMSwarm 36 37 .seealso: DMSwarmSetType(), DMSwarmSetCellDM(), DMSwarmType 38 @*/ 39 PETSC_EXTERN PetscErrorCode DMSwarmSetPointsUniformCoordinates(DM dm,PetscReal min[],PetscReal max[],PetscInt npoints[],InsertMode mode) 40 { 41 PetscErrorCode ierr; 42 PetscReal gmin[] = {PETSC_MAX_REAL ,PETSC_MAX_REAL, PETSC_MAX_REAL}; 43 PetscReal gmax[] = {PETSC_MIN_REAL, PETSC_MIN_REAL, PETSC_MIN_REAL}; 44 PetscInt i,j,k,N,bs,b,n_estimate,n_curr,n_new_est,p,n_found; 45 Vec coorlocal; 46 const PetscScalar *_coor; 47 DM celldm; 48 PetscReal dx[3]; 49 Vec pos; 50 PetscScalar *_pos; 51 PetscReal *swarm_coor; 52 PetscInt *swarm_cellid; 53 PetscSF sfcell = NULL; 54 const PetscSFNode *LA_sfcell; 55 56 PetscFunctionBegin; 57 DMSWARMPICVALID(dm); 58 ierr = DMSwarmGetCellDM(dm,&celldm);CHKERRQ(ierr); 59 ierr = DMGetCoordinatesLocal(celldm,&coorlocal);CHKERRQ(ierr); 60 ierr = VecGetSize(coorlocal,&N);CHKERRQ(ierr); 61 ierr = VecGetBlockSize(coorlocal,&bs);CHKERRQ(ierr); 62 N = N / bs; 63 ierr = VecGetArrayRead(coorlocal,&_coor);CHKERRQ(ierr); 64 for (i=0; i<N; i++) { 65 for (b=0; b<bs; b++) { 66 gmin[b] = PetscMin(gmin[b],_coor[bs*i+b]); 67 gmax[b] = PetscMax(gmax[b],_coor[bs*i+b]); 68 } 69 } 70 ierr = VecRestoreArrayRead(coorlocal,&_coor);CHKERRQ(ierr); 71 72 for (b=0; b<bs; b++) { 73 dx[b] = (max[b] - min[b])/((PetscReal)(npoints[b]-1)); 74 } 75 76 /* determine number of points living in the bounding box */ 77 n_estimate = 0; 78 if (bs == 2) { npoints[2] = 1; } 79 for (k=0; k<npoints[2]; k++) { 80 for (j=0; j<npoints[1]; j++) { 81 for (i=0; i<npoints[0]; i++) { 82 PetscReal xp[] = {0.0,0.0,0.0}; 83 PetscInt ijk[3]; 84 PetscBool point_inside = PETSC_TRUE; 85 86 ijk[0] = i; 87 ijk[1] = j; 88 ijk[2] = k; 89 for (b=0; b<bs; b++) { 90 xp[b] = min[b] + ijk[b] * dx[b]; 91 } 92 for (b=0; b<bs; b++) { 93 if (xp[b] < gmin[b]) { point_inside = PETSC_FALSE; } 94 if (xp[b] > gmax[b]) { point_inside = PETSC_FALSE; } 95 } 96 if (point_inside) { n_estimate++; } 97 } 98 } 99 } 100 101 /* create candidate list */ 102 ierr = VecCreate(PetscObjectComm((PetscObject)dm),&pos);CHKERRQ(ierr); 103 ierr = VecSetSizes(pos,bs*n_estimate,PETSC_DECIDE);CHKERRQ(ierr); 104 ierr = VecSetBlockSize(pos,bs);CHKERRQ(ierr); 105 ierr = VecSetFromOptions(pos);CHKERRQ(ierr); 106 ierr = VecGetArray(pos,&_pos);CHKERRQ(ierr); 107 108 n_estimate = 0; 109 for (k=0; k<npoints[2]; k++) { 110 for (j=0; j<npoints[1]; j++) { 111 for (i=0; i<npoints[0]; i++) { 112 PetscReal xp[] = {0.0,0.0,0.0}; 113 PetscInt ijk[3]; 114 PetscBool point_inside = PETSC_TRUE; 115 116 ijk[0] = i; 117 ijk[1] = j; 118 ijk[2] = k; 119 for (b=0; b<bs; b++) { 120 xp[b] = min[b] + ijk[b] * dx[b]; 121 } 122 for (b=0; b<bs; b++) { 123 if (xp[b] < gmin[b]) { point_inside = PETSC_FALSE; } 124 if (xp[b] > gmax[b]) { point_inside = PETSC_FALSE; } 125 } 126 if (point_inside) { 127 for (b=0; b<bs; b++) { 128 _pos[bs*n_estimate+b] = xp[b]; 129 } 130 n_estimate++; 131 } 132 } 133 } 134 } 135 ierr = VecRestoreArray(pos,&_pos);CHKERRQ(ierr); 136 137 /* locate points */ 138 ierr = DMLocatePoints(celldm,pos,DM_POINTLOCATION_NONE,&sfcell);CHKERRQ(ierr); 139 140 ierr = PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell);CHKERRQ(ierr); 141 n_found = 0; 142 for (p=0; p<n_estimate; p++) { 143 if (LA_sfcell[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) { 144 n_found++; 145 } 146 } 147 148 /* adjust size */ 149 if (mode == ADD_VALUES) { 150 ierr = DMSwarmGetLocalSize(dm,&n_curr);CHKERRQ(ierr); 151 n_new_est = n_curr + n_found; 152 ierr = DMSwarmSetLocalSizes(dm,n_new_est,-1);CHKERRQ(ierr); 153 } 154 if (mode == INSERT_VALUES) { 155 n_curr = 0; 156 n_new_est = n_found; 157 ierr = DMSwarmSetLocalSizes(dm,n_new_est,-1);CHKERRQ(ierr); 158 } 159 160 /* initialize new coords, cell owners, pid */ 161 ierr = VecGetArrayRead(pos,&_coor);CHKERRQ(ierr); 162 ierr = DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr); 163 ierr = DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr); 164 n_found = 0; 165 for (p=0; p<n_estimate; p++) { 166 if (LA_sfcell[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) { 167 for (b=0; b<bs; b++) { 168 swarm_coor[bs*(n_curr + n_found) + b] = _coor[bs*p+b]; 169 } 170 swarm_cellid[n_curr + n_found] = LA_sfcell[p].index; 171 n_found++; 172 } 173 } 174 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr); 175 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr); 176 ierr = VecRestoreArrayRead(pos,&_coor);CHKERRQ(ierr); 177 178 ierr = PetscSFDestroy(&sfcell);CHKERRQ(ierr); 179 ierr = VecDestroy(&pos);CHKERRQ(ierr); 180 181 PetscFunctionReturn(0); 182 } 183 184 /*@C 185 DMSwarmSetPointCoordinates - Set point coordinates in a DMSwarm from a user defined list 186 187 Collective on DM 188 189 Input parameters: 190 + dm - the DMSwarm 191 . npoints - the number of points to insert 192 . coor - the coordinate values 193 . redundant - if set to PETSC_TRUE, it is assumed that npoints and coor[] are only valid on rank 0 and should be broadcast to other ranks 194 - mode - indicates whether to append points to the swarm (ADD_VALUES), or over-ride existing points (INSERT_VALUES) 195 196 Level: beginner 197 198 Notes: 199 If the user has specified redundant = PETSC_FALSE, the cell DM will attempt to locate the coordinates provided by coor[] within 200 its sub-domain. If they any values within coor[] are not located in the sub-domain, they will be ignored and will not get 201 added to the DMSwarm. 202 203 204 .seealso: DMSwarmSetType(), DMSwarmSetCellDM(), DMSwarmType, DMSwarmSetPointsUniformCoordinates() 205 @*/ 206 PETSC_EXTERN PetscErrorCode DMSwarmSetPointCoordinates(DM dm,PetscInt npoints,PetscReal coor[],PetscBool redundant,InsertMode mode) 207 { 208 PetscErrorCode ierr; 209 PetscReal gmin[] = {PETSC_MAX_REAL ,PETSC_MAX_REAL, PETSC_MAX_REAL}; 210 PetscReal gmax[] = {PETSC_MIN_REAL, PETSC_MIN_REAL, PETSC_MIN_REAL}; 211 PetscInt i,N,bs,b,n_estimate,n_curr,n_new_est,p,n_found; 212 Vec coorlocal; 213 const PetscScalar *_coor; 214 DM celldm; 215 Vec pos; 216 PetscScalar *_pos; 217 PetscReal *swarm_coor; 218 PetscInt *swarm_cellid; 219 PetscSF sfcell = NULL; 220 const PetscSFNode *LA_sfcell; 221 PetscReal *my_coor; 222 PetscInt my_npoints; 223 PetscMPIInt rank; 224 MPI_Comm comm; 225 226 PetscFunctionBegin; 227 DMSWARMPICVALID(dm); 228 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 229 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 230 231 ierr = DMSwarmGetCellDM(dm,&celldm);CHKERRQ(ierr); 232 ierr = DMGetCoordinatesLocal(celldm,&coorlocal);CHKERRQ(ierr); 233 ierr = VecGetSize(coorlocal,&N);CHKERRQ(ierr); 234 ierr = VecGetBlockSize(coorlocal,&bs);CHKERRQ(ierr); 235 N = N / bs; 236 ierr = VecGetArrayRead(coorlocal,&_coor);CHKERRQ(ierr); 237 for (i=0; i<N; i++) { 238 for (b=0; b<bs; b++) { 239 gmin[b] = PetscMin(gmin[b],_coor[bs*i+b]); 240 gmax[b] = PetscMax(gmax[b],_coor[bs*i+b]); 241 } 242 } 243 ierr = VecRestoreArrayRead(coorlocal,&_coor);CHKERRQ(ierr); 244 245 /* broadcast points from rank 0 if requested */ 246 if (redundant) { 247 my_npoints = npoints; 248 ierr = MPI_Bcast(&my_npoints,1,MPIU_INT,0,comm);CHKERRQ(ierr); 249 250 if (rank > 0) { /* allocate space */ 251 ierr = PetscMalloc1(my_npoints,&my_coor);CHKERRQ(ierr); 252 } else { 253 my_coor = coor; 254 } 255 ierr = MPI_Bcast(my_coor,bs*my_npoints,MPIU_REAL,0,comm);CHKERRQ(ierr); 256 } else { 257 my_npoints = npoints; 258 my_coor = coor; 259 } 260 261 /* determine the number of points living in the bounding box */ 262 n_estimate = 0; 263 for (i=0; i<my_npoints; i++) { 264 PetscBool point_inside = PETSC_TRUE; 265 266 for (b=0; b<bs; b++) { 267 if (my_coor[bs*i+b] < gmin[b]) { point_inside = PETSC_FALSE; } 268 if (my_coor[bs*i+b] > gmax[b]) { point_inside = PETSC_FALSE; } 269 } 270 if (point_inside) { n_estimate++; } 271 } 272 273 /* create candidate list */ 274 ierr = VecCreate(PetscObjectComm((PetscObject)dm),&pos);CHKERRQ(ierr); 275 ierr = VecSetSizes(pos,bs*n_estimate,PETSC_DECIDE);CHKERRQ(ierr); 276 ierr = VecSetBlockSize(pos,bs);CHKERRQ(ierr); 277 ierr = VecSetFromOptions(pos);CHKERRQ(ierr); 278 ierr = VecGetArray(pos,&_pos);CHKERRQ(ierr); 279 280 n_estimate = 0; 281 for (i=0; i<my_npoints; i++) { 282 PetscBool point_inside = PETSC_TRUE; 283 284 for (b=0; b<bs; b++) { 285 if (my_coor[bs*i+b] < gmin[b]) { point_inside = PETSC_FALSE; } 286 if (my_coor[bs*i+b] > gmax[b]) { point_inside = PETSC_FALSE; } 287 } 288 if (point_inside) { 289 for (b=0; b<bs; b++) { 290 _pos[bs*n_estimate+b] = my_coor[bs*i+b]; 291 } 292 n_estimate++; 293 } 294 } 295 ierr = VecRestoreArray(pos,&_pos);CHKERRQ(ierr); 296 297 /* locate points */ 298 ierr = DMLocatePoints(celldm,pos,DM_POINTLOCATION_NONE,&sfcell);CHKERRQ(ierr); 299 300 ierr = PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell);CHKERRQ(ierr); 301 n_found = 0; 302 for (p=0; p<n_estimate; p++) { 303 if (LA_sfcell[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) { 304 n_found++; 305 } 306 } 307 308 /* adjust size */ 309 if (mode == ADD_VALUES) { 310 ierr = DMSwarmGetLocalSize(dm,&n_curr);CHKERRQ(ierr); 311 n_new_est = n_curr + n_found; 312 ierr = DMSwarmSetLocalSizes(dm,n_new_est,-1);CHKERRQ(ierr); 313 } 314 if (mode == INSERT_VALUES) { 315 n_curr = 0; 316 n_new_est = n_found; 317 ierr = DMSwarmSetLocalSizes(dm,n_new_est,-1);CHKERRQ(ierr); 318 } 319 320 /* initialize new coords, cell owners, pid */ 321 ierr = VecGetArrayRead(pos,&_coor);CHKERRQ(ierr); 322 ierr = DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr); 323 ierr = DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr); 324 n_found = 0; 325 for (p=0; p<n_estimate; p++) { 326 if (LA_sfcell[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) { 327 for (b=0; b<bs; b++) { 328 swarm_coor[bs*(n_curr + n_found) + b] = _coor[bs*p+b]; 329 } 330 swarm_cellid[n_curr + n_found] = LA_sfcell[p].index; 331 n_found++; 332 } 333 } 334 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr); 335 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr); 336 ierr = VecRestoreArrayRead(pos,&_coor);CHKERRQ(ierr); 337 338 if (redundant) { 339 if (rank > 0) { 340 ierr = PetscFree(my_coor);CHKERRQ(ierr); 341 } 342 } 343 ierr = PetscSFDestroy(&sfcell);CHKERRQ(ierr); 344 ierr = VecDestroy(&pos);CHKERRQ(ierr); 345 346 PetscFunctionReturn(0); 347 } 348 349 extern PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_DA(DM,DM,DMSwarmPICLayoutType,PetscInt); 350 extern PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX(DM,DM,DMSwarmPICLayoutType,PetscInt); 351 352 /*@C 353 DMSwarmInsertPointsUsingCellDM - Insert point coordinates within each cell 354 355 Not collective 356 357 Input parameters: 358 + dm - the DMSwarm 359 . layout_type - method used to fill each cell with the cell DM 360 - fill_param - parameter controlling how many points per cell are added (the meaning of this parameter is dependent on the layout type) 361 362 Level: beginner 363 364 Notes: 365 The insert method will reset any previous defined points within the DMSwarm 366 367 .seealso: DMSwarmPICLayoutType, DMSwarmSetType(), DMSwarmSetCellDM(), DMSwarmType 368 @*/ 369 PETSC_EXTERN PetscErrorCode DMSwarmInsertPointsUsingCellDM(DM dm,DMSwarmPICLayoutType layout_type,PetscInt fill_param) 370 { 371 PetscErrorCode ierr; 372 DM celldm; 373 PetscBool isDA,isPLEX; 374 375 PetscFunctionBegin; 376 DMSWARMPICVALID(dm); 377 ierr = DMSwarmGetCellDM(dm,&celldm);CHKERRQ(ierr); 378 ierr = PetscObjectTypeCompare((PetscObject)celldm,DMDA,&isDA);CHKERRQ(ierr); 379 ierr = PetscObjectTypeCompare((PetscObject)celldm,DMPLEX,&isPLEX);CHKERRQ(ierr); 380 if (isDA) { 381 ierr = private_DMSwarmInsertPointsUsingCellDM_DA(dm,celldm,layout_type,fill_param);CHKERRQ(ierr); 382 } else if (isPLEX) { 383 ierr = private_DMSwarmInsertPointsUsingCellDM_PLEX(dm,celldm,layout_type,fill_param);CHKERRQ(ierr); 384 } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only supported for cell DMs of type DMDA and DMPLEX"); 385 386 PetscFunctionReturn(0); 387 } 388 389 /* 390 PETSC_EXTERN PetscErrorCode DMSwarmAddPointCoordinatesCellWise(DM dm,PetscInt cell,PetscInt npoints,PetscReal xi[],PetscBool proximity_initialization) 391 { 392 PetscFunctionBegin; 393 PetscFunctionReturn(0); 394 } 395 */ 396 397 /* Field projection API */ 398 /* 399 PETSC_EXTERN PetscErrorCode DMSwarmProjectFields(DM dm,PetscInt project_type,PetscInt nfields,const char *fieldnames[],Vec *fields) 400 { 401 PetscFunctionBegin; 402 PetscFunctionReturn(0); 403 } 404 */ 405 406 /*@C 407 DMSwarmCreatePointPerCellCount - Count the number of points (particles) per cell 408 409 Not collective 410 411 Input parameter: 412 . dm - the DMSwarm 413 414 Output parameters: 415 + ncells - the number of cells in the cell DM (optional argument, pass NULL to ignore) 416 - count - array of length ncells containing the number of particles per cell 417 418 Level: beginner 419 420 Notes: 421 The array count is allocated internally and must be free'd by the user. 422 423 .seealso: DMSwarmSetType(), DMSwarmSetCellDM(), DMSwarmType 424 @*/ 425 /* 426 PETSC_EXTERN PetscErrorCode DMSwarmCreatePointPerCellCount(DM dm,PetscInt *ncells,PetscInt **count) 427 { 428 PetscFunctionBegin; 429 PetscFunctionReturn(0); 430 } 431 */ 432