1 #define PETSCDM_DLL 2 #include <petsc/private/dmswarmimpl.h> /*I "petscdmswarm.h" I*/ 3 #include <petscsf.h> 4 #include <petscdmda.h> 5 #include <petscdmplex.h> 6 #include "../src/dm/impls/swarm/data_bucket.h" 7 8 /* 9 Error checking to ensure the swarm type is correct and that a cell DM has been set 10 */ 11 #define DMSWARMPICVALID(dm) \ 12 { \ 13 DM_Swarm *_swarm = (DM_Swarm*)(dm)->data; \ 14 if (_swarm->swarm_type != DMSWARM_PIC) SETERRQ(PetscObjectComm((PetscObject)(dm)),PETSC_ERR_SUP,"Valid only for DMSwarm-PIC. You must call DMSwarmSetType(dm,DMSWARM_PIC)"); \ 15 else \ 16 if (!_swarm->dmcell) SETERRQ(PetscObjectComm((PetscObject)(dm)),PETSC_ERR_SUP,"Valid only for DMSwarmPIC if the cell DM is set. You must call DMSwarmSetCellDM(dm,celldm)"); \ 17 } 18 19 /* Coordinate insertition/addition API */ 20 /*@C 21 DMSwarmSetPointsUniformCoordinates - Set point coordinates in a DMSwarm on a regular (ijk) grid 22 23 Collective on dm 24 25 Input parameters: 26 + dm - the DMSwarm 27 . min - minimum coordinate values in the x, y, z directions (array of length dim) 28 . max - maximum coordinate values in the x, y, z directions (array of length dim) 29 . npoints - number of points in each spatial direction (array of length dim) 30 - mode - indicates whether to append points to the swarm (ADD_VALUES), or over-ride existing points (INSERT_VALUES) 31 32 Level: beginner 33 34 Notes: 35 When using mode = INSERT_VALUES, this method will reset the number of particles in the DMSwarm 36 to be npoints[0]*npoints[1] (2D) or npoints[0]*npoints[1]*npoints[2] (3D). When using mode = ADD_VALUES, 37 new points will be appended to any already existing in the DMSwarm 38 39 .seealso: DMSwarmSetType(), DMSwarmSetCellDM(), DMSwarmType 40 @*/ 41 PETSC_EXTERN PetscErrorCode DMSwarmSetPointsUniformCoordinates(DM dm,PetscReal min[],PetscReal max[],PetscInt npoints[],InsertMode mode) 42 { 43 PetscErrorCode ierr; 44 PetscReal gmin[] = {PETSC_MAX_REAL ,PETSC_MAX_REAL, PETSC_MAX_REAL}; 45 PetscReal gmax[] = {PETSC_MIN_REAL, PETSC_MIN_REAL, PETSC_MIN_REAL}; 46 PetscInt i,j,k,N,bs,b,n_estimate,n_curr,n_new_est,p,n_found; 47 Vec coorlocal; 48 const PetscScalar *_coor; 49 DM celldm; 50 PetscReal dx[3]; 51 PetscInt _npoints[] = { 0, 0, 1 }; 52 Vec pos; 53 PetscScalar *_pos; 54 PetscReal *swarm_coor; 55 PetscInt *swarm_cellid; 56 PetscSF sfcell = NULL; 57 const PetscSFNode *LA_sfcell; 58 59 PetscFunctionBegin; 60 DMSWARMPICVALID(dm); 61 ierr = DMSwarmGetCellDM(dm,&celldm);CHKERRQ(ierr); 62 ierr = DMGetCoordinatesLocal(celldm,&coorlocal);CHKERRQ(ierr); 63 ierr = VecGetSize(coorlocal,&N);CHKERRQ(ierr); 64 ierr = VecGetBlockSize(coorlocal,&bs);CHKERRQ(ierr); 65 N = N / bs; 66 ierr = VecGetArrayRead(coorlocal,&_coor);CHKERRQ(ierr); 67 for (i=0; i<N; i++) { 68 for (b=0; b<bs; b++) { 69 gmin[b] = PetscMin(gmin[b],PetscRealPart(_coor[bs*i+b])); 70 gmax[b] = PetscMax(gmax[b],PetscRealPart(_coor[bs*i+b])); 71 } 72 } 73 ierr = VecRestoreArrayRead(coorlocal,&_coor);CHKERRQ(ierr); 74 75 for (b=0; b<bs; b++) { 76 if (npoints[b] > 1) { 77 dx[b] = (max[b] - min[b])/((PetscReal)(npoints[b]-1)); 78 } else { 79 dx[b] = 0.0; 80 } 81 _npoints[b] = npoints[b]; 82 } 83 84 /* determine number of points living in the bounding box */ 85 n_estimate = 0; 86 for (k=0; k<_npoints[2]; k++) { 87 for (j=0; j<_npoints[1]; j++) { 88 for (i=0; i<_npoints[0]; i++) { 89 PetscReal xp[] = {0.0,0.0,0.0}; 90 PetscInt ijk[3]; 91 PetscBool point_inside = PETSC_TRUE; 92 93 ijk[0] = i; 94 ijk[1] = j; 95 ijk[2] = k; 96 for (b=0; b<bs; b++) { 97 xp[b] = min[b] + ijk[b] * dx[b]; 98 } 99 for (b=0; b<bs; b++) { 100 if (xp[b] < gmin[b]) { point_inside = PETSC_FALSE; } 101 if (xp[b] > gmax[b]) { point_inside = PETSC_FALSE; } 102 } 103 if (point_inside) { n_estimate++; } 104 } 105 } 106 } 107 108 /* create candidate list */ 109 ierr = VecCreate(PETSC_COMM_SELF,&pos);CHKERRQ(ierr); 110 ierr = VecSetSizes(pos,bs*n_estimate,PETSC_DECIDE);CHKERRQ(ierr); 111 ierr = VecSetBlockSize(pos,bs);CHKERRQ(ierr); 112 ierr = VecSetFromOptions(pos);CHKERRQ(ierr); 113 ierr = VecGetArray(pos,&_pos);CHKERRQ(ierr); 114 115 n_estimate = 0; 116 for (k=0; k<_npoints[2]; k++) { 117 for (j=0; j<_npoints[1]; j++) { 118 for (i=0; i<_npoints[0]; i++) { 119 PetscReal xp[] = {0.0,0.0,0.0}; 120 PetscInt ijk[3]; 121 PetscBool point_inside = PETSC_TRUE; 122 123 ijk[0] = i; 124 ijk[1] = j; 125 ijk[2] = k; 126 for (b=0; b<bs; b++) { 127 xp[b] = min[b] + ijk[b] * dx[b]; 128 } 129 for (b=0; b<bs; b++) { 130 if (xp[b] < gmin[b]) { point_inside = PETSC_FALSE; } 131 if (xp[b] > gmax[b]) { point_inside = PETSC_FALSE; } 132 } 133 if (point_inside) { 134 for (b=0; b<bs; b++) { 135 _pos[bs*n_estimate+b] = xp[b]; 136 } 137 n_estimate++; 138 } 139 } 140 } 141 } 142 ierr = VecRestoreArray(pos,&_pos);CHKERRQ(ierr); 143 144 /* locate points */ 145 ierr = DMLocatePoints(celldm,pos,DM_POINTLOCATION_NONE,&sfcell);CHKERRQ(ierr); 146 ierr = PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell);CHKERRQ(ierr); 147 n_found = 0; 148 for (p=0; p<n_estimate; p++) { 149 if (LA_sfcell[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) { 150 n_found++; 151 } 152 } 153 154 /* adjust size */ 155 if (mode == ADD_VALUES) { 156 ierr = DMSwarmGetLocalSize(dm,&n_curr);CHKERRQ(ierr); 157 n_new_est = n_curr + n_found; 158 ierr = DMSwarmSetLocalSizes(dm,n_new_est,-1);CHKERRQ(ierr); 159 } 160 if (mode == INSERT_VALUES) { 161 n_curr = 0; 162 n_new_est = n_found; 163 ierr = DMSwarmSetLocalSizes(dm,n_new_est,-1);CHKERRQ(ierr); 164 } 165 166 /* initialize new coords, cell owners, pid */ 167 ierr = VecGetArrayRead(pos,&_coor);CHKERRQ(ierr); 168 ierr = DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr); 169 ierr = DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr); 170 n_found = 0; 171 for (p=0; p<n_estimate; p++) { 172 if (LA_sfcell[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) { 173 for (b=0; b<bs; b++) { 174 swarm_coor[bs*(n_curr + n_found) + b] = PetscRealPart(_coor[bs*p+b]); 175 } 176 swarm_cellid[n_curr + n_found] = LA_sfcell[p].index; 177 n_found++; 178 } 179 } 180 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr); 181 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr); 182 ierr = VecRestoreArrayRead(pos,&_coor);CHKERRQ(ierr); 183 184 ierr = PetscSFDestroy(&sfcell);CHKERRQ(ierr); 185 ierr = VecDestroy(&pos);CHKERRQ(ierr); 186 PetscFunctionReturn(0); 187 } 188 189 /*@C 190 DMSwarmSetPointCoordinates - Set point coordinates in a DMSwarm from a user defined list 191 192 Collective on dm 193 194 Input parameters: 195 + dm - the DMSwarm 196 . npoints - the number of points to insert 197 . coor - the coordinate values 198 . 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 199 - mode - indicates whether to append points to the swarm (ADD_VALUES), or over-ride existing points (INSERT_VALUES) 200 201 Level: beginner 202 203 Notes: 204 If the user has specified redundant = PETSC_FALSE, the cell DM will attempt to locate the coordinates provided by coor[] within 205 its sub-domain. If they any values within coor[] are not located in the sub-domain, they will be ignored and will not get 206 added to the DMSwarm. 207 208 .seealso: DMSwarmSetType(), DMSwarmSetCellDM(), DMSwarmType, DMSwarmSetPointsUniformCoordinates() 209 @*/ 210 PETSC_EXTERN PetscErrorCode DMSwarmSetPointCoordinates(DM dm,PetscInt npoints,PetscReal coor[],PetscBool redundant,InsertMode mode) 211 { 212 PetscErrorCode ierr; 213 PetscReal gmin[] = {PETSC_MAX_REAL ,PETSC_MAX_REAL, PETSC_MAX_REAL}; 214 PetscReal gmax[] = {PETSC_MIN_REAL, PETSC_MIN_REAL, PETSC_MIN_REAL}; 215 PetscInt i,N,bs,b,n_estimate,n_curr,n_new_est,p,n_found; 216 Vec coorlocal; 217 const PetscScalar *_coor; 218 DM celldm; 219 Vec pos; 220 PetscScalar *_pos; 221 PetscReal *swarm_coor; 222 PetscInt *swarm_cellid; 223 PetscSF sfcell = NULL; 224 const PetscSFNode *LA_sfcell; 225 PetscReal *my_coor; 226 PetscInt my_npoints; 227 PetscMPIInt rank; 228 MPI_Comm comm; 229 230 PetscFunctionBegin; 231 DMSWARMPICVALID(dm); 232 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 233 ierr = MPI_Comm_rank(comm,&rank);CHKERRMPI(ierr); 234 235 ierr = DMSwarmGetCellDM(dm,&celldm);CHKERRQ(ierr); 236 ierr = DMGetCoordinatesLocal(celldm,&coorlocal);CHKERRQ(ierr); 237 ierr = VecGetSize(coorlocal,&N);CHKERRQ(ierr); 238 ierr = VecGetBlockSize(coorlocal,&bs);CHKERRQ(ierr); 239 N = N / bs; 240 ierr = VecGetArrayRead(coorlocal,&_coor);CHKERRQ(ierr); 241 for (i=0; i<N; i++) { 242 for (b=0; b<bs; b++) { 243 gmin[b] = PetscMin(gmin[b],PetscRealPart(_coor[bs*i+b])); 244 gmax[b] = PetscMax(gmax[b],PetscRealPart(_coor[bs*i+b])); 245 } 246 } 247 ierr = VecRestoreArrayRead(coorlocal,&_coor);CHKERRQ(ierr); 248 249 /* broadcast points from rank 0 if requested */ 250 if (redundant) { 251 my_npoints = npoints; 252 ierr = MPI_Bcast(&my_npoints,1,MPIU_INT,0,comm);CHKERRMPI(ierr); 253 254 if (rank > 0) { /* allocate space */ 255 ierr = PetscMalloc1(bs*my_npoints,&my_coor);CHKERRQ(ierr); 256 } else { 257 my_coor = coor; 258 } 259 ierr = MPI_Bcast(my_coor,bs*my_npoints,MPIU_REAL,0,comm);CHKERRMPI(ierr); 260 } else { 261 my_npoints = npoints; 262 my_coor = coor; 263 } 264 265 /* determine the number of points living in the bounding box */ 266 n_estimate = 0; 267 for (i=0; i<my_npoints; i++) { 268 PetscBool point_inside = PETSC_TRUE; 269 270 for (b=0; b<bs; b++) { 271 if (my_coor[bs*i+b] < gmin[b]) { point_inside = PETSC_FALSE; } 272 if (my_coor[bs*i+b] > gmax[b]) { point_inside = PETSC_FALSE; } 273 } 274 if (point_inside) { n_estimate++; } 275 } 276 277 /* create candidate list */ 278 ierr = VecCreate(PETSC_COMM_SELF,&pos);CHKERRQ(ierr); 279 ierr = VecSetSizes(pos,bs*n_estimate,PETSC_DECIDE);CHKERRQ(ierr); 280 ierr = VecSetBlockSize(pos,bs);CHKERRQ(ierr); 281 ierr = VecSetFromOptions(pos);CHKERRQ(ierr); 282 ierr = VecGetArray(pos,&_pos);CHKERRQ(ierr); 283 284 n_estimate = 0; 285 for (i=0; i<my_npoints; i++) { 286 PetscBool point_inside = PETSC_TRUE; 287 288 for (b=0; b<bs; b++) { 289 if (my_coor[bs*i+b] < gmin[b]) { point_inside = PETSC_FALSE; } 290 if (my_coor[bs*i+b] > gmax[b]) { point_inside = PETSC_FALSE; } 291 } 292 if (point_inside) { 293 for (b=0; b<bs; b++) { 294 _pos[bs*n_estimate+b] = my_coor[bs*i+b]; 295 } 296 n_estimate++; 297 } 298 } 299 ierr = VecRestoreArray(pos,&_pos);CHKERRQ(ierr); 300 301 /* locate points */ 302 ierr = DMLocatePoints(celldm,pos,DM_POINTLOCATION_NONE,&sfcell);CHKERRQ(ierr); 303 304 ierr = PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell);CHKERRQ(ierr); 305 n_found = 0; 306 for (p=0; p<n_estimate; p++) { 307 if (LA_sfcell[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) { 308 n_found++; 309 } 310 } 311 312 /* adjust size */ 313 if (mode == ADD_VALUES) { 314 ierr = DMSwarmGetLocalSize(dm,&n_curr);CHKERRQ(ierr); 315 n_new_est = n_curr + n_found; 316 ierr = DMSwarmSetLocalSizes(dm,n_new_est,-1);CHKERRQ(ierr); 317 } 318 if (mode == INSERT_VALUES) { 319 n_curr = 0; 320 n_new_est = n_found; 321 ierr = DMSwarmSetLocalSizes(dm,n_new_est,-1);CHKERRQ(ierr); 322 } 323 324 /* initialize new coords, cell owners, pid */ 325 ierr = VecGetArrayRead(pos,&_coor);CHKERRQ(ierr); 326 ierr = DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr); 327 ierr = DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr); 328 n_found = 0; 329 for (p=0; p<n_estimate; p++) { 330 if (LA_sfcell[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) { 331 for (b=0; b<bs; b++) { 332 swarm_coor[bs*(n_curr + n_found) + b] = PetscRealPart(_coor[bs*p+b]); 333 } 334 swarm_cellid[n_curr + n_found] = LA_sfcell[p].index; 335 n_found++; 336 } 337 } 338 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr); 339 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr); 340 ierr = VecRestoreArrayRead(pos,&_coor);CHKERRQ(ierr); 341 342 if (redundant) { 343 if (rank > 0) { 344 ierr = PetscFree(my_coor);CHKERRQ(ierr); 345 } 346 } 347 ierr = PetscSFDestroy(&sfcell);CHKERRQ(ierr); 348 ierr = VecDestroy(&pos);CHKERRQ(ierr); 349 PetscFunctionReturn(0); 350 } 351 352 extern PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_DA(DM,DM,DMSwarmPICLayoutType,PetscInt); 353 extern PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX(DM,DM,DMSwarmPICLayoutType,PetscInt); 354 355 /*@C 356 DMSwarmInsertPointsUsingCellDM - Insert point coordinates within each cell 357 358 Not collective 359 360 Input parameters: 361 + dm - the DMSwarm 362 . layout_type - method used to fill each cell with the cell DM 363 - fill_param - parameter controlling how many points per cell are added (the meaning of this parameter is dependent on the layout type) 364 365 Level: beginner 366 367 Notes: 368 369 The insert method will reset any previous defined points within the DMSwarm. 370 371 When using a DMDA both 2D and 3D are supported for all layout types provided you are using DMDA_ELEMENT_Q1. 372 373 When using a DMPLEX the following case are supported: 374 (i) DMSWARMPIC_LAYOUT_REGULAR: 2D (triangle), 375 (ii) DMSWARMPIC_LAYOUT_GAUSS: 2D and 3D provided the cell is a tri/tet or a quad/hex, 376 (iii) DMSWARMPIC_LAYOUT_SUBDIVISION: 2D and 3D for quad/hex and 2D tri. 377 378 .seealso: DMSwarmPICLayoutType, DMSwarmSetType(), DMSwarmSetCellDM(), DMSwarmType 379 @*/ 380 PETSC_EXTERN PetscErrorCode DMSwarmInsertPointsUsingCellDM(DM dm,DMSwarmPICLayoutType layout_type,PetscInt fill_param) 381 { 382 PetscErrorCode ierr; 383 DM celldm; 384 PetscBool isDA,isPLEX; 385 386 PetscFunctionBegin; 387 DMSWARMPICVALID(dm); 388 ierr = DMSwarmGetCellDM(dm,&celldm);CHKERRQ(ierr); 389 ierr = PetscObjectTypeCompare((PetscObject)celldm,DMDA,&isDA);CHKERRQ(ierr); 390 ierr = PetscObjectTypeCompare((PetscObject)celldm,DMPLEX,&isPLEX);CHKERRQ(ierr); 391 if (isDA) { 392 ierr = private_DMSwarmInsertPointsUsingCellDM_DA(dm,celldm,layout_type,fill_param);CHKERRQ(ierr); 393 } else if (isPLEX) { 394 ierr = private_DMSwarmInsertPointsUsingCellDM_PLEX(dm,celldm,layout_type,fill_param);CHKERRQ(ierr); 395 } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only supported for cell DMs of type DMDA and DMPLEX"); 396 PetscFunctionReturn(0); 397 } 398 399 extern PetscErrorCode private_DMSwarmSetPointCoordinatesCellwise_PLEX(DM,DM,PetscInt,PetscReal*); 400 401 /*@C 402 DMSwarmSetPointCoordinatesCellwise - Insert point coordinates (defined over the reference cell) within each cell 403 404 Not collective 405 406 Input parameters: 407 + dm - the DMSwarm 408 . celldm - the cell DM 409 . npoints - the number of points to insert in each cell 410 - xi - the coordinates (defined in the local coordinate system for each cell) to insert 411 412 Level: beginner 413 414 Notes: 415 The method will reset any previous defined points within the DMSwarm. 416 Only supported for DMPLEX. If you are using a DMDA it is recommended to either use 417 DMSwarmInsertPointsUsingCellDM(), or extract and set the coordinates yourself the following code 418 419 $ PetscReal *coor; 420 $ DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&coor); 421 $ // user code to define the coordinates here 422 $ DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&coor); 423 424 .seealso: DMSwarmSetCellDM(), DMSwarmInsertPointsUsingCellDM() 425 @*/ 426 PETSC_EXTERN PetscErrorCode DMSwarmSetPointCoordinatesCellwise(DM dm,PetscInt npoints,PetscReal xi[]) 427 { 428 PetscErrorCode ierr; 429 DM celldm; 430 PetscBool isDA,isPLEX; 431 432 PetscFunctionBegin; 433 DMSWARMPICVALID(dm); 434 ierr = DMSwarmGetCellDM(dm,&celldm);CHKERRQ(ierr); 435 ierr = PetscObjectTypeCompare((PetscObject)celldm,DMDA,&isDA);CHKERRQ(ierr); 436 ierr = PetscObjectTypeCompare((PetscObject)celldm,DMPLEX,&isPLEX);CHKERRQ(ierr); 437 if (isDA) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only supported for cell DMs of type DMPLEX. Recommended you use DMSwarmInsertPointsUsingCellDM()"); 438 else if (isPLEX) { 439 ierr = private_DMSwarmSetPointCoordinatesCellwise_PLEX(dm,celldm,npoints,xi);CHKERRQ(ierr); 440 } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only supported for cell DMs of type DMDA and DMPLEX"); 441 PetscFunctionReturn(0); 442 } 443 444 /* Field projection API */ 445 extern PetscErrorCode private_DMSwarmProjectFields_DA(DM swarm,DM celldm,PetscInt project_type,PetscInt nfields,DMSwarmDataField dfield[],Vec vecs[]); 446 extern PetscErrorCode private_DMSwarmProjectFields_PLEX(DM swarm,DM celldm,PetscInt project_type,PetscInt nfields,DMSwarmDataField dfield[],Vec vecs[]); 447 448 /*@C 449 DMSwarmProjectFields - Project a set of swarm fields onto the cell DM 450 451 Collective on dm 452 453 Input parameters: 454 + dm - the DMSwarm 455 . nfields - the number of swarm fields to project 456 . fieldnames - the textual names of the swarm fields to project 457 . fields - an array of Vec's of length nfields 458 - reuse - flag indicating whether the array and contents of fields should be re-used or internally allocated 459 460 Currently, the only available projection method consists of 461 phi_i = \sum_{p=0}^{np} N_i(x_p) phi_p dJ / \sum_{p=0}^{np} N_i(x_p) dJ 462 where phi_p is the swarm field at point p, 463 N_i() is the cell DM basis function at vertex i, 464 dJ is the determinant of the cell Jacobian and 465 phi_i is the projected vertex value of the field phi. 466 467 Level: beginner 468 469 Notes: 470 471 If reuse = PETSC_FALSE, this function will allocate the array of Vec's, and each individual Vec. 472 The user is responsible for destroying both the array and the individual Vec objects. 473 474 Only swarm fields registered with data type = PETSC_REAL can be projected onto the cell DM. 475 476 Only swarm fields of block size = 1 can currently be projected. 477 478 The only projection methods currently only support the DA (2D) and PLEX (triangles 2D). 479 480 .seealso: DMSwarmSetType(), DMSwarmSetCellDM(), DMSwarmType 481 @*/ 482 PETSC_EXTERN PetscErrorCode DMSwarmProjectFields(DM dm,PetscInt nfields,const char *fieldnames[],Vec **fields,PetscBool reuse) 483 { 484 DM_Swarm *swarm = (DM_Swarm*)dm->data; 485 DMSwarmDataField *gfield; 486 DM celldm; 487 PetscBool isDA,isPLEX; 488 Vec *vecs; 489 PetscInt f,nvecs; 490 PetscInt project_type = 0; 491 PetscErrorCode ierr; 492 493 PetscFunctionBegin; 494 DMSWARMPICVALID(dm); 495 ierr = DMSwarmGetCellDM(dm,&celldm);CHKERRQ(ierr); 496 ierr = PetscMalloc1(nfields,&gfield);CHKERRQ(ierr); 497 nvecs = 0; 498 for (f=0; f<nfields; f++) { 499 ierr = DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,fieldnames[f],&gfield[f]);CHKERRQ(ierr); 500 if (gfield[f]->petsc_type != PETSC_REAL) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Projection only valid for fields using a data type = PETSC_REAL"); 501 if (gfield[f]->bs != 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Projection only valid for fields with block size = 1"); 502 nvecs += gfield[f]->bs; 503 } 504 if (!reuse) { 505 ierr = PetscMalloc1(nvecs,&vecs);CHKERRQ(ierr); 506 for (f=0; f<nvecs; f++) { 507 ierr = DMCreateGlobalVector(celldm,&vecs[f]);CHKERRQ(ierr); 508 ierr = PetscObjectSetName((PetscObject)vecs[f],gfield[f]->name);CHKERRQ(ierr); 509 } 510 } else { 511 vecs = *fields; 512 } 513 514 ierr = PetscObjectTypeCompare((PetscObject)celldm,DMDA,&isDA);CHKERRQ(ierr); 515 ierr = PetscObjectTypeCompare((PetscObject)celldm,DMPLEX,&isPLEX);CHKERRQ(ierr); 516 if (isDA) { 517 ierr = private_DMSwarmProjectFields_DA(dm,celldm,project_type,nfields,gfield,vecs);CHKERRQ(ierr); 518 } else if (isPLEX) { 519 ierr = private_DMSwarmProjectFields_PLEX(dm,celldm,project_type,nfields,gfield,vecs);CHKERRQ(ierr); 520 } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only supported for cell DMs of type DMDA and DMPLEX"); 521 522 ierr = PetscFree(gfield);CHKERRQ(ierr); 523 if (!reuse) { 524 *fields = vecs; 525 } 526 PetscFunctionReturn(0); 527 } 528 529 /*@C 530 DMSwarmCreatePointPerCellCount - Count the number of points within all cells in the cell DM 531 532 Not collective 533 534 Input parameter: 535 . dm - the DMSwarm 536 537 Output parameters: 538 + ncells - the number of cells in the cell DM (optional argument, pass NULL to ignore) 539 - count - array of length ncells containing the number of points per cell 540 541 Level: beginner 542 543 Notes: 544 The array count is allocated internally and must be free'd by the user. 545 546 .seealso: DMSwarmSetType(), DMSwarmSetCellDM(), DMSwarmType 547 @*/ 548 PETSC_EXTERN PetscErrorCode DMSwarmCreatePointPerCellCount(DM dm,PetscInt *ncells,PetscInt **count) 549 { 550 PetscErrorCode ierr; 551 PetscBool isvalid; 552 PetscInt nel; 553 PetscInt *sum; 554 555 PetscFunctionBegin; 556 ierr = DMSwarmSortGetIsValid(dm,&isvalid);CHKERRQ(ierr); 557 nel = 0; 558 if (isvalid) { 559 PetscInt e; 560 561 ierr = DMSwarmSortGetSizes(dm,&nel,NULL);CHKERRQ(ierr); 562 563 ierr = PetscMalloc1(nel,&sum);CHKERRQ(ierr); 564 for (e=0; e<nel; e++) { 565 ierr = DMSwarmSortGetNumberOfPointsPerCell(dm,e,&sum[e]);CHKERRQ(ierr); 566 } 567 } else { 568 DM celldm; 569 PetscBool isda,isplex,isshell; 570 PetscInt p,npoints; 571 PetscInt *swarm_cellid; 572 573 /* get the number of cells */ 574 ierr = DMSwarmGetCellDM(dm,&celldm);CHKERRQ(ierr); 575 ierr = PetscObjectTypeCompare((PetscObject)celldm,DMDA,&isda);CHKERRQ(ierr); 576 ierr = PetscObjectTypeCompare((PetscObject)celldm,DMPLEX,&isplex);CHKERRQ(ierr); 577 ierr = PetscObjectTypeCompare((PetscObject)celldm,DMSHELL,&isshell);CHKERRQ(ierr); 578 if (isda) { 579 PetscInt _nel,_npe; 580 const PetscInt *_element; 581 582 ierr = DMDAGetElements(celldm,&_nel,&_npe,&_element);CHKERRQ(ierr); 583 nel = _nel; 584 ierr = DMDARestoreElements(celldm,&_nel,&_npe,&_element);CHKERRQ(ierr); 585 } else if (isplex) { 586 PetscInt ps,pe; 587 588 ierr = DMPlexGetHeightStratum(celldm,0,&ps,&pe);CHKERRQ(ierr); 589 nel = pe - ps; 590 } else if (isshell) { 591 PetscErrorCode (*method_DMShellGetNumberOfCells)(DM,PetscInt*); 592 593 ierr = PetscObjectQueryFunction((PetscObject)celldm,"DMGetNumberOfCells_C",&method_DMShellGetNumberOfCells);CHKERRQ(ierr); 594 if (method_DMShellGetNumberOfCells) { 595 ierr = method_DMShellGetNumberOfCells(celldm,&nel);CHKERRQ(ierr); 596 } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Cannot determine the number of cells for the DMSHELL object. User must provide a method via PetscObjectComposeFunction( (PetscObject)shelldm, \"DMGetNumberOfCells_C\", your_function_to_compute_number_of_cells);"); 597 } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Cannot determine the number of cells for a DM not of type DA, PLEX or SHELL"); 598 599 ierr = PetscMalloc1(nel,&sum);CHKERRQ(ierr); 600 ierr = PetscArrayzero(sum,nel);CHKERRQ(ierr); 601 ierr = DMSwarmGetLocalSize(dm,&npoints);CHKERRQ(ierr); 602 ierr = DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr); 603 for (p=0; p<npoints; p++) { 604 if (swarm_cellid[p] != DMLOCATEPOINT_POINT_NOT_FOUND) { 605 sum[ swarm_cellid[p] ]++; 606 } 607 } 608 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr); 609 } 610 if (ncells) { *ncells = nel; } 611 *count = sum; 612 PetscFunctionReturn(0); 613 } 614