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 PetscErrorCode ierr; 18 19 PetscFunctionBegin; 20 /* Create space */ 21 ierr = PetscSpaceCreate(PetscObjectComm((PetscObject) dm), &P);CHKERRQ(ierr); 22 /*ierr = PetscObjectSetOptionsPrefix((PetscObject) P, prefix);CHKERRQ(ierr);*/ 23 ierr = PetscSpacePolynomialSetTensor(P, tensor);CHKERRQ(ierr); 24 /*ierr = PetscSpaceSetFromOptions(P);CHKERRQ(ierr);*/ 25 ierr = PetscSpaceSetType(P,PETSCSPACEPOLYNOMIAL);CHKERRQ(ierr); 26 ierr = PetscSpaceSetDegree(P,1,PETSC_DETERMINE);CHKERRQ(ierr); 27 ierr = PetscSpaceSetNumComponents(P, Nc);CHKERRQ(ierr); 28 ierr = PetscSpaceSetNumVariables(P, dim);CHKERRQ(ierr); 29 ierr = PetscSpaceSetUp(P);CHKERRQ(ierr); 30 ierr = PetscSpaceGetDegree(P, &order, NULL);CHKERRQ(ierr); 31 ierr = PetscSpacePolynomialGetTensor(P, &tensor);CHKERRQ(ierr); 32 /* Create dual space */ 33 ierr = PetscDualSpaceCreate(PetscObjectComm((PetscObject) dm), &Q);CHKERRQ(ierr); 34 ierr = PetscDualSpaceSetType(Q,PETSCDUALSPACELAGRANGE);CHKERRQ(ierr); 35 /*ierr = PetscObjectSetOptionsPrefix((PetscObject) Q, prefix);CHKERRQ(ierr);*/ 36 ierr = DMPlexCreateReferenceCell(PETSC_COMM_SELF, DMPolytopeTypeSimpleShape(dim, isSimplex), &K);CHKERRQ(ierr); 37 ierr = PetscDualSpaceSetDM(Q, K);CHKERRQ(ierr); 38 ierr = DMDestroy(&K);CHKERRQ(ierr); 39 ierr = PetscDualSpaceSetNumComponents(Q, Nc);CHKERRQ(ierr); 40 ierr = PetscDualSpaceSetOrder(Q, order);CHKERRQ(ierr); 41 ierr = PetscDualSpaceLagrangeSetTensor(Q, tensor);CHKERRQ(ierr); 42 /*ierr = PetscDualSpaceSetFromOptions(Q);CHKERRQ(ierr);*/ 43 ierr = PetscDualSpaceSetType(Q,PETSCDUALSPACELAGRANGE);CHKERRQ(ierr); 44 ierr = PetscDualSpaceSetUp(Q);CHKERRQ(ierr); 45 /* Create element */ 46 ierr = PetscFECreate(PetscObjectComm((PetscObject) dm), fem);CHKERRQ(ierr); 47 /*ierr = PetscObjectSetOptionsPrefix((PetscObject) *fem, prefix);CHKERRQ(ierr);*/ 48 /*ierr = PetscFESetFromOptions(*fem);CHKERRQ(ierr);*/ 49 ierr = PetscFESetType(*fem,PETSCFEBASIC);CHKERRQ(ierr); 50 ierr = PetscFESetBasisSpace(*fem, P);CHKERRQ(ierr); 51 ierr = PetscFESetDualSpace(*fem, Q);CHKERRQ(ierr); 52 ierr = PetscFESetNumComponents(*fem, Nc);CHKERRQ(ierr); 53 ierr = PetscFESetUp(*fem);CHKERRQ(ierr); 54 ierr = PetscSpaceDestroy(&P);CHKERRQ(ierr); 55 ierr = PetscDualSpaceDestroy(&Q);CHKERRQ(ierr); 56 /* Create quadrature (with specified order if given) */ 57 qorder = qorder >= 0 ? qorder : order; 58 quadPointsPerEdge = PetscMax(qorder + 1,1); 59 if (isSimplex) { 60 ierr = PetscDTStroudConicalQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);CHKERRQ(ierr); 61 ierr = PetscDTStroudConicalQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);CHKERRQ(ierr); 62 } 63 else { 64 ierr = PetscDTGaussTensorQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);CHKERRQ(ierr); 65 ierr = PetscDTGaussTensorQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);CHKERRQ(ierr); 66 } 67 ierr = PetscFESetQuadrature(*fem, q);CHKERRQ(ierr); 68 ierr = PetscFESetFaceQuadrature(*fem, fq);CHKERRQ(ierr); 69 ierr = PetscQuadratureDestroy(&q);CHKERRQ(ierr); 70 ierr = PetscQuadratureDestroy(&fq);CHKERRQ(ierr); 71 PetscFunctionReturn(0); 72 } 73 74 PetscErrorCode subdivide_triangle(PetscReal v1[2],PetscReal v2[2],PetscReal v3[2],PetscInt depth,PetscInt max,PetscReal xi[],PetscInt *np) 75 { 76 PetscReal v12[2],v23[2],v31[2]; 77 PetscInt i; 78 PetscErrorCode ierr; 79 80 PetscFunctionBegin; 81 if (depth == max) { 82 PetscReal cx[2]; 83 84 cx[0] = (v1[0] + v2[0] + v3[0])/3.0; 85 cx[1] = (v1[1] + v2[1] + v3[1])/3.0; 86 87 xi[2*(*np)+0] = cx[0]; 88 xi[2*(*np)+1] = cx[1]; 89 (*np)++; 90 PetscFunctionReturn(0); 91 } 92 93 /* calculate midpoints of each side */ 94 for (i = 0; i < 2; i++) { 95 v12[i] = (v1[i]+v2[i])/2.0; 96 v23[i] = (v2[i]+v3[i])/2.0; 97 v31[i] = (v3[i]+v1[i])/2.0; 98 } 99 100 /* recursively subdivide new triangles */ 101 ierr = subdivide_triangle(v1,v12,v31,depth+1,max,xi,np);CHKERRQ(ierr); 102 ierr = subdivide_triangle(v2,v23,v12,depth+1,max,xi,np);CHKERRQ(ierr); 103 ierr = subdivide_triangle(v3,v31,v23,depth+1,max,xi,np);CHKERRQ(ierr); 104 ierr = subdivide_triangle(v12,v23,v31,depth+1,max,xi,np);CHKERRQ(ierr); 105 PetscFunctionReturn(0); 106 } 107 108 PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX2D_SubDivide(DM dm,DM dmc,PetscInt nsub) 109 { 110 PetscErrorCode ierr; 111 const PetscInt dim = 2; 112 PetscInt q,npoints_q,e,nel,npe,pcnt,ps,pe,d,k,depth; 113 PetscReal *xi; 114 PetscReal **basis; 115 Vec coorlocal; 116 PetscSection coordSection; 117 PetscScalar *elcoor = NULL; 118 PetscReal *swarm_coor; 119 PetscInt *swarm_cellid; 120 PetscReal v1[2],v2[2],v3[2]; 121 122 PetscFunctionBegin; 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 PetscFunctionReturn(0); 180 } 181 182 PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX_SubDivide(DM dm,DM dmc,PetscInt nsub) 183 { 184 PetscErrorCode ierr; 185 PetscInt dim,nfaces,nbasis; 186 PetscInt q,npoints_q,e,nel,pcnt,ps,pe,d,k,r; 187 PetscTabulation T; 188 Vec coorlocal; 189 PetscSection coordSection; 190 PetscScalar *elcoor = NULL; 191 PetscReal *swarm_coor; 192 PetscInt *swarm_cellid; 193 const PetscReal *xiq; 194 PetscQuadrature quadrature; 195 PetscFE fe,feRef; 196 PetscBool is_simplex; 197 198 PetscFunctionBegin; 199 ierr = DMGetDimension(dmc,&dim);CHKERRQ(ierr); 200 is_simplex = PETSC_FALSE; 201 ierr = DMPlexGetHeightStratum(dmc,0,&ps,&pe);CHKERRQ(ierr); 202 ierr = DMPlexGetConeSize(dmc, ps, &nfaces);CHKERRQ(ierr); 203 if (nfaces == (dim+1)) { is_simplex = PETSC_TRUE; } 204 205 ierr = private_PetscFECreateDefault_scalar_pk1(dmc, dim, is_simplex, 0, &fe);CHKERRQ(ierr); 206 207 for (r=0; r<nsub; r++) { 208 ierr = PetscFERefine(fe,&feRef);CHKERRQ(ierr); 209 ierr = PetscFECopyQuadrature(feRef,fe);CHKERRQ(ierr); 210 ierr = PetscFEDestroy(&feRef);CHKERRQ(ierr); 211 } 212 213 ierr = PetscFEGetQuadrature(fe,&quadrature);CHKERRQ(ierr); 214 ierr = PetscQuadratureGetData(quadrature, NULL, NULL, &npoints_q, &xiq, NULL);CHKERRQ(ierr); 215 ierr = PetscFEGetDimension(fe,&nbasis);CHKERRQ(ierr); 216 ierr = PetscFEGetCellTabulation(fe, 1, &T);CHKERRQ(ierr); 217 218 /* 0->cell, 1->edge, 2->vert */ 219 ierr = DMPlexGetHeightStratum(dmc,0,&ps,&pe);CHKERRQ(ierr); 220 nel = pe - ps; 221 222 ierr = DMSwarmSetLocalSizes(dm,npoints_q*nel,-1);CHKERRQ(ierr); 223 ierr = DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr); 224 ierr = DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr); 225 226 ierr = DMGetCoordinatesLocal(dmc,&coorlocal);CHKERRQ(ierr); 227 ierr = DMGetCoordinateSection(dmc,&coordSection);CHKERRQ(ierr); 228 229 pcnt = 0; 230 for (e=0; e<nel; e++) { 231 ierr = DMPlexVecGetClosure(dmc,coordSection,coorlocal,ps+e,NULL,&elcoor);CHKERRQ(ierr); 232 233 for (q=0; q<npoints_q; q++) { 234 for (d=0; d<dim; d++) { 235 swarm_coor[dim*pcnt+d] = 0.0; 236 for (k=0; k<nbasis; k++) { 237 swarm_coor[dim*pcnt+d] += T->T[0][q*nbasis + k] * PetscRealPart(elcoor[dim*k+d]); 238 } 239 } 240 swarm_cellid[pcnt] = e; 241 pcnt++; 242 } 243 ierr = DMPlexVecRestoreClosure(dmc,coordSection,coorlocal,ps+e,NULL,&elcoor);CHKERRQ(ierr); 244 } 245 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr); 246 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr); 247 248 ierr = PetscFEDestroy(&fe);CHKERRQ(ierr); 249 PetscFunctionReturn(0); 250 } 251 252 PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX2D_Regular(DM dm,DM dmc,PetscInt npoints) 253 { 254 PetscErrorCode ierr; 255 PetscInt dim; 256 PetscInt ii,jj,q,npoints_q,e,nel,npe,pcnt,ps,pe,d,k,nfaces; 257 PetscReal *xi,ds,ds2; 258 PetscReal **basis; 259 Vec coorlocal; 260 PetscSection coordSection; 261 PetscScalar *elcoor = NULL; 262 PetscReal *swarm_coor; 263 PetscInt *swarm_cellid; 264 PetscBool is_simplex; 265 266 PetscFunctionBegin; 267 ierr = DMGetDimension(dmc,&dim);CHKERRQ(ierr); 268 PetscCheckFalse(dim != 2,PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only 2D is supported"); 269 is_simplex = PETSC_FALSE; 270 ierr = DMPlexGetHeightStratum(dmc,0,&ps,&pe);CHKERRQ(ierr); 271 ierr = DMPlexGetConeSize(dmc, ps, &nfaces);CHKERRQ(ierr); 272 if (nfaces == (dim+1)) { is_simplex = PETSC_TRUE; } 273 PetscCheckFalse(!is_simplex,PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only the simplex is supported"); 274 275 ierr = PetscMalloc1(dim*npoints*npoints,&xi);CHKERRQ(ierr); 276 pcnt = 0; 277 ds = 1.0/((PetscReal)(npoints-1)); 278 ds2 = 1.0/((PetscReal)(npoints)); 279 for (jj = 0; jj<npoints; jj++) { 280 for (ii=0; ii<npoints-jj; ii++) { 281 xi[dim*pcnt+0] = ii * ds; 282 xi[dim*pcnt+1] = jj * ds; 283 284 xi[dim*pcnt+0] *= (1.0 - 1.2*ds2); 285 xi[dim*pcnt+1] *= (1.0 - 1.2*ds2); 286 287 xi[dim*pcnt+0] += 0.35*ds2; 288 xi[dim*pcnt+1] += 0.35*ds2; 289 pcnt++; 290 } 291 } 292 npoints_q = pcnt; 293 294 npe = 3; /* nodes per element (triangle) */ 295 ierr = PetscMalloc1(npoints_q,&basis);CHKERRQ(ierr); 296 for (q=0; q<npoints_q; q++) { 297 ierr = PetscMalloc1(npe,&basis[q]);CHKERRQ(ierr); 298 299 basis[q][0] = 1.0 - xi[dim*q+0] - xi[dim*q+1]; 300 basis[q][1] = xi[dim*q+0]; 301 basis[q][2] = xi[dim*q+1]; 302 } 303 304 /* 0->cell, 1->edge, 2->vert */ 305 ierr = DMPlexGetHeightStratum(dmc,0,&ps,&pe);CHKERRQ(ierr); 306 nel = pe - ps; 307 308 ierr = DMSwarmSetLocalSizes(dm,npoints_q*nel,-1);CHKERRQ(ierr); 309 ierr = DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr); 310 ierr = DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr); 311 312 ierr = DMGetCoordinatesLocal(dmc,&coorlocal);CHKERRQ(ierr); 313 ierr = DMGetCoordinateSection(dmc,&coordSection);CHKERRQ(ierr); 314 315 pcnt = 0; 316 for (e=0; e<nel; e++) { 317 ierr = DMPlexVecGetClosure(dmc,coordSection,coorlocal,e,NULL,&elcoor);CHKERRQ(ierr); 318 319 for (q=0; q<npoints_q; q++) { 320 for (d=0; d<dim; d++) { 321 swarm_coor[dim*pcnt+d] = 0.0; 322 for (k=0; k<npe; k++) { 323 swarm_coor[dim*pcnt+d] += basis[q][k] * PetscRealPart(elcoor[dim*k+d]); 324 } 325 } 326 swarm_cellid[pcnt] = e; 327 pcnt++; 328 } 329 ierr = DMPlexVecRestoreClosure(dmc,coordSection,coorlocal,e,NULL,&elcoor);CHKERRQ(ierr); 330 } 331 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr); 332 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr); 333 334 ierr = PetscFree(xi);CHKERRQ(ierr); 335 for (q=0; q<npoints_q; q++) { 336 ierr = PetscFree(basis[q]);CHKERRQ(ierr); 337 } 338 ierr = PetscFree(basis);CHKERRQ(ierr); 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 PetscCheckFalse(dim == 3,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 = PetscDTStroudConicalQuadrature(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 PetscCheckFalse(!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); 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 PetscCheckFalse(!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); 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 ierr = PetscPrintf(PETSC_COMM_SELF,"[Error] xi,eta = %+1.8e, %+1.8e\n",(double)xi_p[0],(double)xi_p[1]);CHKERRQ(ierr); 593 ierr = 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]));CHKERRQ(ierr); 594 } 595 PetscCheckFalse(!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); 596 597 for (k=0; k<3; k++) { 598 Ni[k] = Ni[k] * dJ; 599 elfield[k] = Ni[k] * swarm_field[p]; 600 } 601 ierr = DMPlexVecRestoreClosure(dm,coordSection,coor_l,e,NULL,&elcoor);CHKERRQ(ierr); 602 603 ierr = DMPlexVecSetClosure(dm, NULL,v_field_l, e, elfield, ADD_VALUES);CHKERRQ(ierr); 604 ierr = DMPlexVecSetClosure(dm, NULL,denom_l, e, Ni, ADD_VALUES);CHKERRQ(ierr); 605 } 606 607 ierr = DMSwarmRestoreField(swarm,DMSwarmPICField_cellid,NULL,NULL,(void**)&mpfield_cell);CHKERRQ(ierr); 608 ierr = DMSwarmRestoreField(swarm,DMSwarmPICField_coor,NULL,NULL,(void**)&mpfield_coor);CHKERRQ(ierr); 609 610 ierr = DMLocalToGlobalBegin(dm,v_field_l,ADD_VALUES,v_field);CHKERRQ(ierr); 611 ierr = DMLocalToGlobalEnd(dm,v_field_l,ADD_VALUES,v_field);CHKERRQ(ierr); 612 ierr = DMLocalToGlobalBegin(dm,denom_l,ADD_VALUES,denom);CHKERRQ(ierr); 613 ierr = DMLocalToGlobalEnd(dm,denom_l,ADD_VALUES,denom);CHKERRQ(ierr); 614 615 ierr = VecPointwiseDivide(v_field,v_field,denom);CHKERRQ(ierr); 616 617 ierr = DMRestoreLocalVector(dm,&v_field_l);CHKERRQ(ierr); 618 ierr = DMRestoreLocalVector(dm,&denom_l);CHKERRQ(ierr); 619 ierr = DMRestoreGlobalVector(dm,&denom);CHKERRQ(ierr); 620 PetscFunctionReturn(0); 621 } 622 623 PetscErrorCode private_DMSwarmProjectFields_PLEX(DM swarm,DM celldm,PetscInt project_type,PetscInt nfields,DMSwarmDataField dfield[],Vec vecs[]) 624 { 625 PetscErrorCode ierr; 626 PetscInt f,dim; 627 628 PetscFunctionBegin; 629 ierr = DMGetDimension(swarm,&dim);CHKERRQ(ierr); 630 switch (dim) { 631 case 2: 632 for (f=0; f<nfields; f++) { 633 PetscReal *swarm_field; 634 635 ierr = DMSwarmDataFieldGetEntries(dfield[f],(void**)&swarm_field);CHKERRQ(ierr); 636 ierr = DMSwarmProjectField_ApproxP1_PLEX_2D(swarm,swarm_field,celldm,vecs[f]);CHKERRQ(ierr); 637 } 638 break; 639 case 3: 640 SETERRQ(PetscObjectComm((PetscObject)swarm),PETSC_ERR_SUP,"No support for 3D"); 641 default: 642 break; 643 } 644 PetscFunctionReturn(0); 645 } 646 647 PetscErrorCode private_DMSwarmSetPointCoordinatesCellwise_PLEX(DM dm,DM dmc,PetscInt npoints,PetscReal xi[]) 648 { 649 PetscBool is_simplex,is_tensorcell; 650 PetscErrorCode ierr; 651 PetscInt dim,nfaces,ps,pe,p,d,nbasis,pcnt,e,k,nel; 652 PetscFE fe; 653 PetscQuadrature quadrature; 654 PetscTabulation T; 655 PetscReal *xiq; 656 Vec coorlocal; 657 PetscSection coordSection; 658 PetscScalar *elcoor = NULL; 659 PetscReal *swarm_coor; 660 PetscInt *swarm_cellid; 661 662 PetscFunctionBegin; 663 ierr = DMGetDimension(dmc,&dim);CHKERRQ(ierr); 664 665 is_simplex = PETSC_FALSE; 666 is_tensorcell = PETSC_FALSE; 667 ierr = DMPlexGetHeightStratum(dmc,0,&ps,&pe);CHKERRQ(ierr); 668 ierr = DMPlexGetConeSize(dmc, ps, &nfaces);CHKERRQ(ierr); 669 670 if (nfaces == (dim+1)) { is_simplex = PETSC_TRUE; } 671 672 switch (dim) { 673 case 2: 674 if (nfaces == 4) { is_tensorcell = PETSC_TRUE; } 675 break; 676 case 3: 677 if (nfaces == 6) { is_tensorcell = PETSC_TRUE; } 678 break; 679 default: 680 SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only support for 2D, 3D"); 681 } 682 683 /* check points provided fail inside the reference cell */ 684 if (is_simplex) { 685 for (p=0; p<npoints; p++) { 686 PetscReal sum; 687 for (d=0; d<dim; d++) { 688 PetscCheckFalse(xi[dim*p+d] < -1.0,PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Points do not fail inside the simplex domain"); 689 } 690 sum = 0.0; 691 for (d=0; d<dim; d++) { 692 sum += xi[dim*p+d]; 693 } 694 PetscCheckFalse(sum > 0.0,PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Points do not fail inside the simplex domain"); 695 } 696 } else if (is_tensorcell) { 697 for (p=0; p<npoints; p++) { 698 for (d=0; d<dim; d++) { 699 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"); 700 } 701 } 702 } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only support for d-simplex and d-tensorcell"); 703 704 ierr = PetscQuadratureCreate(PetscObjectComm((PetscObject)dm),&quadrature);CHKERRQ(ierr); 705 ierr = PetscMalloc1(npoints*dim,&xiq);CHKERRQ(ierr); 706 ierr = PetscArraycpy(xiq,xi,npoints*dim);CHKERRQ(ierr); 707 ierr = PetscQuadratureSetData(quadrature,dim,1,npoints,(const PetscReal*)xiq,NULL);CHKERRQ(ierr); 708 ierr = private_PetscFECreateDefault_scalar_pk1(dmc, dim, is_simplex, 0, &fe);CHKERRQ(ierr); 709 ierr = PetscFESetQuadrature(fe,quadrature);CHKERRQ(ierr); 710 ierr = PetscFEGetDimension(fe,&nbasis);CHKERRQ(ierr); 711 ierr = PetscFEGetCellTabulation(fe, 1, &T);CHKERRQ(ierr); 712 713 /* for each cell, interpolate coordaintes and insert the interpolated points coordinates into swarm */ 714 /* 0->cell, 1->edge, 2->vert */ 715 ierr = DMPlexGetHeightStratum(dmc,0,&ps,&pe);CHKERRQ(ierr); 716 nel = pe - ps; 717 718 ierr = DMSwarmSetLocalSizes(dm,npoints*nel,-1);CHKERRQ(ierr); 719 ierr = DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr); 720 ierr = DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr); 721 722 ierr = DMGetCoordinatesLocal(dmc,&coorlocal);CHKERRQ(ierr); 723 ierr = DMGetCoordinateSection(dmc,&coordSection);CHKERRQ(ierr); 724 725 pcnt = 0; 726 for (e=0; e<nel; e++) { 727 ierr = DMPlexVecGetClosure(dmc,coordSection,coorlocal,ps+e,NULL,&elcoor);CHKERRQ(ierr); 728 729 for (p=0; p<npoints; p++) { 730 for (d=0; d<dim; d++) { 731 swarm_coor[dim*pcnt+d] = 0.0; 732 for (k=0; k<nbasis; k++) { 733 swarm_coor[dim*pcnt+d] += T->T[0][p*nbasis + k] * PetscRealPart(elcoor[dim*k+d]); 734 } 735 } 736 swarm_cellid[pcnt] = e; 737 pcnt++; 738 } 739 ierr = DMPlexVecRestoreClosure(dmc,coordSection,coorlocal,ps+e,NULL,&elcoor);CHKERRQ(ierr); 740 } 741 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr); 742 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr); 743 744 ierr = PetscQuadratureDestroy(&quadrature);CHKERRQ(ierr); 745 ierr = PetscFEDestroy(&fe);CHKERRQ(ierr); 746 PetscFunctionReturn(0); 747 } 748