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