xref: /petsc/src/dm/impls/swarm/swarmpic_plex.c (revision e3cd5995fa23de2668d243c125fd8a043e34f9ce)
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