1 #include <petscdm.h> 2 #include <petscdmplex.h> 3 #include <petscdmswarm.h> 4 #include "data_bucket.h" 5 6 PetscErrorCode subdivide_triangle(PetscReal v1[2],PetscReal v2[2],PetscReal v3[3],PetscInt depth,PetscInt max,PetscReal xi[],PetscInt *np) 7 { 8 PetscReal v12[2],v23[2],v31[2]; 9 PetscInt i; 10 PetscErrorCode ierr; 11 12 PetscFunctionBegin; 13 if (depth == max) { 14 PetscReal cx[2]; 15 16 cx[0] = (v1[0] + v2[0] + v3[0])/3.0; 17 cx[1] = (v1[1] + v2[1] + v3[1])/3.0; 18 19 xi[2*(*np)+0] = cx[0]; 20 xi[2*(*np)+1] = cx[1]; 21 (*np)++; 22 PetscFunctionReturn(0); 23 } 24 25 /* calculate midpoints of each side */ 26 for (i = 0; i < 2; i++) { 27 v12[i] = (v1[i]+v2[i])/2.0; 28 v23[i] = (v2[i]+v3[i])/2.0; 29 v31[i] = (v3[i]+v1[i])/2.0; 30 } 31 32 /* recursively subdivide new triangles */ 33 ierr = subdivide_triangle(v1,v12,v31,depth+1,max,xi,np);CHKERRQ(ierr); 34 ierr = subdivide_triangle(v2,v23,v12,depth+1,max,xi,np);CHKERRQ(ierr); 35 ierr = subdivide_triangle(v3,v31,v23,depth+1,max,xi,np);CHKERRQ(ierr); 36 ierr = subdivide_triangle(v12,v23,v31,depth+1,max,xi,np);CHKERRQ(ierr); 37 PetscFunctionReturn(0); 38 } 39 40 41 PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX2D_SubDivide(DM dm,DM dmc,PetscInt nsub) 42 { 43 PetscErrorCode ierr; 44 const PetscInt dim = 2; 45 PetscInt q,npoints_q,e,nel,npe,pcnt,ps,pe,d,k,depth; 46 PetscReal *xi; 47 PetscReal **basis; 48 Vec coorlocal; 49 PetscSection coordSection; 50 PetscScalar *elcoor = NULL; 51 PetscReal *swarm_coor; 52 PetscInt *swarm_cellid; 53 PetscReal v1[2],v2[2],v3[2]; 54 55 PetscFunctionBegin; 56 57 npoints_q = 1; 58 for (d=0; d<nsub; d++) { npoints_q *= 4; } 59 ierr = PetscMalloc1(dim*npoints_q,&xi);CHKERRQ(ierr); 60 61 v1[0] = 0.0; v1[1] = 0.0; 62 v2[0] = 1.0; v2[1] = 0.0; 63 v3[0] = 0.0; v3[1] = 1.0; 64 depth = 0; 65 pcnt = 0; 66 ierr = subdivide_triangle(v1,v2,v3,depth,nsub,xi,&pcnt);CHKERRQ(ierr); 67 68 npe = 3; /* nodes per element (triangle) */ 69 ierr = PetscMalloc1(npoints_q,&basis);CHKERRQ(ierr); 70 for (q=0; q<npoints_q; q++) { 71 ierr = PetscMalloc1(npe,&basis[q]);CHKERRQ(ierr); 72 73 basis[q][0] = 1.0 - xi[dim*q+0] - xi[dim*q+1]; 74 basis[q][1] = xi[dim*q+0]; 75 basis[q][2] = xi[dim*q+1]; 76 } 77 78 /* 0->cell, 1->edge, 2->vert */ 79 ierr = DMPlexGetHeightStratum(dmc,0,&ps,&pe);CHKERRQ(ierr); 80 nel = pe - ps; 81 82 ierr = DMSwarmSetLocalSizes(dm,npoints_q*nel,-1);CHKERRQ(ierr); 83 ierr = DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr); 84 ierr = DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr); 85 86 ierr = DMGetCoordinatesLocal(dmc,&coorlocal);CHKERRQ(ierr); 87 ierr = DMGetCoordinateSection(dmc,&coordSection);CHKERRQ(ierr); 88 89 pcnt = 0; 90 for (e=0; e<nel; e++) { 91 ierr = DMPlexVecGetClosure(dmc,coordSection,coorlocal,e,NULL,&elcoor);CHKERRQ(ierr); 92 93 for (q=0; q<npoints_q; q++) { 94 for (d=0; d<dim; d++) { 95 swarm_coor[dim*pcnt+d] = 0.0; 96 for (k=0; k<npe; k++) { 97 swarm_coor[dim*pcnt+d] += basis[q][k] * PetscRealPart(elcoor[dim*k+d]); 98 } 99 } 100 swarm_cellid[pcnt] = e; 101 pcnt++; 102 } 103 ierr = DMPlexVecRestoreClosure(dmc,coordSection,coorlocal,e,NULL,&elcoor);CHKERRQ(ierr); 104 } 105 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr); 106 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr); 107 108 ierr = PetscFree(xi);CHKERRQ(ierr); 109 for (q=0; q<npoints_q; q++) { 110 ierr = PetscFree(basis[q]);CHKERRQ(ierr); 111 } 112 ierr = PetscFree(basis);CHKERRQ(ierr); 113 114 PetscFunctionReturn(0); 115 } 116 117 PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX_SubDivide(DM dm,DM dmc,PetscInt nsub) 118 { 119 PetscErrorCode ierr; 120 PetscInt dim,nfaces,nbasis; 121 PetscInt q,npoints_q,e,nel,pcnt,ps,pe,d,k,r; 122 PetscReal *B; 123 Vec coorlocal; 124 PetscSection coordSection; 125 PetscScalar *elcoor = NULL; 126 PetscReal *swarm_coor; 127 PetscInt *swarm_cellid; 128 const PetscReal *xiq; 129 PetscQuadrature quadrature; 130 PetscFE fe,feRef; 131 PetscBool is_simplex; 132 133 PetscFunctionBegin; 134 135 ierr = DMGetDimension(dmc,&dim);CHKERRQ(ierr); 136 137 is_simplex = PETSC_FALSE; 138 ierr = DMPlexGetHeightStratum(dmc,0,&ps,&pe);CHKERRQ(ierr); 139 ierr = DMPlexGetConeSize(dmc, ps, &nfaces);CHKERRQ(ierr); 140 if (nfaces == (dim+1)) { is_simplex = PETSC_TRUE; } 141 142 ierr = PetscFECreateDefault(dmc, dim, 1, is_simplex, NULL, 0, &fe);CHKERRQ(ierr); 143 144 for (r=0; r<nsub; r++) { 145 ierr = PetscFERefine(fe,&feRef);CHKERRQ(ierr); 146 ierr = PetscFEGetQuadrature(feRef,&quadrature);CHKERRQ(ierr); 147 ierr = PetscFESetQuadrature(fe,quadrature);CHKERRQ(ierr); 148 ierr = PetscFEDestroy(&feRef);CHKERRQ(ierr); 149 } 150 151 ierr = PetscFEGetQuadrature(fe,&quadrature);CHKERRQ(ierr); 152 ierr = PetscQuadratureGetData(quadrature, NULL, NULL, &npoints_q, &xiq, NULL);CHKERRQ(ierr); 153 ierr = PetscFEGetDimension(fe,&nbasis);CHKERRQ(ierr); 154 ierr = PetscFEGetDefaultTabulation(fe, &B, NULL, NULL);CHKERRQ(ierr); 155 156 /* 0->cell, 1->edge, 2->vert */ 157 ierr = DMPlexGetHeightStratum(dmc,0,&ps,&pe);CHKERRQ(ierr); 158 nel = pe - ps; 159 160 ierr = DMSwarmSetLocalSizes(dm,npoints_q*nel,-1);CHKERRQ(ierr); 161 ierr = DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr); 162 ierr = DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr); 163 164 ierr = DMGetCoordinatesLocal(dmc,&coorlocal);CHKERRQ(ierr); 165 ierr = DMGetCoordinateSection(dmc,&coordSection);CHKERRQ(ierr); 166 167 pcnt = 0; 168 for (e=0; e<nel; e++) { 169 ierr = DMPlexVecGetClosure(dmc,coordSection,coorlocal,ps+e,NULL,&elcoor);CHKERRQ(ierr); 170 171 for (q=0; q<npoints_q; q++) { 172 for (d=0; d<dim; d++) { 173 swarm_coor[dim*pcnt+d] = 0.0; 174 for (k=0; k<nbasis; k++) { 175 swarm_coor[dim*pcnt+d] += B[q*nbasis + k] * PetscRealPart(elcoor[dim*k+d]); 176 } 177 } 178 swarm_cellid[pcnt] = e; 179 pcnt++; 180 } 181 ierr = DMPlexVecRestoreClosure(dmc,coordSection,coorlocal,e,NULL,&elcoor);CHKERRQ(ierr); 182 } 183 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr); 184 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr); 185 186 ierr = PetscFEDestroy(&fe);CHKERRQ(ierr); 187 PetscFunctionReturn(0); 188 } 189 190 PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX2D_Regular(DM dm,DM dmc,PetscInt npoints) 191 { 192 PetscErrorCode ierr; 193 const PetscInt dim = 2; 194 PetscInt ii,jj,q,npoints_q,e,nel,npe,pcnt,ps,pe,d,k; 195 PetscReal *xi,ds,ds2; 196 PetscReal **basis; 197 Vec coorlocal; 198 PetscSection coordSection; 199 PetscScalar *elcoor = NULL; 200 PetscReal *swarm_coor; 201 PetscInt *swarm_cellid; 202 203 PetscFunctionBegin; 204 ierr = PetscMalloc1(dim*npoints*npoints,&xi);CHKERRQ(ierr); 205 pcnt = 0; 206 ds = 1.0/((PetscReal)(npoints-1)); 207 ds2 = 1.0/((PetscReal)(npoints)); 208 for (jj = 0; jj<npoints; jj++) { 209 for (ii=0; ii<npoints-jj; ii++) { 210 xi[dim*pcnt+0] = ii * ds; 211 xi[dim*pcnt+1] = jj * ds; 212 213 xi[dim*pcnt+0] *= (1.0 - 1.2*ds2); 214 xi[dim*pcnt+1] *= (1.0 - 1.2*ds2); 215 216 xi[dim*pcnt+0] += 0.35*ds2; 217 xi[dim*pcnt+1] += 0.35*ds2; 218 219 pcnt++; 220 } 221 } 222 npoints_q = pcnt; 223 224 npe = 3; /* nodes per element (triangle) */ 225 ierr = PetscMalloc1(npoints_q,&basis);CHKERRQ(ierr); 226 for (q=0; q<npoints_q; q++) { 227 ierr = PetscMalloc1(npe,&basis[q]);CHKERRQ(ierr); 228 229 basis[q][0] = 1.0 - xi[dim*q+0] - xi[dim*q+1]; 230 basis[q][1] = xi[dim*q+0]; 231 basis[q][2] = xi[dim*q+1]; 232 } 233 234 /* 0->cell, 1->edge, 2->vert */ 235 ierr = DMPlexGetHeightStratum(dmc,0,&ps,&pe);CHKERRQ(ierr); 236 nel = pe - ps; 237 238 ierr = DMSwarmSetLocalSizes(dm,npoints_q*nel,-1);CHKERRQ(ierr); 239 ierr = DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr); 240 ierr = DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr); 241 242 ierr = DMGetCoordinatesLocal(dmc,&coorlocal);CHKERRQ(ierr); 243 ierr = DMGetCoordinateSection(dmc,&coordSection);CHKERRQ(ierr); 244 245 pcnt = 0; 246 for (e=0; e<nel; e++) { 247 ierr = DMPlexVecGetClosure(dmc,coordSection,coorlocal,e,NULL,&elcoor);CHKERRQ(ierr); 248 249 for (q=0; q<npoints_q; q++) { 250 for (d=0; d<dim; d++) { 251 swarm_coor[dim*pcnt+d] = 0.0; 252 for (k=0; k<npe; k++) { 253 swarm_coor[dim*pcnt+d] += basis[q][k] * PetscRealPart(elcoor[dim*k+d]); 254 } 255 } 256 swarm_cellid[pcnt] = e; 257 pcnt++; 258 } 259 ierr = DMPlexVecRestoreClosure(dmc,coordSection,coorlocal,e,NULL,&elcoor);CHKERRQ(ierr); 260 } 261 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr); 262 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr); 263 264 ierr = PetscFree(xi);CHKERRQ(ierr); 265 for (q=0; q<npoints_q; q++) { 266 ierr = PetscFree(basis[q]);CHKERRQ(ierr); 267 } 268 ierr = PetscFree(basis);CHKERRQ(ierr); 269 270 PetscFunctionReturn(0); 271 } 272 273 PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX(DM dm,DM celldm,DMSwarmPICLayoutType layout,PetscInt layout_param) 274 { 275 PetscErrorCode ierr; 276 PetscInt dim; 277 278 PetscFunctionBegin; 279 ierr = DMGetDimension(celldm,&dim);CHKERRQ(ierr); 280 switch (layout) { 281 case DMSWARMPIC_LAYOUT_REGULAR: 282 if (dim == 3) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No 3D support for REGULAR+PLEX"); 283 ierr = private_DMSwarmInsertPointsUsingCellDM_PLEX2D_Regular(dm,celldm,layout_param);CHKERRQ(ierr); 284 break; 285 case DMSWARMPIC_LAYOUT_GAUSS: 286 SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Gauss layout not supported for PLEX"); 287 break; 288 case DMSWARMPIC_LAYOUT_SUBDIVISION: 289 ierr = private_DMSwarmInsertPointsUsingCellDM_PLEX_SubDivide(dm,celldm,layout_param);CHKERRQ(ierr); 290 break; 291 } 292 PetscFunctionReturn(0); 293 } 294 295 /* 296 typedef struct { 297 PetscReal x,y; 298 } Point2d; 299 300 static PetscErrorCode signp2d(Point2d p1, Point2d p2, Point2d p3,PetscReal *s) 301 { 302 PetscFunctionBegin; 303 *s = (p1.x - p3.x) * (p2.y - p3.y) - (p2.x - p3.x) * (p1.y - p3.y); 304 PetscFunctionReturn(0); 305 } 306 */ 307 /* 308 static PetscErrorCode PointInTriangle(Point2d pt, Point2d v1, Point2d v2, Point2d v3,PetscBool *v) 309 { 310 PetscReal s1,s2,s3; 311 PetscBool b1, b2, b3; 312 313 PetscFunctionBegin; 314 signp2d(pt, v1, v2,&s1); b1 = s1 < 0.0f; 315 signp2d(pt, v2, v3,&s2); b2 = s2 < 0.0f; 316 signp2d(pt, v3, v1,&s3); b3 = s3 < 0.0f; 317 318 *v = ((b1 == b2) && (b2 == b3)); 319 PetscFunctionReturn(0); 320 } 321 */ 322 /* 323 static PetscErrorCode _ComputeLocalCoordinateAffine2d(PetscReal xp[],PetscReal coords[],PetscReal xip[],PetscReal *dJ) 324 { 325 PetscReal x1,y1,x2,y2,x3,y3; 326 PetscReal c,b[2],A[2][2],inv[2][2],detJ,od; 327 328 PetscFunctionBegin; 329 x1 = coords[2*0+0]; 330 x2 = coords[2*1+0]; 331 x3 = coords[2*2+0]; 332 333 y1 = coords[2*0+1]; 334 y2 = coords[2*1+1]; 335 y3 = coords[2*2+1]; 336 337 c = x1 - 0.5*x1 - 0.5*x1 + 0.5*x2 + 0.5*x3; 338 b[0] = xp[0] - c; 339 c = y1 - 0.5*y1 - 0.5*y1 + 0.5*y2 + 0.5*y3; 340 b[1] = xp[1] - c; 341 342 A[0][0] = -0.5*x1 + 0.5*x2; A[0][1] = -0.5*x1 + 0.5*x3; 343 A[1][0] = -0.5*y1 + 0.5*y2; A[1][1] = -0.5*y1 + 0.5*y3; 344 345 detJ = A[0][0]*A[1][1] - A[0][1]*A[1][0]; 346 *dJ = PetscAbsReal(detJ); 347 od = 1.0/detJ; 348 349 inv[0][0] = A[1][1] * od; 350 inv[0][1] = -A[0][1] * od; 351 inv[1][0] = -A[1][0] * od; 352 inv[1][1] = A[0][0] * od; 353 354 xip[0] = inv[0][0]*b[0] + inv[0][1]*b[1]; 355 xip[1] = inv[1][0]*b[0] + inv[1][1]*b[1]; 356 PetscFunctionReturn(0); 357 } 358 */ 359 360 static PetscErrorCode ComputeLocalCoordinateAffine2d(PetscReal xp[],PetscScalar coords[],PetscReal xip[],PetscReal *dJ) 361 { 362 PetscReal x1,y1,x2,y2,x3,y3; 363 PetscReal b[2],A[2][2],inv[2][2],detJ,od; 364 365 PetscFunctionBegin; 366 x1 = PetscRealPart(coords[2*0+0]); 367 x2 = PetscRealPart(coords[2*1+0]); 368 x3 = PetscRealPart(coords[2*2+0]); 369 370 y1 = PetscRealPart(coords[2*0+1]); 371 y2 = PetscRealPart(coords[2*1+1]); 372 y3 = PetscRealPart(coords[2*2+1]); 373 374 b[0] = xp[0] - x1; 375 b[1] = xp[1] - y1; 376 377 A[0][0] = x2-x1; A[0][1] = x3-x1; 378 A[1][0] = y2-y1; A[1][1] = y3-y1; 379 380 detJ = A[0][0]*A[1][1] - A[0][1]*A[1][0]; 381 *dJ = PetscAbsReal(detJ); 382 od = 1.0/detJ; 383 384 inv[0][0] = A[1][1] * od; 385 inv[0][1] = -A[0][1] * od; 386 inv[1][0] = -A[1][0] * od; 387 inv[1][1] = A[0][0] * od; 388 389 xip[0] = inv[0][0]*b[0] + inv[0][1]*b[1]; 390 xip[1] = inv[1][0]*b[0] + inv[1][1]*b[1]; 391 PetscFunctionReturn(0); 392 } 393 394 PetscErrorCode DMSwarmProjectField_ApproxP1_PLEX_2D(DM swarm,PetscReal *swarm_field,DM dm,Vec v_field) 395 { 396 PetscErrorCode ierr; 397 const PetscReal PLEX_C_EPS = 1.0e-8; 398 Vec v_field_l,denom_l,coor_l,denom; 399 PetscInt k,p,e,npoints; 400 PetscInt *mpfield_cell; 401 PetscReal *mpfield_coor; 402 PetscReal xi_p[2]; 403 PetscScalar Ni[3]; 404 PetscSection coordSection; 405 PetscScalar *elcoor = NULL; 406 407 PetscFunctionBegin; 408 ierr = VecZeroEntries(v_field);CHKERRQ(ierr); 409 410 ierr = DMGetLocalVector(dm,&v_field_l);CHKERRQ(ierr); 411 ierr = DMGetGlobalVector(dm,&denom);CHKERRQ(ierr); 412 ierr = DMGetLocalVector(dm,&denom_l);CHKERRQ(ierr); 413 ierr = VecZeroEntries(v_field_l);CHKERRQ(ierr); 414 ierr = VecZeroEntries(denom);CHKERRQ(ierr); 415 ierr = VecZeroEntries(denom_l);CHKERRQ(ierr); 416 417 ierr = DMGetCoordinatesLocal(dm,&coor_l);CHKERRQ(ierr); 418 ierr = DMGetCoordinateSection(dm,&coordSection);CHKERRQ(ierr); 419 420 ierr = DMSwarmGetLocalSize(swarm,&npoints);CHKERRQ(ierr); 421 ierr = DMSwarmGetField(swarm,DMSwarmPICField_coor,NULL,NULL,(void**)&mpfield_coor);CHKERRQ(ierr); 422 ierr = DMSwarmGetField(swarm,DMSwarmPICField_cellid,NULL,NULL,(void**)&mpfield_cell);CHKERRQ(ierr); 423 424 for (p=0; p<npoints; p++) { 425 PetscReal *coor_p,dJ; 426 PetscScalar elfield[3]; 427 PetscBool point_located; 428 429 e = mpfield_cell[p]; 430 coor_p = &mpfield_coor[2*p]; 431 432 ierr = DMPlexVecGetClosure(dm,coordSection,coor_l,e,NULL,&elcoor);CHKERRQ(ierr); 433 434 /* 435 while (!point_located && (failed_counter < 25)) { 436 ierr = PointInTriangle(point, coords[0], coords[1], coords[2], &point_located);CHKERRQ(ierr); 437 point.x = coor_p[0]; 438 point.y = coor_p[1]; 439 point.x += 1.0e-10 * (2.0 * rand()/((double)RAND_MAX)-1.0); 440 point.y += 1.0e-10 * (2.0 * rand()/((double)RAND_MAX)-1.0); 441 failed_counter++; 442 } 443 444 if (!point_located){ 445 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); 446 } 447 448 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); 449 else { 450 ierr = _ComputeLocalCoordinateAffine2d(coor_p,elcoor,xi_p,&dJ);CHKERRQ(ierr); 451 xi_p[0] = 0.5*(xi_p[0] + 1.0); 452 xi_p[1] = 0.5*(xi_p[1] + 1.0); 453 454 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]); 455 456 } 457 */ 458 459 ierr = ComputeLocalCoordinateAffine2d(coor_p,elcoor,xi_p,&dJ);CHKERRQ(ierr); 460 /* 461 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]); 462 */ 463 /* 464 point_located = PETSC_TRUE; 465 if (xi_p[0] < 0.0) { 466 if (xi_p[0] > -PLEX_C_EPS) { 467 xi_p[0] = 0.0; 468 } else { 469 point_located = PETSC_FALSE; 470 } 471 } 472 if (xi_p[1] < 0.0) { 473 if (xi_p[1] > -PLEX_C_EPS) { 474 xi_p[1] = 0.0; 475 } else { 476 point_located = PETSC_FALSE; 477 } 478 } 479 if (xi_p[1] > (1.0-xi_p[0])) { 480 if ((xi_p[1] - 1.0 + xi_p[0]) < PLEX_C_EPS) { 481 xi_p[1] = 1.0 - xi_p[0]; 482 } else { 483 point_located = PETSC_FALSE; 484 } 485 } 486 if (!point_located){ 487 PetscPrintf(PETSC_COMM_SELF,"[Error] xi,eta = %+1.8e, %+1.8e\n",xi_p[0],xi_p[1]); 488 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]); 489 } 490 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); 491 */ 492 493 Ni[0] = 1.0 - xi_p[0] - xi_p[1]; 494 Ni[1] = xi_p[0]; 495 Ni[2] = xi_p[1]; 496 497 point_located = PETSC_TRUE; 498 for (k=0; k<3; k++) { 499 if (PetscRealPart(Ni[k]) < -PLEX_C_EPS) point_located = PETSC_FALSE; 500 if (PetscRealPart(Ni[k]) > (1.0+PLEX_C_EPS)) point_located = PETSC_FALSE; 501 } 502 if (!point_located){ 503 PetscPrintf(PETSC_COMM_SELF,"[Error] xi,eta = %+1.8e, %+1.8e\n",(double)xi_p[0],(double)xi_p[1]); 504 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])); 505 } 506 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); 507 508 for (k=0; k<3; k++) { 509 Ni[k] = Ni[k] * dJ; 510 elfield[k] = Ni[k] * swarm_field[p]; 511 } 512 513 ierr = DMPlexVecRestoreClosure(dm,coordSection,coor_l,e,NULL,&elcoor);CHKERRQ(ierr); 514 515 ierr = DMPlexVecSetClosure(dm, NULL,v_field_l, e, elfield, ADD_VALUES);CHKERRQ(ierr); 516 ierr = DMPlexVecSetClosure(dm, NULL,denom_l, e, Ni, ADD_VALUES);CHKERRQ(ierr); 517 } 518 519 ierr = DMSwarmRestoreField(swarm,DMSwarmPICField_cellid,NULL,NULL,(void**)&mpfield_cell);CHKERRQ(ierr); 520 ierr = DMSwarmRestoreField(swarm,DMSwarmPICField_coor,NULL,NULL,(void**)&mpfield_coor);CHKERRQ(ierr); 521 522 ierr = DMLocalToGlobalBegin(dm,v_field_l,ADD_VALUES,v_field);CHKERRQ(ierr); 523 ierr = DMLocalToGlobalEnd(dm,v_field_l,ADD_VALUES,v_field);CHKERRQ(ierr); 524 ierr = DMLocalToGlobalBegin(dm,denom_l,ADD_VALUES,denom);CHKERRQ(ierr); 525 ierr = DMLocalToGlobalEnd(dm,denom_l,ADD_VALUES,denom);CHKERRQ(ierr); 526 527 ierr = VecPointwiseDivide(v_field,v_field,denom);CHKERRQ(ierr); 528 529 ierr = DMRestoreLocalVector(dm,&v_field_l);CHKERRQ(ierr); 530 ierr = DMRestoreLocalVector(dm,&denom_l);CHKERRQ(ierr); 531 ierr = DMRestoreGlobalVector(dm,&denom);CHKERRQ(ierr); 532 533 PetscFunctionReturn(0); 534 } 535 536 PetscErrorCode private_DMSwarmProjectFields_PLEX(DM swarm,DM celldm,PetscInt project_type,PetscInt nfields,DataField dfield[],Vec vecs[]) 537 { 538 PetscErrorCode ierr; 539 PetscInt f,dim; 540 541 PetscFunctionBegin; 542 ierr = DMGetDimension(swarm,&dim);CHKERRQ(ierr); 543 switch (dim) { 544 case 2: 545 for (f=0; f<nfields; f++) { 546 PetscReal *swarm_field; 547 548 ierr = DataFieldGetEntries(dfield[f],(void**)&swarm_field);CHKERRQ(ierr); 549 ierr = DMSwarmProjectField_ApproxP1_PLEX_2D(swarm,swarm_field,celldm,vecs[f]);CHKERRQ(ierr); 550 } 551 break; 552 case 3: 553 SETERRQ(PetscObjectComm((PetscObject)swarm),PETSC_ERR_SUP,"No support for 3D"); 554 break; 555 default: 556 break; 557 } 558 559 PetscFunctionReturn(0); 560 } 561