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