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