10e2ec84fSDave May 20e2ec84fSDave May #include <petsc.h> 30e2ec84fSDave May #include <petscdm.h> 40e2ec84fSDave May #include <petscdmplex.h> 50e2ec84fSDave May #include <petscdmswarm.h> 6*e3cd5995SDave May #include "data_bucket.h" 70e2ec84fSDave May 80e2ec84fSDave May PetscErrorCode subdivide_triangle(PetscReal v1[2],PetscReal v2[2],PetscReal v3[3],PetscInt depth,PetscInt max,PetscReal xi[],PetscInt *np) 90e2ec84fSDave May { 100e2ec84fSDave May PetscReal v12[2],v23[2],v31[2]; 110e2ec84fSDave May PetscInt i; 120e2ec84fSDave May PetscErrorCode ierr; 130e2ec84fSDave May 140e2ec84fSDave May PetscFunctionBegin; 150e2ec84fSDave May if (depth == max) { 160e2ec84fSDave May PetscReal cx[2]; 170e2ec84fSDave May 180e2ec84fSDave May //printf("centroid \n"); 190e2ec84fSDave May cx[0] = (v1[0] + v2[0] + v3[0])/3.0; 200e2ec84fSDave May cx[1] = (v1[1] + v2[1] + v3[1])/3.0; 210e2ec84fSDave May 220e2ec84fSDave May xi[2*(*np)+0] = cx[0]; 230e2ec84fSDave May xi[2*(*np)+1] = cx[1]; 240e2ec84fSDave May (*np)++; 250e2ec84fSDave May PetscFunctionReturn(0); 260e2ec84fSDave May } 270e2ec84fSDave May 280e2ec84fSDave May /* calculate midpoints of each side */ 290e2ec84fSDave May for (i = 0; i < 2; i++) { 300e2ec84fSDave May v12[i] = (v1[i]+v2[i])/2.0; 310e2ec84fSDave May v23[i] = (v2[i]+v3[i])/2.0; 320e2ec84fSDave May v31[i] = (v3[i]+v1[i])/2.0; 330e2ec84fSDave May } 340e2ec84fSDave May 350e2ec84fSDave May /* recursively subdivide new triangles */ 360e2ec84fSDave May ierr = subdivide_triangle(v1,v12,v31,depth+1,max,xi,np);CHKERRQ(ierr); 370e2ec84fSDave May ierr = subdivide_triangle(v2,v23,v12,depth+1,max,xi,np);CHKERRQ(ierr); 380e2ec84fSDave May ierr = subdivide_triangle(v3,v31,v23,depth+1,max,xi,np);CHKERRQ(ierr); 390e2ec84fSDave May ierr = subdivide_triangle(v12,v23,v31,depth+1,max,xi,np);CHKERRQ(ierr); 400e2ec84fSDave May PetscFunctionReturn(0); 410e2ec84fSDave May } 420e2ec84fSDave May 430e2ec84fSDave May 440e2ec84fSDave May PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX2D_SubDivide(DM dm,DM dmc,PetscInt nsub) 450e2ec84fSDave May { 460e2ec84fSDave May PetscErrorCode ierr; 470e2ec84fSDave May const PetscInt dim = 2; 480e2ec84fSDave May PetscInt q,npoints_q,e,nel,npe,pcnt,ps,pe,d,k,depth; 490e2ec84fSDave May PetscReal *xi; 500e2ec84fSDave May PetscReal **basis; 510e2ec84fSDave May Vec coorlocal; 520e2ec84fSDave May PetscSection coordSection; 530e2ec84fSDave May PetscScalar *elcoor = NULL; 540e2ec84fSDave May PetscReal *swarm_coor; 550e2ec84fSDave May PetscInt *swarm_cellid; 560e2ec84fSDave May PetscReal v1[2],v2[2],v3[2]; 570e2ec84fSDave May 580e2ec84fSDave May PetscFunctionBegin; 590e2ec84fSDave May 600e2ec84fSDave May npoints_q = 1; 610e2ec84fSDave May for (d=0; d<nsub; d++) { npoints_q *= 4; } 620e2ec84fSDave May ierr = PetscMalloc1(dim*npoints_q,&xi);CHKERRQ(ierr); 630e2ec84fSDave May 640e2ec84fSDave May v1[0] = 0.0; v1[1] = 0.0; 650e2ec84fSDave May v2[0] = 1.0; v2[1] = 0.0; 660e2ec84fSDave May v3[0] = 0.0; v3[1] = 1.0; 670e2ec84fSDave May depth = 0; 680e2ec84fSDave May pcnt = 0; 690e2ec84fSDave May ierr = subdivide_triangle(v1,v2,v3,depth,nsub,xi,&pcnt);CHKERRQ(ierr); 700e2ec84fSDave May 710e2ec84fSDave May npe = 3; /* nodes per element (triangle) */ 720e2ec84fSDave May ierr = PetscMalloc1(npoints_q,&basis);CHKERRQ(ierr); 730e2ec84fSDave May for (q=0; q<npoints_q; q++) { 740e2ec84fSDave May ierr = PetscMalloc1(npe,&basis[q]);CHKERRQ(ierr); 750e2ec84fSDave May 760e2ec84fSDave May basis[q][0] = 1.0 - xi[dim*q+0] - xi[dim*q+1]; 770e2ec84fSDave May basis[q][1] = xi[dim*q+0]; 780e2ec84fSDave May basis[q][2] = xi[dim*q+1]; 790e2ec84fSDave May } 800e2ec84fSDave May 810e2ec84fSDave May /* 0->cell, 1->edge, 2->vert */ 820e2ec84fSDave May ierr = DMPlexGetHeightStratum(dmc,0,&ps,&pe);CHKERRQ(ierr); 830e2ec84fSDave May nel = pe - ps; 840e2ec84fSDave May 850e2ec84fSDave May ierr = DMSwarmSetLocalSizes(dm,npoints_q*nel,-1);CHKERRQ(ierr); 860e2ec84fSDave May ierr = DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr); 870e2ec84fSDave May ierr = DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr); 880e2ec84fSDave May 890e2ec84fSDave May ierr = DMGetCoordinatesLocal(dmc,&coorlocal);CHKERRQ(ierr); 900e2ec84fSDave May ierr = DMGetCoordinateSection(dmc,&coordSection);CHKERRQ(ierr); 910e2ec84fSDave May 920e2ec84fSDave May pcnt = 0; 930e2ec84fSDave May for (e=0; e<nel; e++) { 940e2ec84fSDave May ierr = DMPlexVecGetClosure(dmc,coordSection,coorlocal,e,NULL,&elcoor);CHKERRQ(ierr); 950e2ec84fSDave May 960e2ec84fSDave May for (q=0; q<npoints_q; q++) { 970e2ec84fSDave May for (d=0; d<dim; d++) { 980e2ec84fSDave May swarm_coor[dim*pcnt+d] = 0.0; 990e2ec84fSDave May for (k=0; k<npe; k++) { 1000e2ec84fSDave May swarm_coor[dim*pcnt+d] += basis[q][k] * elcoor[dim*k+d]; 1010e2ec84fSDave May } 1020e2ec84fSDave May } 1030e2ec84fSDave May swarm_cellid[pcnt] = e; 1040e2ec84fSDave May pcnt++; 1050e2ec84fSDave May } 1060e2ec84fSDave May ierr = DMPlexVecRestoreClosure(dmc,coordSection,coorlocal,e,NULL,&elcoor);CHKERRQ(ierr); 1070e2ec84fSDave May } 1080e2ec84fSDave May ierr = DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr); 1090e2ec84fSDave May ierr = DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr); 1100e2ec84fSDave May 1110e2ec84fSDave May ierr = PetscFree(xi);CHKERRQ(ierr); 1120e2ec84fSDave May for (q=0; q<npoints_q; q++) { 1130e2ec84fSDave May ierr = PetscFree(basis[q]);CHKERRQ(ierr); 1140e2ec84fSDave May } 1150e2ec84fSDave May ierr = PetscFree(basis);CHKERRQ(ierr); 1160e2ec84fSDave May 1170e2ec84fSDave May PetscFunctionReturn(0); 1180e2ec84fSDave May } 1190e2ec84fSDave May 1200e2ec84fSDave May PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX2D_Regular(DM dm,DM dmc,PetscInt npoints) 1210e2ec84fSDave May { 1220e2ec84fSDave May PetscErrorCode ierr; 1230e2ec84fSDave May const PetscInt dim = 2; 1240e2ec84fSDave May PetscInt ii,jj,q,npoints_q,e,nel,npe,pcnt,ps,pe,d,k; 1250e2ec84fSDave May PetscReal *xi,ds,ds2; 1260e2ec84fSDave May PetscReal **basis; 1270e2ec84fSDave May Vec coorlocal; 1280e2ec84fSDave May PetscSection coordSection; 1290e2ec84fSDave May PetscScalar *elcoor = NULL; 1300e2ec84fSDave May PetscReal *swarm_coor; 1310e2ec84fSDave May PetscInt *swarm_cellid; 1320e2ec84fSDave May 1330e2ec84fSDave May PetscFunctionBegin; 1340e2ec84fSDave May ierr = PetscMalloc1(dim*npoints*npoints,&xi);CHKERRQ(ierr); 1350e2ec84fSDave May pcnt = 0; 1360e2ec84fSDave May ds = 1.0/((PetscReal)(npoints-1)); 1370e2ec84fSDave May ds2 = 1.0/((PetscReal)(npoints)); 1380e2ec84fSDave May for (jj = 0; jj<npoints; jj++) { 1390e2ec84fSDave May for (ii=0; ii<npoints-jj; ii++) { 1400e2ec84fSDave May xi[dim*pcnt+0] = ii * ds; 1410e2ec84fSDave May xi[dim*pcnt+1] = jj * ds; 1420e2ec84fSDave May 1430e2ec84fSDave May xi[dim*pcnt+0] *= (1.0 - 1.2*ds2); 1440e2ec84fSDave May xi[dim*pcnt+1] *= (1.0 - 1.2*ds2); 1450e2ec84fSDave May 1460e2ec84fSDave May xi[dim*pcnt+0] += 0.35*ds2; 1470e2ec84fSDave May xi[dim*pcnt+1] += 0.35*ds2; 1480e2ec84fSDave May 1490e2ec84fSDave May pcnt++; 1500e2ec84fSDave May } 1510e2ec84fSDave May } 1520e2ec84fSDave May npoints_q = pcnt; 1530e2ec84fSDave May 1540e2ec84fSDave May npe = 3; /* nodes per element (triangle) */ 1550e2ec84fSDave May ierr = PetscMalloc1(npoints_q,&basis);CHKERRQ(ierr); 1560e2ec84fSDave May for (q=0; q<npoints_q; q++) { 1570e2ec84fSDave May ierr = PetscMalloc1(npe,&basis[q]);CHKERRQ(ierr); 1580e2ec84fSDave May 1590e2ec84fSDave May basis[q][0] = 1.0 - xi[dim*q+0] - xi[dim*q+1]; 1600e2ec84fSDave May basis[q][1] = xi[dim*q+0]; 1610e2ec84fSDave May basis[q][2] = xi[dim*q+1]; 1620e2ec84fSDave May } 1630e2ec84fSDave May 1640e2ec84fSDave May /* 0->cell, 1->edge, 2->vert */ 1650e2ec84fSDave May ierr = DMPlexGetHeightStratum(dmc,0,&ps,&pe);CHKERRQ(ierr); 1660e2ec84fSDave May nel = pe - ps; 1670e2ec84fSDave May 1680e2ec84fSDave May ierr = DMSwarmSetLocalSizes(dm,npoints_q*nel,-1);CHKERRQ(ierr); 1690e2ec84fSDave May ierr = DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr); 1700e2ec84fSDave May ierr = DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr); 1710e2ec84fSDave May 1720e2ec84fSDave May ierr = DMGetCoordinatesLocal(dmc,&coorlocal);CHKERRQ(ierr); 1730e2ec84fSDave May ierr = DMGetCoordinateSection(dmc,&coordSection);CHKERRQ(ierr); 1740e2ec84fSDave May 1750e2ec84fSDave May pcnt = 0; 1760e2ec84fSDave May for (e=0; e<nel; e++) { 1770e2ec84fSDave May ierr = DMPlexVecGetClosure(dmc,coordSection,coorlocal,e,NULL,&elcoor);CHKERRQ(ierr); 1780e2ec84fSDave May 1790e2ec84fSDave May for (q=0; q<npoints_q; q++) { 1800e2ec84fSDave May for (d=0; d<dim; d++) { 1810e2ec84fSDave May swarm_coor[dim*pcnt+d] = 0.0; 1820e2ec84fSDave May for (k=0; k<npe; k++) { 1830e2ec84fSDave May swarm_coor[dim*pcnt+d] += basis[q][k] * elcoor[dim*k+d]; 1840e2ec84fSDave May } 1850e2ec84fSDave May } 1860e2ec84fSDave May swarm_cellid[pcnt] = e; 1870e2ec84fSDave May pcnt++; 1880e2ec84fSDave May } 1890e2ec84fSDave May ierr = DMPlexVecRestoreClosure(dmc,coordSection,coorlocal,e,NULL,&elcoor);CHKERRQ(ierr); 1900e2ec84fSDave May } 1910e2ec84fSDave May ierr = DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr); 1920e2ec84fSDave May ierr = DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr); 1930e2ec84fSDave May 1940e2ec84fSDave May ierr = PetscFree(xi);CHKERRQ(ierr); 1950e2ec84fSDave May for (q=0; q<npoints_q; q++) { 1960e2ec84fSDave May ierr = PetscFree(basis[q]);CHKERRQ(ierr); 1970e2ec84fSDave May } 1980e2ec84fSDave May ierr = PetscFree(basis);CHKERRQ(ierr); 1990e2ec84fSDave May 2000e2ec84fSDave May PetscFunctionReturn(0); 2010e2ec84fSDave May } 2020e2ec84fSDave May 2030e2ec84fSDave May PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX(DM dm,DM celldm,DMSwarmPICLayoutType layout,PetscInt layout_param) 2040e2ec84fSDave May { 2050e2ec84fSDave May PetscErrorCode ierr; 2060e2ec84fSDave May PetscInt dim; 2070e2ec84fSDave May 2080e2ec84fSDave May PetscFunctionBegin; 2090e2ec84fSDave May ierr = DMGetDimension(celldm,&dim);CHKERRQ(ierr); 2100e2ec84fSDave May if (dim == 3) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No 3D support for PLEX"); 2110e2ec84fSDave May switch (layout) { 2120e2ec84fSDave May case DMSWARMPIC_LAYOUT_REGULAR: 2130e2ec84fSDave May ierr = private_DMSwarmInsertPointsUsingCellDM_PLEX2D_Regular(dm,celldm,layout_param);CHKERRQ(ierr); 2140e2ec84fSDave May break; 2150e2ec84fSDave May case DMSWARMPIC_LAYOUT_GAUSS: 2160e2ec84fSDave May SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Gauss layout not supported for PLEX"); 2170e2ec84fSDave May break; 2180e2ec84fSDave May case DMSWARMPIC_LAYOUT_SUBDIVISION: 2190e2ec84fSDave May ierr = private_DMSwarmInsertPointsUsingCellDM_PLEX2D_SubDivide(dm,celldm,layout_param);CHKERRQ(ierr); 2200e2ec84fSDave May break; 2210e2ec84fSDave May } 2220e2ec84fSDave May PetscFunctionReturn(0); 2230e2ec84fSDave May } 224*e3cd5995SDave May 225*e3cd5995SDave May /* 226*e3cd5995SDave May typedef struct { 227*e3cd5995SDave May PetscReal x,y; 228*e3cd5995SDave May } Point2d; 229*e3cd5995SDave May 230*e3cd5995SDave May static PetscErrorCode signp2d(Point2d p1, Point2d p2, Point2d p3,PetscReal *s) 231*e3cd5995SDave May { 232*e3cd5995SDave May PetscFunctionBegin; 233*e3cd5995SDave May *s = (p1.x - p3.x) * (p2.y - p3.y) - (p2.x - p3.x) * (p1.y - p3.y); 234*e3cd5995SDave May PetscFunctionReturn(0); 235*e3cd5995SDave May } 236*e3cd5995SDave May */ 237*e3cd5995SDave May /* 238*e3cd5995SDave May static PetscErrorCode PointInTriangle(Point2d pt, Point2d v1, Point2d v2, Point2d v3,PetscBool *v) 239*e3cd5995SDave May { 240*e3cd5995SDave May PetscReal s1,s2,s3; 241*e3cd5995SDave May PetscBool b1, b2, b3; 242*e3cd5995SDave May 243*e3cd5995SDave May PetscFunctionBegin; 244*e3cd5995SDave May signp2d(pt, v1, v2,&s1); b1 = s1 < 0.0f; 245*e3cd5995SDave May signp2d(pt, v2, v3,&s2); b2 = s2 < 0.0f; 246*e3cd5995SDave May signp2d(pt, v3, v1,&s3); b3 = s3 < 0.0f; 247*e3cd5995SDave May 248*e3cd5995SDave May *v = ((b1 == b2) && (b2 == b3)); 249*e3cd5995SDave May PetscFunctionReturn(0); 250*e3cd5995SDave May } 251*e3cd5995SDave May */ 252*e3cd5995SDave May /* 253*e3cd5995SDave May static PetscErrorCode _ComputeLocalCoordinateAffine2d(PetscReal xp[],PetscReal coords[],PetscReal xip[],PetscReal *dJ) 254*e3cd5995SDave May { 255*e3cd5995SDave May PetscReal x1,y1,x2,y2,x3,y3; 256*e3cd5995SDave May PetscReal c,b[2],A[2][2],inv[2][2],detJ,od; 257*e3cd5995SDave May 258*e3cd5995SDave May PetscFunctionBegin; 259*e3cd5995SDave May x1 = coords[2*0+0]; 260*e3cd5995SDave May x2 = coords[2*1+0]; 261*e3cd5995SDave May x3 = coords[2*2+0]; 262*e3cd5995SDave May 263*e3cd5995SDave May y1 = coords[2*0+1]; 264*e3cd5995SDave May y2 = coords[2*1+1]; 265*e3cd5995SDave May y3 = coords[2*2+1]; 266*e3cd5995SDave May 267*e3cd5995SDave May c = x1 - 0.5*x1 - 0.5*x1 + 0.5*x2 + 0.5*x3; 268*e3cd5995SDave May b[0] = xp[0] - c; 269*e3cd5995SDave May c = y1 - 0.5*y1 - 0.5*y1 + 0.5*y2 + 0.5*y3; 270*e3cd5995SDave May b[1] = xp[1] - c; 271*e3cd5995SDave May 272*e3cd5995SDave May A[0][0] = -0.5*x1 + 0.5*x2; A[0][1] = -0.5*x1 + 0.5*x3; 273*e3cd5995SDave May A[1][0] = -0.5*y1 + 0.5*y2; A[1][1] = -0.5*y1 + 0.5*y3; 274*e3cd5995SDave May 275*e3cd5995SDave May detJ = A[0][0]*A[1][1] - A[0][1]*A[1][0]; 276*e3cd5995SDave May *dJ = PetscAbsReal(detJ); 277*e3cd5995SDave May od = 1.0/detJ; 278*e3cd5995SDave May 279*e3cd5995SDave May inv[0][0] = A[1][1] * od; 280*e3cd5995SDave May inv[0][1] = -A[0][1] * od; 281*e3cd5995SDave May inv[1][0] = -A[1][0] * od; 282*e3cd5995SDave May inv[1][1] = A[0][0] * od; 283*e3cd5995SDave May 284*e3cd5995SDave May xip[0] = inv[0][0]*b[0] + inv[0][1]*b[1]; 285*e3cd5995SDave May xip[1] = inv[1][0]*b[0] + inv[1][1]*b[1]; 286*e3cd5995SDave May PetscFunctionReturn(0); 287*e3cd5995SDave May } 288*e3cd5995SDave May */ 289*e3cd5995SDave May 290*e3cd5995SDave May static PetscErrorCode ComputeLocalCoordinateAffine2d(PetscReal xp[],PetscReal coords[],PetscReal xip[],PetscReal *dJ) 291*e3cd5995SDave May { 292*e3cd5995SDave May PetscReal x1,y1,x2,y2,x3,y3; 293*e3cd5995SDave May PetscReal b[2],A[2][2],inv[2][2],detJ,od; 294*e3cd5995SDave May 295*e3cd5995SDave May PetscFunctionBegin; 296*e3cd5995SDave May x1 = coords[2*0+0]; 297*e3cd5995SDave May x2 = coords[2*1+0]; 298*e3cd5995SDave May x3 = coords[2*2+0]; 299*e3cd5995SDave May 300*e3cd5995SDave May y1 = coords[2*0+1]; 301*e3cd5995SDave May y2 = coords[2*1+1]; 302*e3cd5995SDave May y3 = coords[2*2+1]; 303*e3cd5995SDave May 304*e3cd5995SDave May b[0] = xp[0] - x1; 305*e3cd5995SDave May b[1] = xp[1] - y1; 306*e3cd5995SDave May 307*e3cd5995SDave May A[0][0] = x2-x1; A[0][1] = x3-x1; 308*e3cd5995SDave May A[1][0] = y2-y1; A[1][1] = y3-y1; 309*e3cd5995SDave May 310*e3cd5995SDave May detJ = A[0][0]*A[1][1] - A[0][1]*A[1][0]; 311*e3cd5995SDave May *dJ = PetscAbsReal(detJ); 312*e3cd5995SDave May od = 1.0/detJ; 313*e3cd5995SDave May 314*e3cd5995SDave May inv[0][0] = A[1][1] * od; 315*e3cd5995SDave May inv[0][1] = -A[0][1] * od; 316*e3cd5995SDave May inv[1][0] = -A[1][0] * od; 317*e3cd5995SDave May inv[1][1] = A[0][0] * od; 318*e3cd5995SDave May 319*e3cd5995SDave May xip[0] = inv[0][0]*b[0] + inv[0][1]*b[1]; 320*e3cd5995SDave May xip[1] = inv[1][0]*b[0] + inv[1][1]*b[1]; 321*e3cd5995SDave May PetscFunctionReturn(0); 322*e3cd5995SDave May } 323*e3cd5995SDave May 324*e3cd5995SDave May PetscErrorCode DMSwarmProjectField_ApproxP1_PLEX_2D(DM swarm,PetscReal *swarm_field,DM dm,Vec v_field) 325*e3cd5995SDave May { 326*e3cd5995SDave May PetscErrorCode ierr; 327*e3cd5995SDave May const PetscReal PLEX_C_EPS = 1.0e-8; 328*e3cd5995SDave May Vec v_field_l,denom_l,coor_l,denom; 329*e3cd5995SDave May PetscInt k,p,e,npoints; 330*e3cd5995SDave May PetscInt *mpfield_cell; 331*e3cd5995SDave May PetscReal *mpfield_coor; 332*e3cd5995SDave May PetscReal xi_p[2],Ni[3]; 333*e3cd5995SDave May PetscSection coordSection; 334*e3cd5995SDave May PetscScalar *elcoor = NULL; 335*e3cd5995SDave May 336*e3cd5995SDave May PetscFunctionBegin; 337*e3cd5995SDave May ierr = VecZeroEntries(v_field);CHKERRQ(ierr); 338*e3cd5995SDave May 339*e3cd5995SDave May ierr = DMGetLocalVector(dm,&v_field_l);CHKERRQ(ierr); 340*e3cd5995SDave May ierr = DMGetGlobalVector(dm,&denom);CHKERRQ(ierr); 341*e3cd5995SDave May ierr = DMGetLocalVector(dm,&denom_l);CHKERRQ(ierr); 342*e3cd5995SDave May ierr = VecZeroEntries(v_field_l);CHKERRQ(ierr); 343*e3cd5995SDave May ierr = VecZeroEntries(denom);CHKERRQ(ierr); 344*e3cd5995SDave May ierr = VecZeroEntries(denom_l);CHKERRQ(ierr); 345*e3cd5995SDave May 346*e3cd5995SDave May ierr = DMGetCoordinatesLocal(dm,&coor_l);CHKERRQ(ierr); 347*e3cd5995SDave May ierr = DMGetCoordinateSection(dm,&coordSection);CHKERRQ(ierr); 348*e3cd5995SDave May 349*e3cd5995SDave May ierr = DMSwarmGetLocalSize(swarm,&npoints);CHKERRQ(ierr); 350*e3cd5995SDave May ierr = DMSwarmGetField(swarm,DMSwarmPICField_coor,NULL,NULL,(void**)&mpfield_coor);CHKERRQ(ierr); 351*e3cd5995SDave May ierr = DMSwarmGetField(swarm,DMSwarmPICField_cellid,NULL,NULL,(void**)&mpfield_cell);CHKERRQ(ierr); 352*e3cd5995SDave May 353*e3cd5995SDave May for (p=0; p<npoints; p++) { 354*e3cd5995SDave May PetscReal *coor_p; 355*e3cd5995SDave May PetscReal elfield[3],dJ; 356*e3cd5995SDave May PetscBool point_located; 357*e3cd5995SDave May 358*e3cd5995SDave May e = mpfield_cell[p]; 359*e3cd5995SDave May coor_p = &mpfield_coor[2*p]; 360*e3cd5995SDave May 361*e3cd5995SDave May ierr = DMPlexVecGetClosure(dm,coordSection,coor_l,e,NULL,&elcoor);CHKERRQ(ierr); 362*e3cd5995SDave May 363*e3cd5995SDave May /* 364*e3cd5995SDave May while (!point_located && (failed_counter < 25)) { 365*e3cd5995SDave May ierr = PointInTriangle(point, coords[0], coords[1], coords[2], &point_located);CHKERRQ(ierr); 366*e3cd5995SDave May point.x = coor_p[0]; 367*e3cd5995SDave May point.y = coor_p[1]; 368*e3cd5995SDave May point.x += 1.0e-10 * (2.0 * rand()/((double)RAND_MAX)-1.0); 369*e3cd5995SDave May point.y += 1.0e-10 * (2.0 * rand()/((double)RAND_MAX)-1.0); 370*e3cd5995SDave May failed_counter++; 371*e3cd5995SDave May } 372*e3cd5995SDave May 373*e3cd5995SDave May if (!point_located){ 374*e3cd5995SDave May 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*e3cd5995SDave May } 376*e3cd5995SDave May 377*e3cd5995SDave May 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*e3cd5995SDave May else { 379*e3cd5995SDave May ierr = _ComputeLocalCoordinateAffine2d(coor_p,elcoor,xi_p,&dJ);CHKERRQ(ierr); 380*e3cd5995SDave May // rescale to [0,1] // 381*e3cd5995SDave May xi_p[0] = 0.5*(xi_p[0] + 1.0); 382*e3cd5995SDave May xi_p[1] = 0.5*(xi_p[1] + 1.0); 383*e3cd5995SDave May 384*e3cd5995SDave May 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*e3cd5995SDave May 386*e3cd5995SDave May } 387*e3cd5995SDave May */ 388*e3cd5995SDave May 389*e3cd5995SDave May ierr = ComputeLocalCoordinateAffine2d(coor_p,elcoor,xi_p,&dJ);CHKERRQ(ierr); 390*e3cd5995SDave May /* 391*e3cd5995SDave May 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*e3cd5995SDave May */ 393*e3cd5995SDave May /* 394*e3cd5995SDave May point_located = PETSC_TRUE; 395*e3cd5995SDave May if (xi_p[0] < 0.0) { 396*e3cd5995SDave May if (xi_p[0] > -PLEX_C_EPS) { 397*e3cd5995SDave May xi_p[0] = 0.0; 398*e3cd5995SDave May } else { 399*e3cd5995SDave May point_located = PETSC_FALSE; 400*e3cd5995SDave May } 401*e3cd5995SDave May } 402*e3cd5995SDave May if (xi_p[1] < 0.0) { 403*e3cd5995SDave May if (xi_p[1] > -PLEX_C_EPS) { 404*e3cd5995SDave May xi_p[1] = 0.0; 405*e3cd5995SDave May } else { 406*e3cd5995SDave May point_located = PETSC_FALSE; 407*e3cd5995SDave May } 408*e3cd5995SDave May } 409*e3cd5995SDave May if (xi_p[1] > (1.0-xi_p[0])) { 410*e3cd5995SDave May if ((xi_p[1] - 1.0 + xi_p[0]) < PLEX_C_EPS) { 411*e3cd5995SDave May xi_p[1] = 1.0 - xi_p[0]; 412*e3cd5995SDave May } else { 413*e3cd5995SDave May point_located = PETSC_FALSE; 414*e3cd5995SDave May } 415*e3cd5995SDave May } 416*e3cd5995SDave May if (!point_located){ 417*e3cd5995SDave May PetscPrintf(PETSC_COMM_SELF,"[Error] xi,eta = %+1.8e, %+1.8e\n",xi_p[0],xi_p[1]); 418*e3cd5995SDave May 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*e3cd5995SDave May } 420*e3cd5995SDave May 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*e3cd5995SDave May */ 422*e3cd5995SDave May 423*e3cd5995SDave May Ni[0] = 1.0 - xi_p[0] - xi_p[1]; 424*e3cd5995SDave May Ni[1] = xi_p[0]; 425*e3cd5995SDave May Ni[2] = xi_p[1]; 426*e3cd5995SDave May 427*e3cd5995SDave May point_located = PETSC_TRUE; 428*e3cd5995SDave May for (k=0; k<3; k++) { 429*e3cd5995SDave May if (Ni[k] < -PLEX_C_EPS) point_located = PETSC_FALSE; 430*e3cd5995SDave May if (Ni[k] > (1.0+PLEX_C_EPS)) point_located = PETSC_FALSE; 431*e3cd5995SDave May } 432*e3cd5995SDave May if (!point_located){ 433*e3cd5995SDave May PetscPrintf(PETSC_COMM_SELF,"[Error] xi,eta = %+1.8e, %+1.8e\n",xi_p[0],xi_p[1]); 434*e3cd5995SDave May 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*e3cd5995SDave May } 436*e3cd5995SDave May 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*e3cd5995SDave May 438*e3cd5995SDave May for (k=0; k<3; k++) { 439*e3cd5995SDave May Ni[k] = Ni[k] * dJ; 440*e3cd5995SDave May elfield[k] = Ni[k] * swarm_field[p]; 441*e3cd5995SDave May } 442*e3cd5995SDave May 443*e3cd5995SDave May ierr = DMPlexVecRestoreClosure(dm,coordSection,coor_l,e,NULL,&elcoor);CHKERRQ(ierr); 444*e3cd5995SDave May 445*e3cd5995SDave May ierr = DMPlexVecSetClosure(dm, NULL,v_field_l, e, elfield, ADD_VALUES);CHKERRQ(ierr); 446*e3cd5995SDave May ierr = DMPlexVecSetClosure(dm, NULL,denom_l, e, Ni, ADD_VALUES);CHKERRQ(ierr); 447*e3cd5995SDave May } 448*e3cd5995SDave May 449*e3cd5995SDave May ierr = DMSwarmRestoreField(swarm,DMSwarmPICField_cellid,NULL,NULL,(void**)&mpfield_cell);CHKERRQ(ierr); 450*e3cd5995SDave May ierr = DMSwarmRestoreField(swarm,DMSwarmPICField_coor,NULL,NULL,(void**)&mpfield_coor);CHKERRQ(ierr); 451*e3cd5995SDave May 452*e3cd5995SDave May ierr = DMLocalToGlobalBegin(dm,v_field_l,ADD_VALUES,v_field);CHKERRQ(ierr); 453*e3cd5995SDave May ierr = DMLocalToGlobalEnd(dm,v_field_l,ADD_VALUES,v_field);CHKERRQ(ierr); 454*e3cd5995SDave May ierr = DMLocalToGlobalBegin(dm,denom_l,ADD_VALUES,denom);CHKERRQ(ierr); 455*e3cd5995SDave May ierr = DMLocalToGlobalEnd(dm,denom_l,ADD_VALUES,denom);CHKERRQ(ierr); 456*e3cd5995SDave May 457*e3cd5995SDave May ierr = VecPointwiseDivide(v_field,v_field,denom);CHKERRQ(ierr); 458*e3cd5995SDave May 459*e3cd5995SDave May ierr = DMRestoreLocalVector(dm,&v_field_l);CHKERRQ(ierr); 460*e3cd5995SDave May ierr = DMRestoreLocalVector(dm,&denom_l);CHKERRQ(ierr); 461*e3cd5995SDave May ierr = DMRestoreGlobalVector(dm,&denom);CHKERRQ(ierr); 462*e3cd5995SDave May 463*e3cd5995SDave May PetscFunctionReturn(0); 464*e3cd5995SDave May } 465*e3cd5995SDave May 466*e3cd5995SDave May PetscErrorCode private_DMSwarmProjectFields_PLEX(DM swarm,DM celldm,PetscInt project_type,PetscInt nfields,DataField dfield[],Vec vecs[]) 467*e3cd5995SDave May { 468*e3cd5995SDave May PetscErrorCode ierr; 469*e3cd5995SDave May PetscInt f,dim; 470*e3cd5995SDave May 471*e3cd5995SDave May PetscFunctionBegin; 472*e3cd5995SDave May ierr = DMGetDimension(swarm,&dim);CHKERRQ(ierr); 473*e3cd5995SDave May switch (dim) { 474*e3cd5995SDave May case 2: 475*e3cd5995SDave May for (f=0; f<nfields; f++) { 476*e3cd5995SDave May PetscReal *swarm_field; 477*e3cd5995SDave May 478*e3cd5995SDave May ierr = DataFieldGetEntries(dfield[f],(void**)&swarm_field);CHKERRQ(ierr); 479*e3cd5995SDave May ierr = DMSwarmProjectField_ApproxP1_PLEX_2D(swarm,swarm_field,celldm,vecs[f]);CHKERRQ(ierr); 480*e3cd5995SDave May } 481*e3cd5995SDave May break; 482*e3cd5995SDave May case 3: 483*e3cd5995SDave May SETERRQ(PetscObjectComm((PetscObject)swarm),PETSC_ERR_SUP,"No support for 3D"); 484*e3cd5995SDave May break; 485*e3cd5995SDave May default: 486*e3cd5995SDave May break; 487*e3cd5995SDave May } 488*e3cd5995SDave May 489*e3cd5995SDave May PetscFunctionReturn(0); 490*e3cd5995SDave May } 491