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