1 #include <petscdm.h> 2 #include <petscdmplex.h> 3 #include <petscdmswarm.h> 4 #include "../src/dm/impls/swarm/data_bucket.h" 5 6 PetscBool SwarmProjcite = PETSC_FALSE; 7 const char SwarmProjCitation[] = "@article{PusztayKnepleyAdams2022,\n" 8 "title = {Conservative Projection Between FEM and Particle Bases},\n" 9 "author = {Joseph V. Pusztay and Matthew G. Knepley and Mark F. Adams},\n" 10 "journal = {SIAM Journal on Scientific Computing},\n" 11 "volume = {44},\n" 12 "number = {4},\n" 13 "pages = {C310--C319},\n" 14 "doi = {10.1137/21M145407},\n" 15 "year = {2022}\n}\n"; 16 17 PetscErrorCode private_DMSwarmSetPointCoordinatesCellwise_PLEX(DM, DM, PetscInt, PetscReal *xi); 18 19 static PetscErrorCode private_PetscFECreateDefault_scalar_pk1(DM dm, PetscInt dim, PetscBool isSimplex, PetscInt qorder, PetscFE *fem) 20 { 21 const PetscInt Nc = 1; 22 PetscQuadrature q, fq; 23 DM K; 24 PetscSpace P; 25 PetscDualSpace Q; 26 PetscInt order, quadPointsPerEdge; 27 PetscBool tensor = isSimplex ? PETSC_FALSE : PETSC_TRUE; 28 29 PetscFunctionBegin; 30 /* Create space */ 31 PetscCall(PetscSpaceCreate(PetscObjectComm((PetscObject)dm), &P)); 32 /* PetscCall(PetscObjectSetOptionsPrefix((PetscObject) P, prefix)); */ 33 PetscCall(PetscSpacePolynomialSetTensor(P, tensor)); 34 /* PetscCall(PetscSpaceSetFromOptions(P)); */ 35 PetscCall(PetscSpaceSetType(P, PETSCSPACEPOLYNOMIAL)); 36 PetscCall(PetscSpaceSetDegree(P, 1, PETSC_DETERMINE)); 37 PetscCall(PetscSpaceSetNumComponents(P, Nc)); 38 PetscCall(PetscSpaceSetNumVariables(P, dim)); 39 PetscCall(PetscSpaceSetUp(P)); 40 PetscCall(PetscSpaceGetDegree(P, &order, NULL)); 41 PetscCall(PetscSpacePolynomialGetTensor(P, &tensor)); 42 /* Create dual space */ 43 PetscCall(PetscDualSpaceCreate(PetscObjectComm((PetscObject)dm), &Q)); 44 PetscCall(PetscDualSpaceSetType(Q, PETSCDUALSPACELAGRANGE)); 45 /* PetscCall(PetscObjectSetOptionsPrefix((PetscObject) Q, prefix)); */ 46 PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, DMPolytopeTypeSimpleShape(dim, isSimplex), &K)); 47 PetscCall(PetscDualSpaceSetDM(Q, K)); 48 PetscCall(DMDestroy(&K)); 49 PetscCall(PetscDualSpaceSetNumComponents(Q, Nc)); 50 PetscCall(PetscDualSpaceSetOrder(Q, order)); 51 PetscCall(PetscDualSpaceLagrangeSetTensor(Q, tensor)); 52 /* PetscCall(PetscDualSpaceSetFromOptions(Q)); */ 53 PetscCall(PetscDualSpaceSetType(Q, PETSCDUALSPACELAGRANGE)); 54 PetscCall(PetscDualSpaceSetUp(Q)); 55 /* Create element */ 56 PetscCall(PetscFECreate(PetscObjectComm((PetscObject)dm), fem)); 57 /* PetscCall(PetscObjectSetOptionsPrefix((PetscObject) *fem, prefix)); */ 58 /* PetscCall(PetscFESetFromOptions(*fem)); */ 59 PetscCall(PetscFESetType(*fem, PETSCFEBASIC)); 60 PetscCall(PetscFESetBasisSpace(*fem, P)); 61 PetscCall(PetscFESetDualSpace(*fem, Q)); 62 PetscCall(PetscFESetNumComponents(*fem, Nc)); 63 PetscCall(PetscFESetUp(*fem)); 64 PetscCall(PetscSpaceDestroy(&P)); 65 PetscCall(PetscDualSpaceDestroy(&Q)); 66 /* Create quadrature (with specified order if given) */ 67 qorder = qorder >= 0 ? qorder : order; 68 quadPointsPerEdge = PetscMax(qorder + 1, 1); 69 if (isSimplex) { 70 PetscCall(PetscDTStroudConicalQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q)); 71 PetscCall(PetscDTStroudConicalQuadrature(dim - 1, 1, quadPointsPerEdge, -1.0, 1.0, &fq)); 72 } else { 73 PetscCall(PetscDTGaussTensorQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q)); 74 PetscCall(PetscDTGaussTensorQuadrature(dim - 1, 1, quadPointsPerEdge, -1.0, 1.0, &fq)); 75 } 76 PetscCall(PetscFESetQuadrature(*fem, q)); 77 PetscCall(PetscFESetFaceQuadrature(*fem, fq)); 78 PetscCall(PetscQuadratureDestroy(&q)); 79 PetscCall(PetscQuadratureDestroy(&fq)); 80 PetscFunctionReturn(PETSC_SUCCESS); 81 } 82 83 static PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX_SubDivide(DM dm, DM dmc, PetscInt nsub) 84 { 85 PetscInt dim, nfaces, nbasis; 86 PetscInt q, npoints_q, e, nel, pcnt, ps, pe, d, k, r; 87 PetscTabulation T; 88 Vec coorlocal; 89 PetscSection coordSection; 90 PetscScalar *elcoor = NULL; 91 PetscReal *swarm_coor; 92 PetscInt *swarm_cellid; 93 const PetscReal *xiq; 94 PetscQuadrature quadrature; 95 PetscFE fe, feRef; 96 PetscBool is_simplex; 97 98 PetscFunctionBegin; 99 PetscCall(DMGetDimension(dmc, &dim)); 100 is_simplex = PETSC_FALSE; 101 PetscCall(DMPlexGetHeightStratum(dmc, 0, &ps, &pe)); 102 PetscCall(DMPlexGetConeSize(dmc, ps, &nfaces)); 103 if (nfaces == (dim + 1)) is_simplex = PETSC_TRUE; 104 105 PetscCall(private_PetscFECreateDefault_scalar_pk1(dmc, dim, is_simplex, 0, &fe)); 106 107 for (r = 0; r < nsub; r++) { 108 PetscCall(PetscFERefine(fe, &feRef)); 109 PetscCall(PetscFECopyQuadrature(feRef, fe)); 110 PetscCall(PetscFEDestroy(&feRef)); 111 } 112 113 PetscCall(PetscFEGetQuadrature(fe, &quadrature)); 114 PetscCall(PetscQuadratureGetData(quadrature, NULL, NULL, &npoints_q, &xiq, NULL)); 115 PetscCall(PetscFEGetDimension(fe, &nbasis)); 116 PetscCall(PetscFEGetCellTabulation(fe, 1, &T)); 117 118 /* 0->cell, 1->edge, 2->vert */ 119 PetscCall(DMPlexGetHeightStratum(dmc, 0, &ps, &pe)); 120 nel = pe - ps; 121 122 PetscCall(DMSwarmSetLocalSizes(dm, npoints_q * nel, -1)); 123 PetscCall(DMSwarmGetField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor)); 124 PetscCall(DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid)); 125 126 PetscCall(DMGetCoordinatesLocal(dmc, &coorlocal)); 127 PetscCall(DMGetCoordinateSection(dmc, &coordSection)); 128 129 pcnt = 0; 130 for (e = 0; e < nel; e++) { 131 PetscCall(DMPlexVecGetClosure(dmc, coordSection, coorlocal, ps + e, NULL, &elcoor)); 132 133 for (q = 0; q < npoints_q; q++) { 134 for (d = 0; d < dim; d++) { 135 swarm_coor[dim * pcnt + d] = 0.0; 136 for (k = 0; k < nbasis; k++) swarm_coor[dim * pcnt + d] += T->T[0][q * nbasis + k] * PetscRealPart(elcoor[dim * k + d]); 137 } 138 swarm_cellid[pcnt] = e; 139 pcnt++; 140 } 141 PetscCall(DMPlexVecRestoreClosure(dmc, coordSection, coorlocal, ps + e, NULL, &elcoor)); 142 } 143 PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid)); 144 PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor)); 145 146 PetscCall(PetscFEDestroy(&fe)); 147 PetscFunctionReturn(PETSC_SUCCESS); 148 } 149 150 static PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX2D_Regular(DM dm, DM dmc, PetscInt npoints) 151 { 152 PetscInt dim; 153 PetscInt ii, jj, q, npoints_q, e, nel, npe, pcnt, ps, pe, d, k, nfaces; 154 PetscReal *xi, ds, ds2; 155 PetscReal **basis; 156 Vec coorlocal; 157 PetscSection coordSection; 158 PetscScalar *elcoor = NULL; 159 PetscReal *swarm_coor; 160 PetscInt *swarm_cellid; 161 PetscBool is_simplex; 162 163 PetscFunctionBegin; 164 PetscCall(DMGetDimension(dmc, &dim)); 165 PetscCheck(dim == 2, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only 2D is supported"); 166 is_simplex = PETSC_FALSE; 167 PetscCall(DMPlexGetHeightStratum(dmc, 0, &ps, &pe)); 168 PetscCall(DMPlexGetConeSize(dmc, ps, &nfaces)); 169 if (nfaces == (dim + 1)) is_simplex = PETSC_TRUE; 170 PetscCheck(is_simplex, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only the simplex is supported"); 171 172 PetscCall(PetscMalloc1(dim * npoints * npoints, &xi)); 173 pcnt = 0; 174 ds = 1.0 / ((PetscReal)(npoints - 1)); 175 ds2 = 1.0 / ((PetscReal)(npoints)); 176 for (jj = 0; jj < npoints; jj++) { 177 for (ii = 0; ii < npoints - jj; ii++) { 178 xi[dim * pcnt + 0] = ii * ds; 179 xi[dim * pcnt + 1] = jj * ds; 180 181 xi[dim * pcnt + 0] *= (1.0 - 1.2 * ds2); 182 xi[dim * pcnt + 1] *= (1.0 - 1.2 * ds2); 183 184 xi[dim * pcnt + 0] += 0.35 * ds2; 185 xi[dim * pcnt + 1] += 0.35 * ds2; 186 pcnt++; 187 } 188 } 189 npoints_q = pcnt; 190 191 npe = 3; /* nodes per element (triangle) */ 192 PetscCall(PetscMalloc1(npoints_q, &basis)); 193 for (q = 0; q < npoints_q; q++) { 194 PetscCall(PetscMalloc1(npe, &basis[q])); 195 196 basis[q][0] = 1.0 - xi[dim * q + 0] - xi[dim * q + 1]; 197 basis[q][1] = xi[dim * q + 0]; 198 basis[q][2] = xi[dim * q + 1]; 199 } 200 201 /* 0->cell, 1->edge, 2->vert */ 202 PetscCall(DMPlexGetHeightStratum(dmc, 0, &ps, &pe)); 203 nel = pe - ps; 204 205 PetscCall(DMSwarmSetLocalSizes(dm, npoints_q * nel, -1)); 206 PetscCall(DMSwarmGetField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor)); 207 PetscCall(DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid)); 208 209 PetscCall(DMGetCoordinatesLocal(dmc, &coorlocal)); 210 PetscCall(DMGetCoordinateSection(dmc, &coordSection)); 211 212 pcnt = 0; 213 for (e = 0; e < nel; e++) { 214 PetscCall(DMPlexVecGetClosure(dmc, coordSection, coorlocal, e, NULL, &elcoor)); 215 216 for (q = 0; q < npoints_q; q++) { 217 for (d = 0; d < dim; d++) { 218 swarm_coor[dim * pcnt + d] = 0.0; 219 for (k = 0; k < npe; k++) swarm_coor[dim * pcnt + d] += basis[q][k] * PetscRealPart(elcoor[dim * k + d]); 220 } 221 swarm_cellid[pcnt] = e; 222 pcnt++; 223 } 224 PetscCall(DMPlexVecRestoreClosure(dmc, coordSection, coorlocal, e, NULL, &elcoor)); 225 } 226 PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid)); 227 PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor)); 228 229 PetscCall(PetscFree(xi)); 230 for (q = 0; q < npoints_q; q++) PetscCall(PetscFree(basis[q])); 231 PetscCall(PetscFree(basis)); 232 PetscFunctionReturn(PETSC_SUCCESS); 233 } 234 235 PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX(DM dm, DM celldm, DMSwarmPICLayoutType layout, PetscInt layout_param) 236 { 237 PetscInt dim; 238 239 PetscFunctionBegin; 240 PetscCall(DMGetDimension(celldm, &dim)); 241 switch (layout) { 242 case DMSWARMPIC_LAYOUT_REGULAR: 243 PetscCheck(dim != 3, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "No 3D support for REGULAR+PLEX"); 244 PetscCall(private_DMSwarmInsertPointsUsingCellDM_PLEX2D_Regular(dm, celldm, layout_param)); 245 break; 246 case DMSWARMPIC_LAYOUT_GAUSS: { 247 PetscInt npoints, npoints1, ps, pe, nfaces; 248 const PetscReal *xi; 249 PetscBool is_simplex; 250 PetscQuadrature quadrature; 251 252 is_simplex = PETSC_FALSE; 253 PetscCall(DMPlexGetHeightStratum(celldm, 0, &ps, &pe)); 254 PetscCall(DMPlexGetConeSize(celldm, ps, &nfaces)); 255 if (nfaces == (dim + 1)) is_simplex = PETSC_TRUE; 256 257 npoints1 = layout_param; 258 if (is_simplex) { 259 PetscCall(PetscDTStroudConicalQuadrature(dim, 1, npoints1, -1.0, 1.0, &quadrature)); 260 } else { 261 PetscCall(PetscDTGaussTensorQuadrature(dim, 1, npoints1, -1.0, 1.0, &quadrature)); 262 } 263 PetscCall(PetscQuadratureGetData(quadrature, NULL, NULL, &npoints, &xi, NULL)); 264 PetscCall(private_DMSwarmSetPointCoordinatesCellwise_PLEX(dm, celldm, npoints, (PetscReal *)xi)); 265 PetscCall(PetscQuadratureDestroy(&quadrature)); 266 } break; 267 case DMSWARMPIC_LAYOUT_SUBDIVISION: 268 PetscCall(private_DMSwarmInsertPointsUsingCellDM_PLEX_SubDivide(dm, celldm, layout_param)); 269 break; 270 } 271 PetscFunctionReturn(PETSC_SUCCESS); 272 } 273 274 /* 275 typedef struct { 276 PetscReal x,y; 277 } Point2d; 278 279 static PetscErrorCode signp2d(Point2d p1, Point2d p2, Point2d p3,PetscReal *s) 280 { 281 PetscFunctionBegin; 282 *s = (p1.x - p3.x) * (p2.y - p3.y) - (p2.x - p3.x) * (p1.y - p3.y); 283 PetscFunctionReturn(PETSC_SUCCESS); 284 } 285 */ 286 /* 287 static PetscErrorCode PointInTriangle(Point2d pt, Point2d v1, Point2d v2, Point2d v3,PetscBool *v) 288 { 289 PetscReal s1,s2,s3; 290 PetscBool b1, b2, b3; 291 292 PetscFunctionBegin; 293 signp2d(pt, v1, v2,&s1); b1 = s1 < 0.0f; 294 signp2d(pt, v2, v3,&s2); b2 = s2 < 0.0f; 295 signp2d(pt, v3, v1,&s3); b3 = s3 < 0.0f; 296 297 *v = ((b1 == b2) && (b2 == b3)); 298 PetscFunctionReturn(PETSC_SUCCESS); 299 } 300 */ 301 /* 302 static PetscErrorCode _ComputeLocalCoordinateAffine2d(PetscReal xp[],PetscReal coords[],PetscReal xip[],PetscReal *dJ) 303 { 304 PetscReal x1,y1,x2,y2,x3,y3; 305 PetscReal c,b[2],A[2][2],inv[2][2],detJ,od; 306 307 PetscFunctionBegin; 308 x1 = coords[2*0+0]; 309 x2 = coords[2*1+0]; 310 x3 = coords[2*2+0]; 311 312 y1 = coords[2*0+1]; 313 y2 = coords[2*1+1]; 314 y3 = coords[2*2+1]; 315 316 c = x1 - 0.5*x1 - 0.5*x1 + 0.5*x2 + 0.5*x3; 317 b[0] = xp[0] - c; 318 c = y1 - 0.5*y1 - 0.5*y1 + 0.5*y2 + 0.5*y3; 319 b[1] = xp[1] - c; 320 321 A[0][0] = -0.5*x1 + 0.5*x2; A[0][1] = -0.5*x1 + 0.5*x3; 322 A[1][0] = -0.5*y1 + 0.5*y2; A[1][1] = -0.5*y1 + 0.5*y3; 323 324 detJ = A[0][0]*A[1][1] - A[0][1]*A[1][0]; 325 *dJ = PetscAbsReal(detJ); 326 od = 1.0/detJ; 327 328 inv[0][0] = A[1][1] * od; 329 inv[0][1] = -A[0][1] * od; 330 inv[1][0] = -A[1][0] * od; 331 inv[1][1] = A[0][0] * od; 332 333 xip[0] = inv[0][0]*b[0] + inv[0][1]*b[1]; 334 xip[1] = inv[1][0]*b[0] + inv[1][1]*b[1]; 335 PetscFunctionReturn(PETSC_SUCCESS); 336 } 337 */ 338 339 static PetscErrorCode ComputeLocalCoordinateAffine2d(PetscReal xp[], PetscScalar coords[], PetscReal xip[], PetscReal *dJ) 340 { 341 PetscReal x1, y1, x2, y2, x3, y3; 342 PetscReal b[2], A[2][2], inv[2][2], detJ, od; 343 344 PetscFunctionBegin; 345 x1 = PetscRealPart(coords[2 * 0 + 0]); 346 x2 = PetscRealPart(coords[2 * 1 + 0]); 347 x3 = PetscRealPart(coords[2 * 2 + 0]); 348 349 y1 = PetscRealPart(coords[2 * 0 + 1]); 350 y2 = PetscRealPart(coords[2 * 1 + 1]); 351 y3 = PetscRealPart(coords[2 * 2 + 1]); 352 353 b[0] = xp[0] - x1; 354 b[1] = xp[1] - y1; 355 356 A[0][0] = x2 - x1; 357 A[0][1] = x3 - x1; 358 A[1][0] = y2 - y1; 359 A[1][1] = y3 - y1; 360 361 detJ = A[0][0] * A[1][1] - A[0][1] * A[1][0]; 362 *dJ = PetscAbsReal(detJ); 363 od = 1.0 / detJ; 364 365 inv[0][0] = A[1][1] * od; 366 inv[0][1] = -A[0][1] * od; 367 inv[1][0] = -A[1][0] * od; 368 inv[1][1] = A[0][0] * od; 369 370 xip[0] = inv[0][0] * b[0] + inv[0][1] * b[1]; 371 xip[1] = inv[1][0] * b[0] + inv[1][1] * b[1]; 372 PetscFunctionReturn(PETSC_SUCCESS); 373 } 374 375 static PetscErrorCode DMSwarmProjectField_ApproxP1_PLEX_2D(DM swarm, PetscReal *swarm_field, DM dm, Vec v_field) 376 { 377 const PetscReal PLEX_C_EPS = 1.0e-8; 378 Vec v_field_l, denom_l, coor_l, denom; 379 PetscInt k, p, e, npoints; 380 PetscInt *mpfield_cell; 381 PetscReal *mpfield_coor; 382 PetscReal xi_p[2]; 383 PetscScalar Ni[3]; 384 PetscSection coordSection; 385 PetscScalar *elcoor = NULL; 386 387 PetscFunctionBegin; 388 PetscCall(VecZeroEntries(v_field)); 389 390 PetscCall(DMGetLocalVector(dm, &v_field_l)); 391 PetscCall(DMGetGlobalVector(dm, &denom)); 392 PetscCall(DMGetLocalVector(dm, &denom_l)); 393 PetscCall(VecZeroEntries(v_field_l)); 394 PetscCall(VecZeroEntries(denom)); 395 PetscCall(VecZeroEntries(denom_l)); 396 397 PetscCall(DMGetCoordinatesLocal(dm, &coor_l)); 398 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 399 400 PetscCall(DMSwarmGetLocalSize(swarm, &npoints)); 401 PetscCall(DMSwarmGetField(swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&mpfield_coor)); 402 PetscCall(DMSwarmGetField(swarm, DMSwarmPICField_cellid, NULL, NULL, (void **)&mpfield_cell)); 403 404 for (p = 0; p < npoints; p++) { 405 PetscReal *coor_p, dJ; 406 PetscScalar elfield[3]; 407 PetscBool point_located; 408 409 e = mpfield_cell[p]; 410 coor_p = &mpfield_coor[2 * p]; 411 412 PetscCall(DMPlexVecGetClosure(dm, coordSection, coor_l, e, NULL, &elcoor)); 413 414 /* 415 while (!point_located && (failed_counter < 25)) { 416 PetscCall(PointInTriangle(point, coords[0], coords[1], coords[2], &point_located)); 417 point.x = coor_p[0]; 418 point.y = coor_p[1]; 419 point.x += 1.0e-10 * (2.0 * rand()/((double)RAND_MAX)-1.0); 420 point.y += 1.0e-10 * (2.0 * rand()/((double)RAND_MAX)-1.0); 421 failed_counter++; 422 } 423 424 if (!point_located) { 425 PetscPrintf(PETSC_COMM_SELF,"Failed to locate point (%1.8e,%1.8e) in local mesh (cell %" PetscInt_FMT ") with triangle coords (%1.8e,%1.8e) : (%1.8e,%1.8e) : (%1.8e,%1.8e) in %" PetscInt_FMT " iterations\n",point.x,point.y,e,coords[0].x,coords[0].y,coords[1].x,coords[1].y,coords[2].x,coords[2].y,failed_counter); 426 } 427 428 PetscCheck(point_located,PETSC_COMM_SELF,PETSC_ERR_SUP,"Failed to locate point (%1.8e,%1.8e) in local mesh (cell %" PetscInt_FMT ")",point.x,point.y,e); 429 else { 430 PetscCall(_ComputeLocalCoordinateAffine2d(coor_p,elcoor,xi_p,&dJ)); 431 xi_p[0] = 0.5*(xi_p[0] + 1.0); 432 xi_p[1] = 0.5*(xi_p[1] + 1.0); 433 434 PetscPrintf(PETSC_COMM_SELF,"[p=%" PetscInt_FMT "] x(%+1.4e,%+1.4e) -> mapped to element %" PetscInt_FMT " xi(%+1.4e,%+1.4e)\n",p,point.x,point.y,e,xi_p[0],xi_p[1]); 435 436 } 437 */ 438 439 PetscCall(ComputeLocalCoordinateAffine2d(coor_p, elcoor, xi_p, &dJ)); 440 /* 441 PetscPrintf(PETSC_COMM_SELF,"[p=%" PetscInt_FMT "] x(%+1.4e,%+1.4e) -> mapped to element %" PetscInt_FMT " xi(%+1.4e,%+1.4e)\n",p,point.x,point.y,e,xi_p[0],xi_p[1]); 442 */ 443 /* 444 point_located = PETSC_TRUE; 445 if (xi_p[0] < 0.0) { 446 if (xi_p[0] > -PLEX_C_EPS) { 447 xi_p[0] = 0.0; 448 } else { 449 point_located = PETSC_FALSE; 450 } 451 } 452 if (xi_p[1] < 0.0) { 453 if (xi_p[1] > -PLEX_C_EPS) { 454 xi_p[1] = 0.0; 455 } else { 456 point_located = PETSC_FALSE; 457 } 458 } 459 if (xi_p[1] > (1.0-xi_p[0])) { 460 if ((xi_p[1] - 1.0 + xi_p[0]) < PLEX_C_EPS) { 461 xi_p[1] = 1.0 - xi_p[0]; 462 } else { 463 point_located = PETSC_FALSE; 464 } 465 } 466 if (!point_located) { 467 PetscPrintf(PETSC_COMM_SELF,"[Error] xi,eta = %+1.8e, %+1.8e\n",xi_p[0],xi_p[1]); 468 PetscPrintf(PETSC_COMM_SELF,"[Error] Failed to locate point (%1.8e,%1.8e) in local mesh (cell %" PetscInt_FMT ") with triangle coords (%1.8e,%1.8e) : (%1.8e,%1.8e) : (%1.8e,%1.8e)\n",coor_p[0],coor_p[1],e,elcoor[0],elcoor[1],elcoor[2],elcoor[3],elcoor[4],elcoor[5]); 469 } 470 PetscCheck(point_located,PETSC_COMM_SELF,PETSC_ERR_SUP,"Failed to locate point (%1.8e,%1.8e) in local mesh (cell %" PetscInt_FMT ")",coor_p[0],coor_p[1],e); 471 */ 472 473 Ni[0] = 1.0 - xi_p[0] - xi_p[1]; 474 Ni[1] = xi_p[0]; 475 Ni[2] = xi_p[1]; 476 477 point_located = PETSC_TRUE; 478 for (k = 0; k < 3; k++) { 479 if (PetscRealPart(Ni[k]) < -PLEX_C_EPS) point_located = PETSC_FALSE; 480 if (PetscRealPart(Ni[k]) > (1.0 + PLEX_C_EPS)) point_located = PETSC_FALSE; 481 } 482 if (!point_located) { 483 PetscCall(PetscPrintf(PETSC_COMM_SELF, "[Error] xi,eta = %+1.8e, %+1.8e\n", (double)xi_p[0], (double)xi_p[1])); 484 PetscCall(PetscPrintf(PETSC_COMM_SELF, "[Error] Failed to locate point (%1.8e,%1.8e) in local mesh (cell %" PetscInt_FMT ") with triangle coords (%1.8e,%1.8e) : (%1.8e,%1.8e) : (%1.8e,%1.8e)\n", (double)coor_p[0], (double)coor_p[1], e, (double)PetscRealPart(elcoor[0]), (double)PetscRealPart(elcoor[1]), (double)PetscRealPart(elcoor[2]), (double)PetscRealPart(elcoor[3]), (double)PetscRealPart(elcoor[4]), (double)PetscRealPart(elcoor[5]))); 485 } 486 PetscCheck(point_located, PETSC_COMM_SELF, PETSC_ERR_SUP, "Failed to locate point (%1.8e,%1.8e) in local mesh (cell %" PetscInt_FMT ")", (double)coor_p[0], (double)coor_p[1], e); 487 488 for (k = 0; k < 3; k++) { 489 Ni[k] = Ni[k] * dJ; 490 elfield[k] = Ni[k] * swarm_field[p]; 491 } 492 PetscCall(DMPlexVecRestoreClosure(dm, coordSection, coor_l, e, NULL, &elcoor)); 493 494 PetscCall(DMPlexVecSetClosure(dm, NULL, v_field_l, e, elfield, ADD_VALUES)); 495 PetscCall(DMPlexVecSetClosure(dm, NULL, denom_l, e, Ni, ADD_VALUES)); 496 } 497 498 PetscCall(DMSwarmRestoreField(swarm, DMSwarmPICField_cellid, NULL, NULL, (void **)&mpfield_cell)); 499 PetscCall(DMSwarmRestoreField(swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&mpfield_coor)); 500 501 PetscCall(DMLocalToGlobalBegin(dm, v_field_l, ADD_VALUES, v_field)); 502 PetscCall(DMLocalToGlobalEnd(dm, v_field_l, ADD_VALUES, v_field)); 503 PetscCall(DMLocalToGlobalBegin(dm, denom_l, ADD_VALUES, denom)); 504 PetscCall(DMLocalToGlobalEnd(dm, denom_l, ADD_VALUES, denom)); 505 506 PetscCall(VecPointwiseDivide(v_field, v_field, denom)); 507 508 PetscCall(DMRestoreLocalVector(dm, &v_field_l)); 509 PetscCall(DMRestoreLocalVector(dm, &denom_l)); 510 PetscCall(DMRestoreGlobalVector(dm, &denom)); 511 PetscFunctionReturn(PETSC_SUCCESS); 512 } 513 514 PetscErrorCode private_DMSwarmProjectFields_PLEX(DM swarm, DM celldm, PetscInt project_type, PetscInt nfields, DMSwarmDataField dfield[], Vec vecs[]) 515 { 516 PetscInt f, dim; 517 518 PetscFunctionBegin; 519 PetscCall(DMGetDimension(swarm, &dim)); 520 switch (dim) { 521 case 2: 522 for (f = 0; f < nfields; f++) { 523 PetscReal *swarm_field; 524 525 PetscCall(DMSwarmDataFieldGetEntries(dfield[f], (void **)&swarm_field)); 526 PetscCall(DMSwarmProjectField_ApproxP1_PLEX_2D(swarm, swarm_field, celldm, vecs[f])); 527 } 528 break; 529 case 3: 530 SETERRQ(PetscObjectComm((PetscObject)swarm), PETSC_ERR_SUP, "No support for 3D"); 531 default: 532 break; 533 } 534 PetscFunctionReturn(PETSC_SUCCESS); 535 } 536 537 PetscErrorCode private_DMSwarmSetPointCoordinatesCellwise_PLEX(DM dm, DM dmc, PetscInt npoints, PetscReal xi[]) 538 { 539 PetscBool is_simplex, is_tensorcell; 540 PetscInt dim, nfaces, ps, pe, p, d, nbasis, pcnt, e, k, nel; 541 PetscFE fe; 542 PetscQuadrature quadrature; 543 PetscTabulation T; 544 PetscReal *xiq; 545 Vec coorlocal; 546 PetscSection coordSection; 547 PetscScalar *elcoor = NULL; 548 PetscReal *swarm_coor; 549 PetscInt *swarm_cellid; 550 551 PetscFunctionBegin; 552 PetscCall(DMGetDimension(dmc, &dim)); 553 554 is_simplex = PETSC_FALSE; 555 is_tensorcell = PETSC_FALSE; 556 PetscCall(DMPlexGetHeightStratum(dmc, 0, &ps, &pe)); 557 PetscCall(DMPlexGetConeSize(dmc, ps, &nfaces)); 558 559 if (nfaces == (dim + 1)) is_simplex = PETSC_TRUE; 560 561 switch (dim) { 562 case 2: 563 if (nfaces == 4) is_tensorcell = PETSC_TRUE; 564 break; 565 case 3: 566 if (nfaces == 6) is_tensorcell = PETSC_TRUE; 567 break; 568 default: 569 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only support for 2D, 3D"); 570 } 571 572 /* check points provided fail inside the reference cell */ 573 if (is_simplex) { 574 for (p = 0; p < npoints; p++) { 575 PetscReal sum; 576 for (d = 0; d < dim; d++) PetscCheck(xi[dim * p + d] >= -1.0, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Points do not fail inside the simplex domain"); 577 sum = 0.0; 578 for (d = 0; d < dim; d++) sum += xi[dim * p + d]; 579 PetscCheck(sum <= 0.0, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Points do not fail inside the simplex domain"); 580 } 581 } else if (is_tensorcell) { 582 for (p = 0; p < npoints; p++) { 583 for (d = 0; d < dim; d++) PetscCheck(PetscAbsReal(xi[dim * p + d]) <= 1.0, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Points do not fail inside the tensor domain [-1,1]^d"); 584 } 585 } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only support for d-simplex and d-tensorcell"); 586 587 PetscCall(PetscQuadratureCreate(PetscObjectComm((PetscObject)dm), &quadrature)); 588 PetscCall(PetscMalloc1(npoints * dim, &xiq)); 589 PetscCall(PetscArraycpy(xiq, xi, npoints * dim)); 590 PetscCall(PetscQuadratureSetData(quadrature, dim, 1, npoints, (const PetscReal *)xiq, NULL)); 591 PetscCall(private_PetscFECreateDefault_scalar_pk1(dmc, dim, is_simplex, 0, &fe)); 592 PetscCall(PetscFESetQuadrature(fe, quadrature)); 593 PetscCall(PetscFEGetDimension(fe, &nbasis)); 594 PetscCall(PetscFEGetCellTabulation(fe, 1, &T)); 595 596 /* for each cell, interpolate coordaintes and insert the interpolated points coordinates into swarm */ 597 /* 0->cell, 1->edge, 2->vert */ 598 PetscCall(DMPlexGetHeightStratum(dmc, 0, &ps, &pe)); 599 nel = pe - ps; 600 601 PetscCall(DMSwarmSetLocalSizes(dm, npoints * nel, -1)); 602 PetscCall(DMSwarmGetField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor)); 603 PetscCall(DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid)); 604 605 PetscCall(DMGetCoordinatesLocal(dmc, &coorlocal)); 606 PetscCall(DMGetCoordinateSection(dmc, &coordSection)); 607 608 pcnt = 0; 609 for (e = 0; e < nel; e++) { 610 PetscCall(DMPlexVecGetClosure(dmc, coordSection, coorlocal, ps + e, NULL, &elcoor)); 611 612 for (p = 0; p < npoints; p++) { 613 for (d = 0; d < dim; d++) { 614 swarm_coor[dim * pcnt + d] = 0.0; 615 for (k = 0; k < nbasis; k++) swarm_coor[dim * pcnt + d] += T->T[0][p * nbasis + k] * PetscRealPart(elcoor[dim * k + d]); 616 } 617 swarm_cellid[pcnt] = e; 618 pcnt++; 619 } 620 PetscCall(DMPlexVecRestoreClosure(dmc, coordSection, coorlocal, ps + e, NULL, &elcoor)); 621 } 622 PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid)); 623 PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor)); 624 625 PetscCall(PetscQuadratureDestroy(&quadrature)); 626 PetscCall(PetscFEDestroy(&fe)); 627 PetscFunctionReturn(PETSC_SUCCESS); 628 } 629