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