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