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] * 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] * 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[],PetscReal 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 = coords[2*0+0]; 297 x2 = coords[2*1+0]; 298 x3 = coords[2*2+0]; 299 300 y1 = coords[2*0+1]; 301 y2 = coords[2*1+1]; 302 y3 = 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],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; 355 PetscReal elfield[3],dJ; 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 // rescale to [0,1] // 381 xi_p[0] = 0.5*(xi_p[0] + 1.0); 382 xi_p[1] = 0.5*(xi_p[1] + 1.0); 383 384 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]); 385 386 } 387 */ 388 389 ierr = ComputeLocalCoordinateAffine2d(coor_p,elcoor,xi_p,&dJ);CHKERRQ(ierr); 390 /* 391 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]); 392 */ 393 /* 394 point_located = PETSC_TRUE; 395 if (xi_p[0] < 0.0) { 396 if (xi_p[0] > -PLEX_C_EPS) { 397 xi_p[0] = 0.0; 398 } else { 399 point_located = PETSC_FALSE; 400 } 401 } 402 if (xi_p[1] < 0.0) { 403 if (xi_p[1] > -PLEX_C_EPS) { 404 xi_p[1] = 0.0; 405 } else { 406 point_located = PETSC_FALSE; 407 } 408 } 409 if (xi_p[1] > (1.0-xi_p[0])) { 410 if ((xi_p[1] - 1.0 + xi_p[0]) < PLEX_C_EPS) { 411 xi_p[1] = 1.0 - xi_p[0]; 412 } else { 413 point_located = PETSC_FALSE; 414 } 415 } 416 if (!point_located){ 417 PetscPrintf(PETSC_COMM_SELF,"[Error] xi,eta = %+1.8e, %+1.8e\n",xi_p[0],xi_p[1]); 418 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]); 419 } 420 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); 421 */ 422 423 Ni[0] = 1.0 - xi_p[0] - xi_p[1]; 424 Ni[1] = xi_p[0]; 425 Ni[2] = xi_p[1]; 426 427 point_located = PETSC_TRUE; 428 for (k=0; k<3; k++) { 429 if (Ni[k] < -PLEX_C_EPS) point_located = PETSC_FALSE; 430 if (Ni[k] > (1.0+PLEX_C_EPS)) point_located = PETSC_FALSE; 431 } 432 if (!point_located){ 433 PetscPrintf(PETSC_COMM_SELF,"[Error] xi,eta = %+1.8e, %+1.8e\n",xi_p[0],xi_p[1]); 434 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]); 435 } 436 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); 437 438 for (k=0; k<3; k++) { 439 Ni[k] = Ni[k] * dJ; 440 elfield[k] = Ni[k] * swarm_field[p]; 441 } 442 443 ierr = DMPlexVecRestoreClosure(dm,coordSection,coor_l,e,NULL,&elcoor);CHKERRQ(ierr); 444 445 ierr = DMPlexVecSetClosure(dm, NULL,v_field_l, e, elfield, ADD_VALUES);CHKERRQ(ierr); 446 ierr = DMPlexVecSetClosure(dm, NULL,denom_l, e, Ni, ADD_VALUES);CHKERRQ(ierr); 447 } 448 449 ierr = DMSwarmRestoreField(swarm,DMSwarmPICField_cellid,NULL,NULL,(void**)&mpfield_cell);CHKERRQ(ierr); 450 ierr = DMSwarmRestoreField(swarm,DMSwarmPICField_coor,NULL,NULL,(void**)&mpfield_coor);CHKERRQ(ierr); 451 452 ierr = DMLocalToGlobalBegin(dm,v_field_l,ADD_VALUES,v_field);CHKERRQ(ierr); 453 ierr = DMLocalToGlobalEnd(dm,v_field_l,ADD_VALUES,v_field);CHKERRQ(ierr); 454 ierr = DMLocalToGlobalBegin(dm,denom_l,ADD_VALUES,denom);CHKERRQ(ierr); 455 ierr = DMLocalToGlobalEnd(dm,denom_l,ADD_VALUES,denom);CHKERRQ(ierr); 456 457 ierr = VecPointwiseDivide(v_field,v_field,denom);CHKERRQ(ierr); 458 459 ierr = DMRestoreLocalVector(dm,&v_field_l);CHKERRQ(ierr); 460 ierr = DMRestoreLocalVector(dm,&denom_l);CHKERRQ(ierr); 461 ierr = DMRestoreGlobalVector(dm,&denom);CHKERRQ(ierr); 462 463 PetscFunctionReturn(0); 464 } 465 466 PetscErrorCode private_DMSwarmProjectFields_PLEX(DM swarm,DM celldm,PetscInt project_type,PetscInt nfields,DataField dfield[],Vec vecs[]) 467 { 468 PetscErrorCode ierr; 469 PetscInt f,dim; 470 471 PetscFunctionBegin; 472 ierr = DMGetDimension(swarm,&dim);CHKERRQ(ierr); 473 switch (dim) { 474 case 2: 475 for (f=0; f<nfields; f++) { 476 PetscReal *swarm_field; 477 478 ierr = DataFieldGetEntries(dfield[f],(void**)&swarm_field);CHKERRQ(ierr); 479 ierr = DMSwarmProjectField_ApproxP1_PLEX_2D(swarm,swarm_field,celldm,vecs[f]);CHKERRQ(ierr); 480 } 481 break; 482 case 3: 483 SETERRQ(PetscObjectComm((PetscObject)swarm),PETSC_ERR_SUP,"No support for 3D"); 484 break; 485 default: 486 break; 487 } 488 489 PetscFunctionReturn(0); 490 } 491