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 PetscInt dim; 263 PetscInt ii,jj,q,npoints_q,e,nel,npe,pcnt,ps,pe,d,k,nfaces; 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 PetscBool is_simplex; 272 273 PetscFunctionBegin; 274 ierr = DMGetDimension(dmc,&dim);CHKERRQ(ierr); 275 if (dim != 2) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only 2D is supported"); 276 is_simplex = PETSC_FALSE; 277 ierr = DMPlexGetHeightStratum(dmc,0,&ps,&pe);CHKERRQ(ierr); 278 ierr = DMPlexGetConeSize(dmc, ps, &nfaces);CHKERRQ(ierr); 279 if (nfaces == (dim+1)) { is_simplex = PETSC_TRUE; } 280 if (!is_simplex) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only the simplex is supported"); 281 282 ierr = PetscMalloc1(dim*npoints*npoints,&xi);CHKERRQ(ierr); 283 pcnt = 0; 284 ds = 1.0/((PetscReal)(npoints-1)); 285 ds2 = 1.0/((PetscReal)(npoints)); 286 for (jj = 0; jj<npoints; jj++) { 287 for (ii=0; ii<npoints-jj; ii++) { 288 xi[dim*pcnt+0] = ii * ds; 289 xi[dim*pcnt+1] = jj * ds; 290 291 xi[dim*pcnt+0] *= (1.0 - 1.2*ds2); 292 xi[dim*pcnt+1] *= (1.0 - 1.2*ds2); 293 294 xi[dim*pcnt+0] += 0.35*ds2; 295 xi[dim*pcnt+1] += 0.35*ds2; 296 297 pcnt++; 298 } 299 } 300 npoints_q = pcnt; 301 302 npe = 3; /* nodes per element (triangle) */ 303 ierr = PetscMalloc1(npoints_q,&basis);CHKERRQ(ierr); 304 for (q=0; q<npoints_q; q++) { 305 ierr = PetscMalloc1(npe,&basis[q]);CHKERRQ(ierr); 306 307 basis[q][0] = 1.0 - xi[dim*q+0] - xi[dim*q+1]; 308 basis[q][1] = xi[dim*q+0]; 309 basis[q][2] = xi[dim*q+1]; 310 } 311 312 /* 0->cell, 1->edge, 2->vert */ 313 ierr = DMPlexGetHeightStratum(dmc,0,&ps,&pe);CHKERRQ(ierr); 314 nel = pe - ps; 315 316 ierr = DMSwarmSetLocalSizes(dm,npoints_q*nel,-1);CHKERRQ(ierr); 317 ierr = DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr); 318 ierr = DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr); 319 320 ierr = DMGetCoordinatesLocal(dmc,&coorlocal);CHKERRQ(ierr); 321 ierr = DMGetCoordinateSection(dmc,&coordSection);CHKERRQ(ierr); 322 323 pcnt = 0; 324 for (e=0; e<nel; e++) { 325 ierr = DMPlexVecGetClosure(dmc,coordSection,coorlocal,e,NULL,&elcoor);CHKERRQ(ierr); 326 327 for (q=0; q<npoints_q; q++) { 328 for (d=0; d<dim; d++) { 329 swarm_coor[dim*pcnt+d] = 0.0; 330 for (k=0; k<npe; k++) { 331 swarm_coor[dim*pcnt+d] += basis[q][k] * PetscRealPart(elcoor[dim*k+d]); 332 } 333 } 334 swarm_cellid[pcnt] = e; 335 pcnt++; 336 } 337 ierr = DMPlexVecRestoreClosure(dmc,coordSection,coorlocal,e,NULL,&elcoor);CHKERRQ(ierr); 338 } 339 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr); 340 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr); 341 342 ierr = PetscFree(xi);CHKERRQ(ierr); 343 for (q=0; q<npoints_q; q++) { 344 ierr = PetscFree(basis[q]);CHKERRQ(ierr); 345 } 346 ierr = PetscFree(basis);CHKERRQ(ierr); 347 348 PetscFunctionReturn(0); 349 } 350 351 PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX(DM dm,DM celldm,DMSwarmPICLayoutType layout,PetscInt layout_param) 352 { 353 PetscErrorCode ierr; 354 PetscInt dim; 355 356 PetscFunctionBegin; 357 ierr = DMGetDimension(celldm,&dim);CHKERRQ(ierr); 358 switch (layout) { 359 case DMSWARMPIC_LAYOUT_REGULAR: 360 if (dim == 3) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No 3D support for REGULAR+PLEX"); 361 ierr = private_DMSwarmInsertPointsUsingCellDM_PLEX2D_Regular(dm,celldm,layout_param);CHKERRQ(ierr); 362 break; 363 case DMSWARMPIC_LAYOUT_GAUSS: 364 { 365 PetscInt npoints,npoints1,ps,pe,nfaces; 366 const PetscReal *xi; 367 PetscBool is_simplex; 368 PetscQuadrature quadrature; 369 370 is_simplex = PETSC_FALSE; 371 ierr = DMPlexGetHeightStratum(celldm,0,&ps,&pe);CHKERRQ(ierr); 372 ierr = DMPlexGetConeSize(celldm, ps, &nfaces);CHKERRQ(ierr); 373 if (nfaces == (dim+1)) { is_simplex = PETSC_TRUE; } 374 375 npoints1 = layout_param; 376 if (is_simplex) { 377 ierr = PetscDTGaussJacobiQuadrature(dim,1,npoints1,-1.0,1.0,&quadrature);CHKERRQ(ierr); 378 } else { 379 ierr = PetscDTGaussTensorQuadrature(dim,1,npoints1,-1.0,1.0,&quadrature);CHKERRQ(ierr); 380 } 381 ierr = PetscQuadratureGetData(quadrature,NULL,NULL,&npoints,&xi,NULL);CHKERRQ(ierr); 382 ierr = private_DMSwarmSetPointCoordinatesCellwise_PLEX(dm,celldm,npoints,(PetscReal*)xi);CHKERRQ(ierr); 383 ierr = PetscQuadratureDestroy(&quadrature);CHKERRQ(ierr); 384 } 385 break; 386 case DMSWARMPIC_LAYOUT_SUBDIVISION: 387 ierr = private_DMSwarmInsertPointsUsingCellDM_PLEX_SubDivide(dm,celldm,layout_param);CHKERRQ(ierr); 388 break; 389 } 390 PetscFunctionReturn(0); 391 } 392 393 /* 394 typedef struct { 395 PetscReal x,y; 396 } Point2d; 397 398 static PetscErrorCode signp2d(Point2d p1, Point2d p2, Point2d p3,PetscReal *s) 399 { 400 PetscFunctionBegin; 401 *s = (p1.x - p3.x) * (p2.y - p3.y) - (p2.x - p3.x) * (p1.y - p3.y); 402 PetscFunctionReturn(0); 403 } 404 */ 405 /* 406 static PetscErrorCode PointInTriangle(Point2d pt, Point2d v1, Point2d v2, Point2d v3,PetscBool *v) 407 { 408 PetscReal s1,s2,s3; 409 PetscBool b1, b2, b3; 410 411 PetscFunctionBegin; 412 signp2d(pt, v1, v2,&s1); b1 = s1 < 0.0f; 413 signp2d(pt, v2, v3,&s2); b2 = s2 < 0.0f; 414 signp2d(pt, v3, v1,&s3); b3 = s3 < 0.0f; 415 416 *v = ((b1 == b2) && (b2 == b3)); 417 PetscFunctionReturn(0); 418 } 419 */ 420 /* 421 static PetscErrorCode _ComputeLocalCoordinateAffine2d(PetscReal xp[],PetscReal coords[],PetscReal xip[],PetscReal *dJ) 422 { 423 PetscReal x1,y1,x2,y2,x3,y3; 424 PetscReal c,b[2],A[2][2],inv[2][2],detJ,od; 425 426 PetscFunctionBegin; 427 x1 = coords[2*0+0]; 428 x2 = coords[2*1+0]; 429 x3 = coords[2*2+0]; 430 431 y1 = coords[2*0+1]; 432 y2 = coords[2*1+1]; 433 y3 = coords[2*2+1]; 434 435 c = x1 - 0.5*x1 - 0.5*x1 + 0.5*x2 + 0.5*x3; 436 b[0] = xp[0] - c; 437 c = y1 - 0.5*y1 - 0.5*y1 + 0.5*y2 + 0.5*y3; 438 b[1] = xp[1] - c; 439 440 A[0][0] = -0.5*x1 + 0.5*x2; A[0][1] = -0.5*x1 + 0.5*x3; 441 A[1][0] = -0.5*y1 + 0.5*y2; A[1][1] = -0.5*y1 + 0.5*y3; 442 443 detJ = A[0][0]*A[1][1] - A[0][1]*A[1][0]; 444 *dJ = PetscAbsReal(detJ); 445 od = 1.0/detJ; 446 447 inv[0][0] = A[1][1] * od; 448 inv[0][1] = -A[0][1] * od; 449 inv[1][0] = -A[1][0] * od; 450 inv[1][1] = A[0][0] * od; 451 452 xip[0] = inv[0][0]*b[0] + inv[0][1]*b[1]; 453 xip[1] = inv[1][0]*b[0] + inv[1][1]*b[1]; 454 PetscFunctionReturn(0); 455 } 456 */ 457 458 static PetscErrorCode ComputeLocalCoordinateAffine2d(PetscReal xp[],PetscScalar coords[],PetscReal xip[],PetscReal *dJ) 459 { 460 PetscReal x1,y1,x2,y2,x3,y3; 461 PetscReal b[2],A[2][2],inv[2][2],detJ,od; 462 463 PetscFunctionBegin; 464 x1 = PetscRealPart(coords[2*0+0]); 465 x2 = PetscRealPart(coords[2*1+0]); 466 x3 = PetscRealPart(coords[2*2+0]); 467 468 y1 = PetscRealPart(coords[2*0+1]); 469 y2 = PetscRealPart(coords[2*1+1]); 470 y3 = PetscRealPart(coords[2*2+1]); 471 472 b[0] = xp[0] - x1; 473 b[1] = xp[1] - y1; 474 475 A[0][0] = x2-x1; A[0][1] = x3-x1; 476 A[1][0] = y2-y1; A[1][1] = y3-y1; 477 478 detJ = A[0][0]*A[1][1] - A[0][1]*A[1][0]; 479 *dJ = PetscAbsReal(detJ); 480 od = 1.0/detJ; 481 482 inv[0][0] = A[1][1] * od; 483 inv[0][1] = -A[0][1] * od; 484 inv[1][0] = -A[1][0] * od; 485 inv[1][1] = A[0][0] * od; 486 487 xip[0] = inv[0][0]*b[0] + inv[0][1]*b[1]; 488 xip[1] = inv[1][0]*b[0] + inv[1][1]*b[1]; 489 PetscFunctionReturn(0); 490 } 491 492 PetscErrorCode DMSwarmProjectField_ApproxP1_PLEX_2D(DM swarm,PetscReal *swarm_field,DM dm,Vec v_field) 493 { 494 PetscErrorCode ierr; 495 const PetscReal PLEX_C_EPS = 1.0e-8; 496 Vec v_field_l,denom_l,coor_l,denom; 497 PetscInt k,p,e,npoints; 498 PetscInt *mpfield_cell; 499 PetscReal *mpfield_coor; 500 PetscReal xi_p[2]; 501 PetscScalar Ni[3]; 502 PetscSection coordSection; 503 PetscScalar *elcoor = NULL; 504 505 PetscFunctionBegin; 506 ierr = VecZeroEntries(v_field);CHKERRQ(ierr); 507 508 ierr = DMGetLocalVector(dm,&v_field_l);CHKERRQ(ierr); 509 ierr = DMGetGlobalVector(dm,&denom);CHKERRQ(ierr); 510 ierr = DMGetLocalVector(dm,&denom_l);CHKERRQ(ierr); 511 ierr = VecZeroEntries(v_field_l);CHKERRQ(ierr); 512 ierr = VecZeroEntries(denom);CHKERRQ(ierr); 513 ierr = VecZeroEntries(denom_l);CHKERRQ(ierr); 514 515 ierr = DMGetCoordinatesLocal(dm,&coor_l);CHKERRQ(ierr); 516 ierr = DMGetCoordinateSection(dm,&coordSection);CHKERRQ(ierr); 517 518 ierr = DMSwarmGetLocalSize(swarm,&npoints);CHKERRQ(ierr); 519 ierr = DMSwarmGetField(swarm,DMSwarmPICField_coor,NULL,NULL,(void**)&mpfield_coor);CHKERRQ(ierr); 520 ierr = DMSwarmGetField(swarm,DMSwarmPICField_cellid,NULL,NULL,(void**)&mpfield_cell);CHKERRQ(ierr); 521 522 for (p=0; p<npoints; p++) { 523 PetscReal *coor_p,dJ; 524 PetscScalar elfield[3]; 525 PetscBool point_located; 526 527 e = mpfield_cell[p]; 528 coor_p = &mpfield_coor[2*p]; 529 530 ierr = DMPlexVecGetClosure(dm,coordSection,coor_l,e,NULL,&elcoor);CHKERRQ(ierr); 531 532 /* 533 while (!point_located && (failed_counter < 25)) { 534 ierr = PointInTriangle(point, coords[0], coords[1], coords[2], &point_located);CHKERRQ(ierr); 535 point.x = coor_p[0]; 536 point.y = coor_p[1]; 537 point.x += 1.0e-10 * (2.0 * rand()/((double)RAND_MAX)-1.0); 538 point.y += 1.0e-10 * (2.0 * rand()/((double)RAND_MAX)-1.0); 539 failed_counter++; 540 } 541 542 if (!point_located){ 543 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); 544 } 545 546 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); 547 else { 548 ierr = _ComputeLocalCoordinateAffine2d(coor_p,elcoor,xi_p,&dJ);CHKERRQ(ierr); 549 xi_p[0] = 0.5*(xi_p[0] + 1.0); 550 xi_p[1] = 0.5*(xi_p[1] + 1.0); 551 552 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]); 553 554 } 555 */ 556 557 ierr = ComputeLocalCoordinateAffine2d(coor_p,elcoor,xi_p,&dJ);CHKERRQ(ierr); 558 /* 559 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]); 560 */ 561 /* 562 point_located = PETSC_TRUE; 563 if (xi_p[0] < 0.0) { 564 if (xi_p[0] > -PLEX_C_EPS) { 565 xi_p[0] = 0.0; 566 } else { 567 point_located = PETSC_FALSE; 568 } 569 } 570 if (xi_p[1] < 0.0) { 571 if (xi_p[1] > -PLEX_C_EPS) { 572 xi_p[1] = 0.0; 573 } else { 574 point_located = PETSC_FALSE; 575 } 576 } 577 if (xi_p[1] > (1.0-xi_p[0])) { 578 if ((xi_p[1] - 1.0 + xi_p[0]) < PLEX_C_EPS) { 579 xi_p[1] = 1.0 - xi_p[0]; 580 } else { 581 point_located = PETSC_FALSE; 582 } 583 } 584 if (!point_located){ 585 PetscPrintf(PETSC_COMM_SELF,"[Error] xi,eta = %+1.8e, %+1.8e\n",xi_p[0],xi_p[1]); 586 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]); 587 } 588 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); 589 */ 590 591 Ni[0] = 1.0 - xi_p[0] - xi_p[1]; 592 Ni[1] = xi_p[0]; 593 Ni[2] = xi_p[1]; 594 595 point_located = PETSC_TRUE; 596 for (k=0; k<3; k++) { 597 if (PetscRealPart(Ni[k]) < -PLEX_C_EPS) point_located = PETSC_FALSE; 598 if (PetscRealPart(Ni[k]) > (1.0+PLEX_C_EPS)) point_located = PETSC_FALSE; 599 } 600 if (!point_located){ 601 PetscPrintf(PETSC_COMM_SELF,"[Error] xi,eta = %+1.8e, %+1.8e\n",(double)xi_p[0],(double)xi_p[1]); 602 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])); 603 } 604 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); 605 606 for (k=0; k<3; k++) { 607 Ni[k] = Ni[k] * dJ; 608 elfield[k] = Ni[k] * swarm_field[p]; 609 } 610 611 ierr = DMPlexVecRestoreClosure(dm,coordSection,coor_l,e,NULL,&elcoor);CHKERRQ(ierr); 612 613 ierr = DMPlexVecSetClosure(dm, NULL,v_field_l, e, elfield, ADD_VALUES);CHKERRQ(ierr); 614 ierr = DMPlexVecSetClosure(dm, NULL,denom_l, e, Ni, ADD_VALUES);CHKERRQ(ierr); 615 } 616 617 ierr = DMSwarmRestoreField(swarm,DMSwarmPICField_cellid,NULL,NULL,(void**)&mpfield_cell);CHKERRQ(ierr); 618 ierr = DMSwarmRestoreField(swarm,DMSwarmPICField_coor,NULL,NULL,(void**)&mpfield_coor);CHKERRQ(ierr); 619 620 ierr = DMLocalToGlobalBegin(dm,v_field_l,ADD_VALUES,v_field);CHKERRQ(ierr); 621 ierr = DMLocalToGlobalEnd(dm,v_field_l,ADD_VALUES,v_field);CHKERRQ(ierr); 622 ierr = DMLocalToGlobalBegin(dm,denom_l,ADD_VALUES,denom);CHKERRQ(ierr); 623 ierr = DMLocalToGlobalEnd(dm,denom_l,ADD_VALUES,denom);CHKERRQ(ierr); 624 625 ierr = VecPointwiseDivide(v_field,v_field,denom);CHKERRQ(ierr); 626 627 ierr = DMRestoreLocalVector(dm,&v_field_l);CHKERRQ(ierr); 628 ierr = DMRestoreLocalVector(dm,&denom_l);CHKERRQ(ierr); 629 ierr = DMRestoreGlobalVector(dm,&denom);CHKERRQ(ierr); 630 631 PetscFunctionReturn(0); 632 } 633 634 PetscErrorCode private_DMSwarmProjectFields_PLEX(DM swarm,DM celldm,PetscInt project_type,PetscInt nfields,DMSwarmDataField dfield[],Vec vecs[]) 635 { 636 PetscErrorCode ierr; 637 PetscInt f,dim; 638 639 PetscFunctionBegin; 640 ierr = DMGetDimension(swarm,&dim);CHKERRQ(ierr); 641 switch (dim) { 642 case 2: 643 for (f=0; f<nfields; f++) { 644 PetscReal *swarm_field; 645 646 ierr = DMSwarmDataFieldGetEntries(dfield[f],(void**)&swarm_field);CHKERRQ(ierr); 647 ierr = DMSwarmProjectField_ApproxP1_PLEX_2D(swarm,swarm_field,celldm,vecs[f]);CHKERRQ(ierr); 648 } 649 break; 650 case 3: 651 SETERRQ(PetscObjectComm((PetscObject)swarm),PETSC_ERR_SUP,"No support for 3D"); 652 break; 653 default: 654 break; 655 } 656 657 PetscFunctionReturn(0); 658 } 659 660 PetscErrorCode private_DMSwarmSetPointCoordinatesCellwise_PLEX(DM dm,DM dmc,PetscInt npoints,PetscReal xi[]) 661 { 662 PetscBool is_simplex,is_tensorcell; 663 PetscErrorCode ierr; 664 PetscInt dim,nfaces,ps,pe,p,d,nbasis,pcnt,e,k,nel; 665 PetscFE fe; 666 PetscQuadrature quadrature; 667 PetscReal *B,*xiq; 668 Vec coorlocal; 669 PetscSection coordSection; 670 PetscScalar *elcoor = NULL; 671 PetscReal *swarm_coor; 672 PetscInt *swarm_cellid; 673 674 PetscFunctionBegin; 675 ierr = DMGetDimension(dmc,&dim);CHKERRQ(ierr); 676 677 is_simplex = PETSC_FALSE; 678 is_tensorcell = PETSC_FALSE; 679 ierr = DMPlexGetHeightStratum(dmc,0,&ps,&pe);CHKERRQ(ierr); 680 ierr = DMPlexGetConeSize(dmc, ps, &nfaces);CHKERRQ(ierr); 681 682 if (nfaces == (dim+1)) { is_simplex = PETSC_TRUE; } 683 684 switch (dim) { 685 case 2: 686 if (nfaces == 4) { is_tensorcell = PETSC_TRUE; } 687 break; 688 case 3: 689 if (nfaces == 6) { is_tensorcell = PETSC_TRUE; } 690 break; 691 default: 692 SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only support for 2D, 3D"); 693 break; 694 } 695 696 /* check points provided fail inside the reference cell */ 697 if (is_simplex) { 698 for (p=0; p<npoints; p++) { 699 PetscReal sum; 700 for (d=0; d<dim; d++) { 701 if (xi[dim*p+d] < -1.0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Points do not fail inside the simplex domain"); 702 } 703 sum = 0.0; 704 for (d=0; d<dim; d++) { 705 sum += xi[dim*p+d]; 706 } 707 if (sum > 0.0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Points do not fail inside the simplex domain"); 708 } 709 } else if (is_tensorcell) { 710 for (p=0; p<npoints; p++) { 711 for (d=0; d<dim; d++) { 712 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"); 713 } 714 } 715 } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only support for d-simplex and d-tensorcell"); 716 717 ierr = PetscQuadratureCreate(PetscObjectComm((PetscObject)dm),&quadrature);CHKERRQ(ierr); 718 ierr = PetscMalloc1(npoints*dim,&xiq);CHKERRQ(ierr); 719 ierr = PetscMemcpy(xiq,xi,npoints*dim*sizeof(PetscReal));CHKERRQ(ierr); 720 ierr = PetscQuadratureSetData(quadrature,dim,1,npoints,(const PetscReal*)xiq,NULL);CHKERRQ(ierr); 721 ierr = private_PetscFECreateDefault_scalar_pk1(dmc, dim, is_simplex, 0, &fe);CHKERRQ(ierr); 722 ierr = PetscFESetQuadrature(fe,quadrature);CHKERRQ(ierr); 723 ierr = PetscFEGetDimension(fe,&nbasis);CHKERRQ(ierr); 724 ierr = PetscFEGetDefaultTabulation(fe, &B, NULL, NULL);CHKERRQ(ierr); 725 726 /* for each cell, interpolate coordaintes and insert the interpolated points coordinates into swarm */ 727 /* 0->cell, 1->edge, 2->vert */ 728 ierr = DMPlexGetHeightStratum(dmc,0,&ps,&pe);CHKERRQ(ierr); 729 nel = pe - ps; 730 731 ierr = DMSwarmSetLocalSizes(dm,npoints*nel,-1);CHKERRQ(ierr); 732 ierr = DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr); 733 ierr = DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr); 734 735 ierr = DMGetCoordinatesLocal(dmc,&coorlocal);CHKERRQ(ierr); 736 ierr = DMGetCoordinateSection(dmc,&coordSection);CHKERRQ(ierr); 737 738 pcnt = 0; 739 for (e=0; e<nel; e++) { 740 ierr = DMPlexVecGetClosure(dmc,coordSection,coorlocal,ps+e,NULL,&elcoor);CHKERRQ(ierr); 741 742 for (p=0; p<npoints; p++) { 743 for (d=0; d<dim; d++) { 744 swarm_coor[dim*pcnt+d] = 0.0; 745 for (k=0; k<nbasis; k++) { 746 swarm_coor[dim*pcnt+d] += B[p*nbasis + k] * PetscRealPart(elcoor[dim*k+d]); 747 } 748 } 749 swarm_cellid[pcnt] = e; 750 pcnt++; 751 } 752 ierr = DMPlexVecRestoreClosure(dmc,coordSection,coorlocal,ps+e,NULL,&elcoor);CHKERRQ(ierr); 753 } 754 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr); 755 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr); 756 757 ierr = PetscQuadratureDestroy(&quadrature);CHKERRQ(ierr); 758 ierr = PetscFEDestroy(&fe);CHKERRQ(ierr); 759 760 PetscFunctionReturn(0); 761 } 762