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