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 .seealso: DMSwarmSetType(), DMSwarmSetCellDM(), DMSwarmType, DMSwarmSetPointsUniformCoordinates() 204 @*/ 205 PETSC_EXTERN PetscErrorCode DMSwarmSetPointCoordinates(DM dm,PetscInt npoints,PetscReal coor[],PetscBool redundant,InsertMode mode) 206 { 207 PetscErrorCode ierr; 208 PetscReal gmin[] = {PETSC_MAX_REAL ,PETSC_MAX_REAL, PETSC_MAX_REAL}; 209 PetscReal gmax[] = {PETSC_MIN_REAL, PETSC_MIN_REAL, PETSC_MIN_REAL}; 210 PetscInt i,N,bs,b,n_estimate,n_curr,n_new_est,p,n_found; 211 Vec coorlocal; 212 const PetscScalar *_coor; 213 DM celldm; 214 Vec pos; 215 PetscScalar *_pos; 216 PetscReal *swarm_coor; 217 PetscInt *swarm_cellid; 218 PetscSF sfcell = NULL; 219 const PetscSFNode *LA_sfcell; 220 PetscReal *my_coor; 221 PetscInt my_npoints; 222 PetscMPIInt rank; 223 MPI_Comm comm; 224 225 PetscFunctionBegin; 226 DMSWARMPICVALID(dm); 227 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 228 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 229 230 ierr = DMSwarmGetCellDM(dm,&celldm);CHKERRQ(ierr); 231 ierr = DMGetCoordinatesLocal(celldm,&coorlocal);CHKERRQ(ierr); 232 ierr = VecGetSize(coorlocal,&N);CHKERRQ(ierr); 233 ierr = VecGetBlockSize(coorlocal,&bs);CHKERRQ(ierr); 234 N = N / bs; 235 ierr = VecGetArrayRead(coorlocal,&_coor);CHKERRQ(ierr); 236 for (i=0; i<N; i++) { 237 for (b=0; b<bs; b++) { 238 gmin[b] = PetscMin(gmin[b],_coor[bs*i+b]); 239 gmax[b] = PetscMax(gmax[b],_coor[bs*i+b]); 240 } 241 } 242 ierr = VecRestoreArrayRead(coorlocal,&_coor);CHKERRQ(ierr); 243 244 /* broadcast points from rank 0 if requested */ 245 if (redundant) { 246 my_npoints = npoints; 247 ierr = MPI_Bcast(&my_npoints,1,MPIU_INT,0,comm);CHKERRQ(ierr); 248 249 if (rank > 0) { /* allocate space */ 250 ierr = PetscMalloc1(my_npoints,&my_coor);CHKERRQ(ierr); 251 } else { 252 my_coor = coor; 253 } 254 ierr = MPI_Bcast(my_coor,bs*my_npoints,MPIU_REAL,0,comm);CHKERRQ(ierr); 255 } else { 256 my_npoints = npoints; 257 my_coor = coor; 258 } 259 260 /* determine the number of points living in the bounding box */ 261 n_estimate = 0; 262 for (i=0; i<my_npoints; i++) { 263 PetscBool point_inside = PETSC_TRUE; 264 265 for (b=0; b<bs; b++) { 266 if (my_coor[bs*i+b] < gmin[b]) { point_inside = PETSC_FALSE; } 267 if (my_coor[bs*i+b] > gmax[b]) { point_inside = PETSC_FALSE; } 268 } 269 if (point_inside) { n_estimate++; } 270 } 271 272 /* create candidate list */ 273 ierr = VecCreate(PetscObjectComm((PetscObject)dm),&pos);CHKERRQ(ierr); 274 ierr = VecSetSizes(pos,bs*n_estimate,PETSC_DECIDE);CHKERRQ(ierr); 275 ierr = VecSetBlockSize(pos,bs);CHKERRQ(ierr); 276 ierr = VecSetFromOptions(pos);CHKERRQ(ierr); 277 ierr = VecGetArray(pos,&_pos);CHKERRQ(ierr); 278 279 n_estimate = 0; 280 for (i=0; i<my_npoints; i++) { 281 PetscBool point_inside = PETSC_TRUE; 282 283 for (b=0; b<bs; b++) { 284 if (my_coor[bs*i+b] < gmin[b]) { point_inside = PETSC_FALSE; } 285 if (my_coor[bs*i+b] > gmax[b]) { point_inside = PETSC_FALSE; } 286 } 287 if (point_inside) { 288 for (b=0; b<bs; b++) { 289 _pos[bs*n_estimate+b] = my_coor[bs*i+b]; 290 } 291 n_estimate++; 292 } 293 } 294 ierr = VecRestoreArray(pos,&_pos);CHKERRQ(ierr); 295 296 /* locate points */ 297 ierr = DMLocatePoints(celldm,pos,DM_POINTLOCATION_NONE,&sfcell);CHKERRQ(ierr); 298 299 ierr = PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell);CHKERRQ(ierr); 300 n_found = 0; 301 for (p=0; p<n_estimate; p++) { 302 if (LA_sfcell[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) { 303 n_found++; 304 } 305 } 306 307 /* adjust size */ 308 if (mode == ADD_VALUES) { 309 ierr = DMSwarmGetLocalSize(dm,&n_curr);CHKERRQ(ierr); 310 n_new_est = n_curr + n_found; 311 ierr = DMSwarmSetLocalSizes(dm,n_new_est,-1);CHKERRQ(ierr); 312 } 313 if (mode == INSERT_VALUES) { 314 n_curr = 0; 315 n_new_est = n_found; 316 ierr = DMSwarmSetLocalSizes(dm,n_new_est,-1);CHKERRQ(ierr); 317 } 318 319 /* initialize new coords, cell owners, pid */ 320 ierr = VecGetArrayRead(pos,&_coor);CHKERRQ(ierr); 321 ierr = DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr); 322 ierr = DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr); 323 n_found = 0; 324 for (p=0; p<n_estimate; p++) { 325 if (LA_sfcell[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) { 326 for (b=0; b<bs; b++) { 327 swarm_coor[bs*(n_curr + n_found) + b] = _coor[bs*p+b]; 328 } 329 swarm_cellid[n_curr + n_found] = LA_sfcell[p].index; 330 n_found++; 331 } 332 } 333 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr); 334 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr); 335 ierr = VecRestoreArrayRead(pos,&_coor);CHKERRQ(ierr); 336 337 if (redundant) { 338 if (rank > 0) { 339 ierr = PetscFree(my_coor);CHKERRQ(ierr); 340 } 341 } 342 ierr = PetscSFDestroy(&sfcell);CHKERRQ(ierr); 343 ierr = VecDestroy(&pos);CHKERRQ(ierr); 344 345 PetscFunctionReturn(0); 346 } 347 348 extern PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_DA(DM,DM,DMSwarmPICLayoutType,PetscInt); 349 extern PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX(DM,DM,DMSwarmPICLayoutType,PetscInt); 350 351 /*@C 352 DMSwarmInsertPointsUsingCellDM - Insert point coordinates within each cell 353 354 Not collective 355 356 Input parameters: 357 + dm - the DMSwarm 358 . layout_type - method used to fill each cell with the cell DM 359 - fill_param - parameter controlling how many points per cell are added (the meaning of this parameter is dependent on the layout type) 360 361 Level: beginner 362 363 Notes: 364 The insert method will reset any previous defined points within the DMSwarm 365 366 .seealso: DMSwarmPICLayoutType, DMSwarmSetType(), DMSwarmSetCellDM(), DMSwarmType 367 @*/ 368 PETSC_EXTERN PetscErrorCode DMSwarmInsertPointsUsingCellDM(DM dm,DMSwarmPICLayoutType layout_type,PetscInt fill_param) 369 { 370 PetscErrorCode ierr; 371 DM celldm; 372 PetscBool isDA,isPLEX; 373 374 PetscFunctionBegin; 375 DMSWARMPICVALID(dm); 376 ierr = DMSwarmGetCellDM(dm,&celldm);CHKERRQ(ierr); 377 ierr = PetscObjectTypeCompare((PetscObject)celldm,DMDA,&isDA);CHKERRQ(ierr); 378 ierr = PetscObjectTypeCompare((PetscObject)celldm,DMPLEX,&isPLEX);CHKERRQ(ierr); 379 if (isDA) { 380 ierr = private_DMSwarmInsertPointsUsingCellDM_DA(dm,celldm,layout_type,fill_param);CHKERRQ(ierr); 381 } else if (isPLEX) { 382 ierr = private_DMSwarmInsertPointsUsingCellDM_PLEX(dm,celldm,layout_type,fill_param);CHKERRQ(ierr); 383 } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only supported for cell DMs of type DMDA and DMPLEX"); 384 385 PetscFunctionReturn(0); 386 } 387 388 /* 389 PETSC_EXTERN PetscErrorCode DMSwarmAddPointCoordinatesCellWise(DM dm,PetscInt cell,PetscInt npoints,PetscReal xi[],PetscBool proximity_initialization) 390 { 391 PetscFunctionBegin; 392 PetscFunctionReturn(0); 393 } 394 */ 395 396 /* Field projection API */ 397 /* 398 PETSC_EXTERN PetscErrorCode DMSwarmProjectFields(DM dm,PetscInt project_type,PetscInt nfields,const char *fieldnames[],Vec *fields) 399 { 400 PetscFunctionBegin; 401 PetscFunctionReturn(0); 402 } 403 */ 404 405 /*@C 406 DMSwarmCreatePointPerCellCount - Count the number of points (particles) per cell 407 408 Not collective 409 410 Input parameter: 411 . dm - the DMSwarm 412 413 Output parameters: 414 + ncells - the number of cells in the cell DM (optional argument, pass NULL to ignore) 415 - count - array of length ncells containing the number of particles per cell 416 417 Level: beginner 418 419 Notes: 420 The array count is allocated internally and must be free'd by the user. 421 422 .seealso: DMSwarmSetType(), DMSwarmSetCellDM(), DMSwarmType 423 @*/ 424 /* 425 PETSC_EXTERN PetscErrorCode DMSwarmCreatePointPerCellCount(DM dm,PetscInt *ncells,PetscInt **count) 426 { 427 PetscFunctionBegin; 428 PetscFunctionReturn(0); 429 } 430 */ 431