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