1 #include <petscdm.h> 2 #include <petscdmplex.h> 3 #include <petscdmswarm.h> 4 #include "../src/dm/impls/swarm/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 = PetscSpaceSetDegree(P,1,PETSC_DETERMINE);CHKERRQ(ierr); 28 ierr = PetscSpaceSetNumComponents(P, Nc);CHKERRQ(ierr); 29 ierr = PetscSpaceSetNumVariables(P, dim);CHKERRQ(ierr); 30 ierr = PetscSpaceSetUp(P);CHKERRQ(ierr); 31 ierr = PetscSpaceGetDegree(P, &order, NULL);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 = PetscDTStroudConicalQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);CHKERRQ(ierr); 62 ierr = PetscDTStroudConicalQuadrature(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 PetscTabulation T; 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 = PetscFECopyQuadrature(feRef,fe);CHKERRQ(ierr); 216 ierr = PetscFEDestroy(&feRef);CHKERRQ(ierr); 217 } 218 219 ierr = PetscFEGetQuadrature(fe,&quadrature);CHKERRQ(ierr); 220 ierr = PetscQuadratureGetData(quadrature, NULL, NULL, &npoints_q, &xiq, NULL);CHKERRQ(ierr); 221 ierr = PetscFEGetDimension(fe,&nbasis);CHKERRQ(ierr); 222 ierr = PetscFEGetCellTabulation(fe, &T);CHKERRQ(ierr); 223 224 /* 0->cell, 1->edge, 2->vert */ 225 ierr = DMPlexGetHeightStratum(dmc,0,&ps,&pe);CHKERRQ(ierr); 226 nel = pe - ps; 227 228 ierr = DMSwarmSetLocalSizes(dm,npoints_q*nel,-1);CHKERRQ(ierr); 229 ierr = DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr); 230 ierr = DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr); 231 232 ierr = DMGetCoordinatesLocal(dmc,&coorlocal);CHKERRQ(ierr); 233 ierr = DMGetCoordinateSection(dmc,&coordSection);CHKERRQ(ierr); 234 235 pcnt = 0; 236 for (e=0; e<nel; e++) { 237 ierr = DMPlexVecGetClosure(dmc,coordSection,coorlocal,ps+e,NULL,&elcoor);CHKERRQ(ierr); 238 239 for (q=0; q<npoints_q; q++) { 240 for (d=0; d<dim; d++) { 241 swarm_coor[dim*pcnt+d] = 0.0; 242 for (k=0; k<nbasis; k++) { 243 swarm_coor[dim*pcnt+d] += T->T[0][q*nbasis + k] * PetscRealPart(elcoor[dim*k+d]); 244 } 245 } 246 swarm_cellid[pcnt] = e; 247 pcnt++; 248 } 249 ierr = DMPlexVecRestoreClosure(dmc,coordSection,coorlocal,ps+e,NULL,&elcoor);CHKERRQ(ierr); 250 } 251 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr); 252 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr); 253 254 ierr = PetscFEDestroy(&fe);CHKERRQ(ierr); 255 PetscFunctionReturn(0); 256 } 257 258 PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX2D_Regular(DM dm,DM dmc,PetscInt npoints) 259 { 260 PetscErrorCode ierr; 261 PetscInt dim; 262 PetscInt ii,jj,q,npoints_q,e,nel,npe,pcnt,ps,pe,d,k,nfaces; 263 PetscReal *xi,ds,ds2; 264 PetscReal **basis; 265 Vec coorlocal; 266 PetscSection coordSection; 267 PetscScalar *elcoor = NULL; 268 PetscReal *swarm_coor; 269 PetscInt *swarm_cellid; 270 PetscBool is_simplex; 271 272 PetscFunctionBegin; 273 ierr = DMGetDimension(dmc,&dim);CHKERRQ(ierr); 274 if (dim != 2) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only 2D is supported"); 275 is_simplex = PETSC_FALSE; 276 ierr = DMPlexGetHeightStratum(dmc,0,&ps,&pe);CHKERRQ(ierr); 277 ierr = DMPlexGetConeSize(dmc, ps, &nfaces);CHKERRQ(ierr); 278 if (nfaces == (dim+1)) { is_simplex = PETSC_TRUE; } 279 if (!is_simplex) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only the simplex is supported"); 280 281 ierr = PetscMalloc1(dim*npoints*npoints,&xi);CHKERRQ(ierr); 282 pcnt = 0; 283 ds = 1.0/((PetscReal)(npoints-1)); 284 ds2 = 1.0/((PetscReal)(npoints)); 285 for (jj = 0; jj<npoints; jj++) { 286 for (ii=0; ii<npoints-jj; ii++) { 287 xi[dim*pcnt+0] = ii * ds; 288 xi[dim*pcnt+1] = jj * ds; 289 290 xi[dim*pcnt+0] *= (1.0 - 1.2*ds2); 291 xi[dim*pcnt+1] *= (1.0 - 1.2*ds2); 292 293 xi[dim*pcnt+0] += 0.35*ds2; 294 xi[dim*pcnt+1] += 0.35*ds2; 295 296 pcnt++; 297 } 298 } 299 npoints_q = pcnt; 300 301 npe = 3; /* nodes per element (triangle) */ 302 ierr = PetscMalloc1(npoints_q,&basis);CHKERRQ(ierr); 303 for (q=0; q<npoints_q; q++) { 304 ierr = PetscMalloc1(npe,&basis[q]);CHKERRQ(ierr); 305 306 basis[q][0] = 1.0 - xi[dim*q+0] - xi[dim*q+1]; 307 basis[q][1] = xi[dim*q+0]; 308 basis[q][2] = xi[dim*q+1]; 309 } 310 311 /* 0->cell, 1->edge, 2->vert */ 312 ierr = DMPlexGetHeightStratum(dmc,0,&ps,&pe);CHKERRQ(ierr); 313 nel = pe - ps; 314 315 ierr = DMSwarmSetLocalSizes(dm,npoints_q*nel,-1);CHKERRQ(ierr); 316 ierr = DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr); 317 ierr = DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr); 318 319 ierr = DMGetCoordinatesLocal(dmc,&coorlocal);CHKERRQ(ierr); 320 ierr = DMGetCoordinateSection(dmc,&coordSection);CHKERRQ(ierr); 321 322 pcnt = 0; 323 for (e=0; e<nel; e++) { 324 ierr = DMPlexVecGetClosure(dmc,coordSection,coorlocal,e,NULL,&elcoor);CHKERRQ(ierr); 325 326 for (q=0; q<npoints_q; q++) { 327 for (d=0; d<dim; d++) { 328 swarm_coor[dim*pcnt+d] = 0.0; 329 for (k=0; k<npe; k++) { 330 swarm_coor[dim*pcnt+d] += basis[q][k] * PetscRealPart(elcoor[dim*k+d]); 331 } 332 } 333 swarm_cellid[pcnt] = e; 334 pcnt++; 335 } 336 ierr = DMPlexVecRestoreClosure(dmc,coordSection,coorlocal,e,NULL,&elcoor);CHKERRQ(ierr); 337 } 338 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr); 339 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr); 340 341 ierr = PetscFree(xi);CHKERRQ(ierr); 342 for (q=0; q<npoints_q; q++) { 343 ierr = PetscFree(basis[q]);CHKERRQ(ierr); 344 } 345 ierr = PetscFree(basis);CHKERRQ(ierr); 346 347 PetscFunctionReturn(0); 348 } 349 350 PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX(DM dm,DM celldm,DMSwarmPICLayoutType layout,PetscInt layout_param) 351 { 352 PetscErrorCode ierr; 353 PetscInt dim; 354 355 PetscFunctionBegin; 356 ierr = DMGetDimension(celldm,&dim);CHKERRQ(ierr); 357 switch (layout) { 358 case DMSWARMPIC_LAYOUT_REGULAR: 359 if (dim == 3) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No 3D support for REGULAR+PLEX"); 360 ierr = private_DMSwarmInsertPointsUsingCellDM_PLEX2D_Regular(dm,celldm,layout_param);CHKERRQ(ierr); 361 break; 362 case DMSWARMPIC_LAYOUT_GAUSS: 363 { 364 PetscInt npoints,npoints1,ps,pe,nfaces; 365 const PetscReal *xi; 366 PetscBool is_simplex; 367 PetscQuadrature quadrature; 368 369 is_simplex = PETSC_FALSE; 370 ierr = DMPlexGetHeightStratum(celldm,0,&ps,&pe);CHKERRQ(ierr); 371 ierr = DMPlexGetConeSize(celldm, ps, &nfaces);CHKERRQ(ierr); 372 if (nfaces == (dim+1)) { is_simplex = PETSC_TRUE; } 373 374 npoints1 = layout_param; 375 if (is_simplex) { 376 ierr = PetscDTStroudConicalQuadrature(dim,1,npoints1,-1.0,1.0,&quadrature);CHKERRQ(ierr); 377 } else { 378 ierr = PetscDTGaussTensorQuadrature(dim,1,npoints1,-1.0,1.0,&quadrature);CHKERRQ(ierr); 379 } 380 ierr = PetscQuadratureGetData(quadrature,NULL,NULL,&npoints,&xi,NULL);CHKERRQ(ierr); 381 ierr = private_DMSwarmSetPointCoordinatesCellwise_PLEX(dm,celldm,npoints,(PetscReal*)xi);CHKERRQ(ierr); 382 ierr = PetscQuadratureDestroy(&quadrature);CHKERRQ(ierr); 383 } 384 break; 385 case DMSWARMPIC_LAYOUT_SUBDIVISION: 386 ierr = private_DMSwarmInsertPointsUsingCellDM_PLEX_SubDivide(dm,celldm,layout_param);CHKERRQ(ierr); 387 break; 388 } 389 PetscFunctionReturn(0); 390 } 391 392 /* 393 typedef struct { 394 PetscReal x,y; 395 } Point2d; 396 397 static PetscErrorCode signp2d(Point2d p1, Point2d p2, Point2d p3,PetscReal *s) 398 { 399 PetscFunctionBegin; 400 *s = (p1.x - p3.x) * (p2.y - p3.y) - (p2.x - p3.x) * (p1.y - p3.y); 401 PetscFunctionReturn(0); 402 } 403 */ 404 /* 405 static PetscErrorCode PointInTriangle(Point2d pt, Point2d v1, Point2d v2, Point2d v3,PetscBool *v) 406 { 407 PetscReal s1,s2,s3; 408 PetscBool b1, b2, b3; 409 410 PetscFunctionBegin; 411 signp2d(pt, v1, v2,&s1); b1 = s1 < 0.0f; 412 signp2d(pt, v2, v3,&s2); b2 = s2 < 0.0f; 413 signp2d(pt, v3, v1,&s3); b3 = s3 < 0.0f; 414 415 *v = ((b1 == b2) && (b2 == b3)); 416 PetscFunctionReturn(0); 417 } 418 */ 419 /* 420 static PetscErrorCode _ComputeLocalCoordinateAffine2d(PetscReal xp[],PetscReal coords[],PetscReal xip[],PetscReal *dJ) 421 { 422 PetscReal x1,y1,x2,y2,x3,y3; 423 PetscReal c,b[2],A[2][2],inv[2][2],detJ,od; 424 425 PetscFunctionBegin; 426 x1 = coords[2*0+0]; 427 x2 = coords[2*1+0]; 428 x3 = coords[2*2+0]; 429 430 y1 = coords[2*0+1]; 431 y2 = coords[2*1+1]; 432 y3 = coords[2*2+1]; 433 434 c = x1 - 0.5*x1 - 0.5*x1 + 0.5*x2 + 0.5*x3; 435 b[0] = xp[0] - c; 436 c = y1 - 0.5*y1 - 0.5*y1 + 0.5*y2 + 0.5*y3; 437 b[1] = xp[1] - c; 438 439 A[0][0] = -0.5*x1 + 0.5*x2; A[0][1] = -0.5*x1 + 0.5*x3; 440 A[1][0] = -0.5*y1 + 0.5*y2; A[1][1] = -0.5*y1 + 0.5*y3; 441 442 detJ = A[0][0]*A[1][1] - A[0][1]*A[1][0]; 443 *dJ = PetscAbsReal(detJ); 444 od = 1.0/detJ; 445 446 inv[0][0] = A[1][1] * od; 447 inv[0][1] = -A[0][1] * od; 448 inv[1][0] = -A[1][0] * od; 449 inv[1][1] = A[0][0] * od; 450 451 xip[0] = inv[0][0]*b[0] + inv[0][1]*b[1]; 452 xip[1] = inv[1][0]*b[0] + inv[1][1]*b[1]; 453 PetscFunctionReturn(0); 454 } 455 */ 456 457 static PetscErrorCode ComputeLocalCoordinateAffine2d(PetscReal xp[],PetscScalar coords[],PetscReal xip[],PetscReal *dJ) 458 { 459 PetscReal x1,y1,x2,y2,x3,y3; 460 PetscReal b[2],A[2][2],inv[2][2],detJ,od; 461 462 PetscFunctionBegin; 463 x1 = PetscRealPart(coords[2*0+0]); 464 x2 = PetscRealPart(coords[2*1+0]); 465 x3 = PetscRealPart(coords[2*2+0]); 466 467 y1 = PetscRealPart(coords[2*0+1]); 468 y2 = PetscRealPart(coords[2*1+1]); 469 y3 = PetscRealPart(coords[2*2+1]); 470 471 b[0] = xp[0] - x1; 472 b[1] = xp[1] - y1; 473 474 A[0][0] = x2-x1; A[0][1] = x3-x1; 475 A[1][0] = y2-y1; A[1][1] = y3-y1; 476 477 detJ = A[0][0]*A[1][1] - A[0][1]*A[1][0]; 478 *dJ = PetscAbsReal(detJ); 479 od = 1.0/detJ; 480 481 inv[0][0] = A[1][1] * od; 482 inv[0][1] = -A[0][1] * od; 483 inv[1][0] = -A[1][0] * od; 484 inv[1][1] = A[0][0] * od; 485 486 xip[0] = inv[0][0]*b[0] + inv[0][1]*b[1]; 487 xip[1] = inv[1][0]*b[0] + inv[1][1]*b[1]; 488 PetscFunctionReturn(0); 489 } 490 491 PetscErrorCode DMSwarmProjectField_ApproxP1_PLEX_2D(DM swarm,PetscReal *swarm_field,DM dm,Vec v_field) 492 { 493 PetscErrorCode ierr; 494 const PetscReal PLEX_C_EPS = 1.0e-8; 495 Vec v_field_l,denom_l,coor_l,denom; 496 PetscInt k,p,e,npoints; 497 PetscInt *mpfield_cell; 498 PetscReal *mpfield_coor; 499 PetscReal xi_p[2]; 500 PetscScalar Ni[3]; 501 PetscSection coordSection; 502 PetscScalar *elcoor = NULL; 503 504 PetscFunctionBegin; 505 ierr = VecZeroEntries(v_field);CHKERRQ(ierr); 506 507 ierr = DMGetLocalVector(dm,&v_field_l);CHKERRQ(ierr); 508 ierr = DMGetGlobalVector(dm,&denom);CHKERRQ(ierr); 509 ierr = DMGetLocalVector(dm,&denom_l);CHKERRQ(ierr); 510 ierr = VecZeroEntries(v_field_l);CHKERRQ(ierr); 511 ierr = VecZeroEntries(denom);CHKERRQ(ierr); 512 ierr = VecZeroEntries(denom_l);CHKERRQ(ierr); 513 514 ierr = DMGetCoordinatesLocal(dm,&coor_l);CHKERRQ(ierr); 515 ierr = DMGetCoordinateSection(dm,&coordSection);CHKERRQ(ierr); 516 517 ierr = DMSwarmGetLocalSize(swarm,&npoints);CHKERRQ(ierr); 518 ierr = DMSwarmGetField(swarm,DMSwarmPICField_coor,NULL,NULL,(void**)&mpfield_coor);CHKERRQ(ierr); 519 ierr = DMSwarmGetField(swarm,DMSwarmPICField_cellid,NULL,NULL,(void**)&mpfield_cell);CHKERRQ(ierr); 520 521 for (p=0; p<npoints; p++) { 522 PetscReal *coor_p,dJ; 523 PetscScalar elfield[3]; 524 PetscBool point_located; 525 526 e = mpfield_cell[p]; 527 coor_p = &mpfield_coor[2*p]; 528 529 ierr = DMPlexVecGetClosure(dm,coordSection,coor_l,e,NULL,&elcoor);CHKERRQ(ierr); 530 531 /* 532 while (!point_located && (failed_counter < 25)) { 533 ierr = PointInTriangle(point, coords[0], coords[1], coords[2], &point_located);CHKERRQ(ierr); 534 point.x = coor_p[0]; 535 point.y = coor_p[1]; 536 point.x += 1.0e-10 * (2.0 * rand()/((double)RAND_MAX)-1.0); 537 point.y += 1.0e-10 * (2.0 * rand()/((double)RAND_MAX)-1.0); 538 failed_counter++; 539 } 540 541 if (!point_located){ 542 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); 543 } 544 545 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); 546 else { 547 ierr = _ComputeLocalCoordinateAffine2d(coor_p,elcoor,xi_p,&dJ);CHKERRQ(ierr); 548 xi_p[0] = 0.5*(xi_p[0] + 1.0); 549 xi_p[1] = 0.5*(xi_p[1] + 1.0); 550 551 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]); 552 553 } 554 */ 555 556 ierr = ComputeLocalCoordinateAffine2d(coor_p,elcoor,xi_p,&dJ);CHKERRQ(ierr); 557 /* 558 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]); 559 */ 560 /* 561 point_located = PETSC_TRUE; 562 if (xi_p[0] < 0.0) { 563 if (xi_p[0] > -PLEX_C_EPS) { 564 xi_p[0] = 0.0; 565 } else { 566 point_located = PETSC_FALSE; 567 } 568 } 569 if (xi_p[1] < 0.0) { 570 if (xi_p[1] > -PLEX_C_EPS) { 571 xi_p[1] = 0.0; 572 } else { 573 point_located = PETSC_FALSE; 574 } 575 } 576 if (xi_p[1] > (1.0-xi_p[0])) { 577 if ((xi_p[1] - 1.0 + xi_p[0]) < PLEX_C_EPS) { 578 xi_p[1] = 1.0 - xi_p[0]; 579 } else { 580 point_located = PETSC_FALSE; 581 } 582 } 583 if (!point_located){ 584 PetscPrintf(PETSC_COMM_SELF,"[Error] xi,eta = %+1.8e, %+1.8e\n",xi_p[0],xi_p[1]); 585 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]); 586 } 587 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); 588 */ 589 590 Ni[0] = 1.0 - xi_p[0] - xi_p[1]; 591 Ni[1] = xi_p[0]; 592 Ni[2] = xi_p[1]; 593 594 point_located = PETSC_TRUE; 595 for (k=0; k<3; k++) { 596 if (PetscRealPart(Ni[k]) < -PLEX_C_EPS) point_located = PETSC_FALSE; 597 if (PetscRealPart(Ni[k]) > (1.0+PLEX_C_EPS)) point_located = PETSC_FALSE; 598 } 599 if (!point_located){ 600 PetscPrintf(PETSC_COMM_SELF,"[Error] xi,eta = %+1.8e, %+1.8e\n",(double)xi_p[0],(double)xi_p[1]); 601 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])); 602 } 603 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); 604 605 for (k=0; k<3; k++) { 606 Ni[k] = Ni[k] * dJ; 607 elfield[k] = Ni[k] * swarm_field[p]; 608 } 609 610 ierr = DMPlexVecRestoreClosure(dm,coordSection,coor_l,e,NULL,&elcoor);CHKERRQ(ierr); 611 612 ierr = DMPlexVecSetClosure(dm, NULL,v_field_l, e, elfield, ADD_VALUES);CHKERRQ(ierr); 613 ierr = DMPlexVecSetClosure(dm, NULL,denom_l, e, Ni, ADD_VALUES);CHKERRQ(ierr); 614 } 615 616 ierr = DMSwarmRestoreField(swarm,DMSwarmPICField_cellid,NULL,NULL,(void**)&mpfield_cell);CHKERRQ(ierr); 617 ierr = DMSwarmRestoreField(swarm,DMSwarmPICField_coor,NULL,NULL,(void**)&mpfield_coor);CHKERRQ(ierr); 618 619 ierr = DMLocalToGlobalBegin(dm,v_field_l,ADD_VALUES,v_field);CHKERRQ(ierr); 620 ierr = DMLocalToGlobalEnd(dm,v_field_l,ADD_VALUES,v_field);CHKERRQ(ierr); 621 ierr = DMLocalToGlobalBegin(dm,denom_l,ADD_VALUES,denom);CHKERRQ(ierr); 622 ierr = DMLocalToGlobalEnd(dm,denom_l,ADD_VALUES,denom);CHKERRQ(ierr); 623 624 ierr = VecPointwiseDivide(v_field,v_field,denom);CHKERRQ(ierr); 625 626 ierr = DMRestoreLocalVector(dm,&v_field_l);CHKERRQ(ierr); 627 ierr = DMRestoreLocalVector(dm,&denom_l);CHKERRQ(ierr); 628 ierr = DMRestoreGlobalVector(dm,&denom);CHKERRQ(ierr); 629 630 PetscFunctionReturn(0); 631 } 632 633 PetscErrorCode private_DMSwarmProjectFields_PLEX(DM swarm,DM celldm,PetscInt project_type,PetscInt nfields,DMSwarmDataField dfield[],Vec vecs[]) 634 { 635 PetscErrorCode ierr; 636 PetscInt f,dim; 637 638 PetscFunctionBegin; 639 ierr = DMGetDimension(swarm,&dim);CHKERRQ(ierr); 640 switch (dim) { 641 case 2: 642 for (f=0; f<nfields; f++) { 643 PetscReal *swarm_field; 644 645 ierr = DMSwarmDataFieldGetEntries(dfield[f],(void**)&swarm_field);CHKERRQ(ierr); 646 ierr = DMSwarmProjectField_ApproxP1_PLEX_2D(swarm,swarm_field,celldm,vecs[f]);CHKERRQ(ierr); 647 } 648 break; 649 case 3: 650 SETERRQ(PetscObjectComm((PetscObject)swarm),PETSC_ERR_SUP,"No support for 3D"); 651 break; 652 default: 653 break; 654 } 655 656 PetscFunctionReturn(0); 657 } 658 659 PetscErrorCode private_DMSwarmSetPointCoordinatesCellwise_PLEX(DM dm,DM dmc,PetscInt npoints,PetscReal xi[]) 660 { 661 PetscBool is_simplex,is_tensorcell; 662 PetscErrorCode ierr; 663 PetscInt dim,nfaces,ps,pe,p,d,nbasis,pcnt,e,k,nel; 664 PetscFE fe; 665 PetscQuadrature quadrature; 666 PetscTabulation T; 667 PetscReal *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 = PetscArraycpy(xiq,xi,npoints*dim);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 = PetscFEGetCellTabulation(fe, &T);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] += T->T[0][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