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