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