Actual source code: swarmpic_plex.c

      1: #include <petscdm.h>
      2: #include <petscdmplex.h>
      3: #include <petscdmswarm.h>
      4: #include "../src/dm/impls/swarm/data_bucket.h"

      7: PetscErrorCode private_DMSwarmSetPointCoordinatesCellwise_PLEX(DM,DM,PetscInt,PetscReal*xi);

      9: static PetscErrorCode private_PetscFECreateDefault_scalar_pk1(DM dm, PetscInt dim, PetscBool isSimplex, PetscInt qorder, PetscFE *fem)
     10: {
     11:   const PetscInt Nc = 1;
     12:   PetscQuadrature q, fq;
     13:   DM              K;
     14:   PetscSpace      P;
     15:   PetscDualSpace  Q;
     16:   PetscInt        order, quadPointsPerEdge;
     17:   PetscBool       tensor = isSimplex ? PETSC_FALSE : PETSC_TRUE;
     18:   PetscErrorCode  ierr;

     21:   /* Create space */
     22:   PetscSpaceCreate(PetscObjectComm((PetscObject) dm), &P);
     23:   /*PetscObjectSetOptionsPrefix((PetscObject) P, prefix);*/
     24:   PetscSpacePolynomialSetTensor(P, tensor);
     25:   /*PetscSpaceSetFromOptions(P);*/
     26:   PetscSpaceSetType(P,PETSCSPACEPOLYNOMIAL);
     27:   PetscSpaceSetDegree(P,1,PETSC_DETERMINE);
     28:   PetscSpaceSetNumComponents(P, Nc);
     29:   PetscSpaceSetNumVariables(P, dim);
     30:   PetscSpaceSetUp(P);
     31:   PetscSpaceGetDegree(P, &order, NULL);
     32:   PetscSpacePolynomialGetTensor(P, &tensor);
     33:   /* Create dual space */
     34:   PetscDualSpaceCreate(PetscObjectComm((PetscObject) dm), &Q);
     35:   PetscDualSpaceSetType(Q,PETSCDUALSPACELAGRANGE);
     36:   /*PetscObjectSetOptionsPrefix((PetscObject) Q, prefix);*/
     37:   PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K);
     38:   PetscDualSpaceSetDM(Q, K);
     39:   DMDestroy(&K);
     40:   PetscDualSpaceSetNumComponents(Q, Nc);
     41:   PetscDualSpaceSetOrder(Q, order);
     42:   PetscDualSpaceLagrangeSetTensor(Q, tensor);
     43:   /*PetscDualSpaceSetFromOptions(Q);*/
     44:   PetscDualSpaceSetType(Q,PETSCDUALSPACELAGRANGE);
     45:   PetscDualSpaceSetUp(Q);
     46:   /* Create element */
     47:   PetscFECreate(PetscObjectComm((PetscObject) dm), fem);
     48:   /*PetscObjectSetOptionsPrefix((PetscObject) *fem, prefix);*/
     49:   /*PetscFESetFromOptions(*fem);*/
     50:   PetscFESetType(*fem,PETSCFEBASIC);
     51:   PetscFESetBasisSpace(*fem, P);
     52:   PetscFESetDualSpace(*fem, Q);
     53:   PetscFESetNumComponents(*fem, Nc);
     54:   PetscFESetUp(*fem);
     55:   PetscSpaceDestroy(&P);
     56:   PetscDualSpaceDestroy(&Q);
     57:   /* Create quadrature (with specified order if given) */
     58:   qorder = qorder >= 0 ? qorder : order;
     59:   quadPointsPerEdge = PetscMax(qorder + 1,1);
     60:   if (isSimplex) {
     61:     PetscDTStroudConicalQuadrature(dim,   1, quadPointsPerEdge, -1.0, 1.0, &q);
     62:     PetscDTStroudConicalQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);
     63:   }
     64:   else {
     65:     PetscDTGaussTensorQuadrature(dim,   1, quadPointsPerEdge, -1.0, 1.0, &q);
     66:     PetscDTGaussTensorQuadrature(dim-1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);
     67:   }
     68:   PetscFESetQuadrature(*fem, q);
     69:   PetscFESetFaceQuadrature(*fem, fq);
     70:   PetscQuadratureDestroy(&q);
     71:   PetscQuadratureDestroy(&fq);
     72:   return(0);
     73: }

     75: PetscErrorCode subdivide_triangle(PetscReal v1[2],PetscReal v2[2],PetscReal v3[2],PetscInt depth,PetscInt max,PetscReal xi[],PetscInt *np)
     76: {
     77:   PetscReal v12[2],v23[2],v31[2];
     78:   PetscInt i;

     82:   if (depth == max) {
     83:     PetscReal cx[2];

     85:     cx[0] = (v1[0] + v2[0] + v3[0])/3.0;
     86:     cx[1] = (v1[1] + v2[1] + v3[1])/3.0;

     88:     xi[2*(*np)+0] = cx[0];
     89:     xi[2*(*np)+1] = cx[1];
     90:     (*np)++;
     91:     return(0);
     92:   }

     94:   /* calculate midpoints of each side */
     95:   for (i = 0; i < 2; i++) {
     96:     v12[i] = (v1[i]+v2[i])/2.0;
     97:     v23[i] = (v2[i]+v3[i])/2.0;
     98:     v31[i] = (v3[i]+v1[i])/2.0;
     99:   }

    101:   /* recursively subdivide new triangles */
    102:   subdivide_triangle(v1,v12,v31,depth+1,max,xi,np);
    103:   subdivide_triangle(v2,v23,v12,depth+1,max,xi,np);
    104:   subdivide_triangle(v3,v31,v23,depth+1,max,xi,np);
    105:   subdivide_triangle(v12,v23,v31,depth+1,max,xi,np);
    106:   return(0);
    107: }

    110: PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX2D_SubDivide(DM dm,DM dmc,PetscInt nsub)
    111: {
    113:   const PetscInt dim = 2;
    114:   PetscInt q,npoints_q,e,nel,npe,pcnt,ps,pe,d,k,depth;
    115:   PetscReal *xi;
    116:   PetscReal **basis;
    117:   Vec coorlocal;
    118:   PetscSection coordSection;
    119:   PetscScalar *elcoor = NULL;
    120:   PetscReal *swarm_coor;
    121:   PetscInt *swarm_cellid;
    122:   PetscReal v1[2],v2[2],v3[2];

    126:   npoints_q = 1;
    127:   for (d=0; d<nsub; d++) { npoints_q *= 4; }
    128:   PetscMalloc1(dim*npoints_q,&xi);

    130:   v1[0] = 0.0;  v1[1] = 0.0;
    131:   v2[0] = 1.0;  v2[1] = 0.0;
    132:   v3[0] = 0.0;  v3[1] = 1.0;
    133:   depth = 0;
    134:   pcnt = 0;
    135:   subdivide_triangle(v1,v2,v3,depth,nsub,xi,&pcnt);

    137:   npe = 3; /* nodes per element (triangle) */
    138:   PetscMalloc1(npoints_q,&basis);
    139:   for (q=0; q<npoints_q; q++) {
    140:     PetscMalloc1(npe,&basis[q]);

    142:     basis[q][0] = 1.0 - xi[dim*q+0] - xi[dim*q+1];
    143:     basis[q][1] = xi[dim*q+0];
    144:     basis[q][2] = xi[dim*q+1];
    145:   }

    147:   /* 0->cell, 1->edge, 2->vert */
    148:   DMPlexGetHeightStratum(dmc,0,&ps,&pe);
    149:   nel = pe - ps;

    151:   DMSwarmSetLocalSizes(dm,npoints_q*nel,-1);
    152:   DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
    153:   DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);

    155:   DMGetCoordinatesLocal(dmc,&coorlocal);
    156:   DMGetCoordinateSection(dmc,&coordSection);

    158:   pcnt = 0;
    159:   for (e=0; e<nel; e++) {
    160:     DMPlexVecGetClosure(dmc,coordSection,coorlocal,e,NULL,&elcoor);

    162:     for (q=0; q<npoints_q; q++) {
    163:       for (d=0; d<dim; d++) {
    164:         swarm_coor[dim*pcnt+d] = 0.0;
    165:         for (k=0; k<npe; k++) {
    166:           swarm_coor[dim*pcnt+d] += basis[q][k] * PetscRealPart(elcoor[dim*k+d]);
    167:         }
    168:       }
    169:       swarm_cellid[pcnt] = e;
    170:       pcnt++;
    171:     }
    172:     DMPlexVecRestoreClosure(dmc,coordSection,coorlocal,e,NULL,&elcoor);
    173:   }
    174:   DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
    175:   DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);

    177:   PetscFree(xi);
    178:   for (q=0; q<npoints_q; q++) {
    179:     PetscFree(basis[q]);
    180:   }
    181:   PetscFree(basis);

    183:   return(0);
    184: }

    186: PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX_SubDivide(DM dm,DM dmc,PetscInt nsub)
    187: {
    189:   PetscInt dim,nfaces,nbasis;
    190:   PetscInt q,npoints_q,e,nel,pcnt,ps,pe,d,k,r;
    191:   PetscTabulation T;
    192:   Vec coorlocal;
    193:   PetscSection coordSection;
    194:   PetscScalar *elcoor = NULL;
    195:   PetscReal *swarm_coor;
    196:   PetscInt *swarm_cellid;
    197:   const PetscReal *xiq;
    198:   PetscQuadrature quadrature;
    199:   PetscFE fe,feRef;
    200:   PetscBool is_simplex;

    204:   DMGetDimension(dmc,&dim);

    206:   is_simplex = PETSC_FALSE;
    207:   DMPlexGetHeightStratum(dmc,0,&ps,&pe);
    208:   DMPlexGetConeSize(dmc, ps, &nfaces);
    209:   if (nfaces == (dim+1)) { is_simplex = PETSC_TRUE; }

    211:   private_PetscFECreateDefault_scalar_pk1(dmc, dim, is_simplex, 0, &fe);

    213:   for (r=0; r<nsub; r++) {
    214:     PetscFERefine(fe,&feRef);
    215:     PetscFECopyQuadrature(feRef,fe);
    216:     PetscFEDestroy(&feRef);
    217:   }

    219:   PetscFEGetQuadrature(fe,&quadrature);
    220:   PetscQuadratureGetData(quadrature, NULL, NULL, &npoints_q, &xiq, NULL);
    221:   PetscFEGetDimension(fe,&nbasis);
    222:   PetscFEGetCellTabulation(fe, 1, &T);

    224:   /* 0->cell, 1->edge, 2->vert */
    225:   DMPlexGetHeightStratum(dmc,0,&ps,&pe);
    226:   nel = pe - ps;

    228:   DMSwarmSetLocalSizes(dm,npoints_q*nel,-1);
    229:   DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
    230:   DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);

    232:   DMGetCoordinatesLocal(dmc,&coorlocal);
    233:   DMGetCoordinateSection(dmc,&coordSection);

    235:   pcnt = 0;
    236:   for (e=0; e<nel; e++) {
    237:     DMPlexVecGetClosure(dmc,coordSection,coorlocal,ps+e,NULL,&elcoor);

    239:     for (q=0; q<npoints_q; q++) {
    240:       for (d=0; d<dim; d++) {
    241:         swarm_coor[dim*pcnt+d] = 0.0;
    242:         for (k=0; k<nbasis; k++) {
    243:           swarm_coor[dim*pcnt+d] += T->T[0][q*nbasis + k] * PetscRealPart(elcoor[dim*k+d]);
    244:         }
    245:       }
    246:       swarm_cellid[pcnt] = e;
    247:       pcnt++;
    248:     }
    249:     DMPlexVecRestoreClosure(dmc,coordSection,coorlocal,ps+e,NULL,&elcoor);
    250:   }
    251:   DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
    252:   DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);

    254:   PetscFEDestroy(&fe);
    255:   return(0);
    256: }

    258: PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX2D_Regular(DM dm,DM dmc,PetscInt npoints)
    259: {
    261:   PetscInt dim;
    262:   PetscInt ii,jj,q,npoints_q,e,nel,npe,pcnt,ps,pe,d,k,nfaces;
    263:   PetscReal *xi,ds,ds2;
    264:   PetscReal **basis;
    265:   Vec coorlocal;
    266:   PetscSection coordSection;
    267:   PetscScalar *elcoor = NULL;
    268:   PetscReal *swarm_coor;
    269:   PetscInt *swarm_cellid;
    270:   PetscBool is_simplex;

    273:   DMGetDimension(dmc,&dim);
    274:   if (dim != 2) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only 2D is supported");
    275:   is_simplex = PETSC_FALSE;
    276:   DMPlexGetHeightStratum(dmc,0,&ps,&pe);
    277:   DMPlexGetConeSize(dmc, ps, &nfaces);
    278:   if (nfaces == (dim+1)) { is_simplex = PETSC_TRUE; }
    279:   if (!is_simplex) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only the simplex is supported");

    281:   PetscMalloc1(dim*npoints*npoints,&xi);
    282:   pcnt = 0;
    283:   ds = 1.0/((PetscReal)(npoints-1));
    284:   ds2 = 1.0/((PetscReal)(npoints));
    285:   for (jj = 0; jj<npoints; jj++) {
    286:     for (ii=0; ii<npoints-jj; ii++) {
    287:       xi[dim*pcnt+0] = ii * ds;
    288:       xi[dim*pcnt+1] = jj * ds;

    290:       xi[dim*pcnt+0] *= (1.0 - 1.2*ds2);
    291:       xi[dim*pcnt+1] *= (1.0 - 1.2*ds2);

    293:       xi[dim*pcnt+0] += 0.35*ds2;
    294:       xi[dim*pcnt+1] += 0.35*ds2;

    296:       pcnt++;
    297:     }
    298:   }
    299:   npoints_q = pcnt;

    301:   npe = 3; /* nodes per element (triangle) */
    302:   PetscMalloc1(npoints_q,&basis);
    303:   for (q=0; q<npoints_q; q++) {
    304:     PetscMalloc1(npe,&basis[q]);

    306:     basis[q][0] = 1.0 - xi[dim*q+0] - xi[dim*q+1];
    307:     basis[q][1] = xi[dim*q+0];
    308:     basis[q][2] = xi[dim*q+1];
    309:   }

    311:   /* 0->cell, 1->edge, 2->vert */
    312:   DMPlexGetHeightStratum(dmc,0,&ps,&pe);
    313:   nel = pe - ps;

    315:   DMSwarmSetLocalSizes(dm,npoints_q*nel,-1);
    316:   DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
    317:   DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);

    319:   DMGetCoordinatesLocal(dmc,&coorlocal);
    320:   DMGetCoordinateSection(dmc,&coordSection);

    322:   pcnt = 0;
    323:   for (e=0; e<nel; e++) {
    324:     DMPlexVecGetClosure(dmc,coordSection,coorlocal,e,NULL,&elcoor);

    326:     for (q=0; q<npoints_q; q++) {
    327:       for (d=0; d<dim; d++) {
    328:         swarm_coor[dim*pcnt+d] = 0.0;
    329:         for (k=0; k<npe; k++) {
    330:           swarm_coor[dim*pcnt+d] += basis[q][k] * PetscRealPart(elcoor[dim*k+d]);
    331:         }
    332:       }
    333:       swarm_cellid[pcnt] = e;
    334:       pcnt++;
    335:     }
    336:     DMPlexVecRestoreClosure(dmc,coordSection,coorlocal,e,NULL,&elcoor);
    337:   }
    338:   DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
    339:   DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);

    341:   PetscFree(xi);
    342:   for (q=0; q<npoints_q; q++) {
    343:     PetscFree(basis[q]);
    344:   }
    345:   PetscFree(basis);

    347:   return(0);
    348: }

    350: PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX(DM dm,DM celldm,DMSwarmPICLayoutType layout,PetscInt layout_param)
    351: {
    353:   PetscInt dim;

    356:   DMGetDimension(celldm,&dim);
    357:   switch (layout) {
    358:     case DMSWARMPIC_LAYOUT_REGULAR:
    359:       if (dim == 3) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No 3D support for REGULAR+PLEX");
    360:       private_DMSwarmInsertPointsUsingCellDM_PLEX2D_Regular(dm,celldm,layout_param);
    361:       break;
    362:     case DMSWARMPIC_LAYOUT_GAUSS:
    363:     {
    364:       PetscInt npoints,npoints1,ps,pe,nfaces;
    365:       const PetscReal *xi;
    366:       PetscBool is_simplex;
    367:       PetscQuadrature quadrature;

    369:       is_simplex = PETSC_FALSE;
    370:       DMPlexGetHeightStratum(celldm,0,&ps,&pe);
    371:       DMPlexGetConeSize(celldm, ps, &nfaces);
    372:       if (nfaces == (dim+1)) { is_simplex = PETSC_TRUE; }

    374:       npoints1 = layout_param;
    375:       if (is_simplex) {
    376:         PetscDTStroudConicalQuadrature(dim,1,npoints1,-1.0,1.0,&quadrature);
    377:       } else {
    378:         PetscDTGaussTensorQuadrature(dim,1,npoints1,-1.0,1.0,&quadrature);
    379:       }
    380:       PetscQuadratureGetData(quadrature,NULL,NULL,&npoints,&xi,NULL);
    381:       private_DMSwarmSetPointCoordinatesCellwise_PLEX(dm,celldm,npoints,(PetscReal*)xi);
    382:       PetscQuadratureDestroy(&quadrature);
    383:     }
    384:       break;
    385:     case DMSWARMPIC_LAYOUT_SUBDIVISION:
    386:       private_DMSwarmInsertPointsUsingCellDM_PLEX_SubDivide(dm,celldm,layout_param);
    387:       break;
    388:   }
    389:   return(0);
    390: }

    392: /*
    393: typedef struct {
    394:   PetscReal x,y;
    395: } Point2d;

    397: static PetscErrorCode signp2d(Point2d p1, Point2d p2, Point2d p3,PetscReal *s)
    398: {
    400:   *s = (p1.x - p3.x) * (p2.y - p3.y) - (p2.x - p3.x) * (p1.y - p3.y);
    401:   return(0);
    402: }
    403: */
    404: /*
    405: static PetscErrorCode PointInTriangle(Point2d pt, Point2d v1, Point2d v2, Point2d v3,PetscBool *v)
    406: {
    407:   PetscReal s1,s2,s3;
    408:   PetscBool b1, b2, b3;

    411:   signp2d(pt, v1, v2,&s1); b1 = s1 < 0.0f;
    412:   signp2d(pt, v2, v3,&s2); b2 = s2 < 0.0f;
    413:   signp2d(pt, v3, v1,&s3); b3 = s3 < 0.0f;

    415:   *v = ((b1 == b2) && (b2 == b3));
    416:   return(0);
    417: }
    418: */
    419: /*
    420: static PetscErrorCode _ComputeLocalCoordinateAffine2d(PetscReal xp[],PetscReal coords[],PetscReal xip[],PetscReal *dJ)
    421: {
    422:   PetscReal x1,y1,x2,y2,x3,y3;
    423:   PetscReal c,b[2],A[2][2],inv[2][2],detJ,od;

    426:   x1 = coords[2*0+0];
    427:   x2 = coords[2*1+0];
    428:   x3 = coords[2*2+0];

    430:   y1 = coords[2*0+1];
    431:   y2 = coords[2*1+1];
    432:   y3 = coords[2*2+1];

    434:   c = x1 - 0.5*x1 - 0.5*x1 + 0.5*x2 + 0.5*x3;
    435:   b[0] = xp[0] - c;
    436:   c = y1 - 0.5*y1 - 0.5*y1 + 0.5*y2 + 0.5*y3;
    437:   b[1] = xp[1] - c;

    439:   A[0][0] = -0.5*x1 + 0.5*x2;   A[0][1] = -0.5*x1 + 0.5*x3;
    440:   A[1][0] = -0.5*y1 + 0.5*y2;   A[1][1] = -0.5*y1 + 0.5*y3;

    442:   detJ = A[0][0]*A[1][1] - A[0][1]*A[1][0];
    443:   *dJ = PetscAbsReal(detJ);
    444:   od = 1.0/detJ;

    446:   inv[0][0] =  A[1][1] * od;
    447:   inv[0][1] = -A[0][1] * od;
    448:   inv[1][0] = -A[1][0] * od;
    449:   inv[1][1] =  A[0][0] * od;

    451:   xip[0] = inv[0][0]*b[0] + inv[0][1]*b[1];
    452:   xip[1] = inv[1][0]*b[0] + inv[1][1]*b[1];
    453:   return(0);
    454: }
    455: */

    457: static PetscErrorCode ComputeLocalCoordinateAffine2d(PetscReal xp[],PetscScalar coords[],PetscReal xip[],PetscReal *dJ)
    458: {
    459:   PetscReal x1,y1,x2,y2,x3,y3;
    460:   PetscReal b[2],A[2][2],inv[2][2],detJ,od;

    463:   x1 = PetscRealPart(coords[2*0+0]);
    464:   x2 = PetscRealPart(coords[2*1+0]);
    465:   x3 = PetscRealPart(coords[2*2+0]);

    467:   y1 = PetscRealPart(coords[2*0+1]);
    468:   y2 = PetscRealPart(coords[2*1+1]);
    469:   y3 = PetscRealPart(coords[2*2+1]);

    471:   b[0] = xp[0] - x1;
    472:   b[1] = xp[1] - y1;

    474:   A[0][0] = x2-x1;   A[0][1] = x3-x1;
    475:   A[1][0] = y2-y1;   A[1][1] = y3-y1;

    477:   detJ = A[0][0]*A[1][1] - A[0][1]*A[1][0];
    478:   *dJ = PetscAbsReal(detJ);
    479:   od = 1.0/detJ;

    481:   inv[0][0] =  A[1][1] * od;
    482:   inv[0][1] = -A[0][1] * od;
    483:   inv[1][0] = -A[1][0] * od;
    484:   inv[1][1] =  A[0][0] * od;

    486:   xip[0] = inv[0][0]*b[0] + inv[0][1]*b[1];
    487:   xip[1] = inv[1][0]*b[0] + inv[1][1]*b[1];
    488:   return(0);
    489: }

    491: PetscErrorCode DMSwarmProjectField_ApproxP1_PLEX_2D(DM swarm,PetscReal *swarm_field,DM dm,Vec v_field)
    492: {
    494:   const PetscReal PLEX_C_EPS = 1.0e-8;
    495:   Vec v_field_l,denom_l,coor_l,denom;
    496:   PetscInt k,p,e,npoints;
    497:   PetscInt *mpfield_cell;
    498:   PetscReal *mpfield_coor;
    499:   PetscReal xi_p[2];
    500:   PetscScalar Ni[3];
    501:   PetscSection coordSection;
    502:   PetscScalar *elcoor = NULL;

    505:   VecZeroEntries(v_field);

    507:   DMGetLocalVector(dm,&v_field_l);
    508:   DMGetGlobalVector(dm,&denom);
    509:   DMGetLocalVector(dm,&denom_l);
    510:   VecZeroEntries(v_field_l);
    511:   VecZeroEntries(denom);
    512:   VecZeroEntries(denom_l);

    514:   DMGetCoordinatesLocal(dm,&coor_l);
    515:   DMGetCoordinateSection(dm,&coordSection);

    517:   DMSwarmGetLocalSize(swarm,&npoints);
    518:   DMSwarmGetField(swarm,DMSwarmPICField_coor,NULL,NULL,(void**)&mpfield_coor);
    519:   DMSwarmGetField(swarm,DMSwarmPICField_cellid,NULL,NULL,(void**)&mpfield_cell);

    521:   for (p=0; p<npoints; p++) {
    522:     PetscReal   *coor_p,dJ;
    523:     PetscScalar elfield[3];
    524:     PetscBool   point_located;

    526:     e       = mpfield_cell[p];
    527:     coor_p  = &mpfield_coor[2*p];

    529:     DMPlexVecGetClosure(dm,coordSection,coor_l,e,NULL,&elcoor);

    531: /*
    532:     while (!point_located && (failed_counter < 25)) {
    533:       PointInTriangle(point, coords[0], coords[1], coords[2], &point_located);
    534:       point.x = coor_p[0];
    535:       point.y = coor_p[1];
    536:       point.x += 1.0e-10 * (2.0 * rand()/((double)RAND_MAX)-1.0);
    537:       point.y += 1.0e-10 * (2.0 * rand()/((double)RAND_MAX)-1.0);
    538:       failed_counter++;
    539:     }

    541:     if (!point_located){
    542:         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);
    543:     }

    545:     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);
    546:     else {
    547:       _ComputeLocalCoordinateAffine2d(coor_p,elcoor,xi_p,&dJ);
    548:       xi_p[0] = 0.5*(xi_p[0] + 1.0);
    549:       xi_p[1] = 0.5*(xi_p[1] + 1.0);

    551:       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]);

    553:     }
    554: */

    556:     ComputeLocalCoordinateAffine2d(coor_p,elcoor,xi_p,&dJ);
    557:     /*
    558:     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]);
    559:     */
    560:     /*
    561:      point_located = PETSC_TRUE;
    562:     if (xi_p[0] < 0.0) {
    563:       if (xi_p[0] > -PLEX_C_EPS) {
    564:         xi_p[0] = 0.0;
    565:       } else {
    566:         point_located = PETSC_FALSE;
    567:       }
    568:     }
    569:     if (xi_p[1] < 0.0) {
    570:       if (xi_p[1] > -PLEX_C_EPS) {
    571:         xi_p[1] = 0.0;
    572:       } else {
    573:         point_located = PETSC_FALSE;
    574:       }
    575:     }
    576:     if (xi_p[1] > (1.0-xi_p[0])) {
    577:       if ((xi_p[1] - 1.0 + xi_p[0]) < PLEX_C_EPS) {
    578:         xi_p[1] = 1.0 - xi_p[0];
    579:       } else {
    580:         point_located = PETSC_FALSE;
    581:       }
    582:     }
    583:     if (!point_located){
    584:       PetscPrintf(PETSC_COMM_SELF,"[Error] xi,eta = %+1.8e, %+1.8e\n",xi_p[0],xi_p[1]);
    585:       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]);
    586:     }
    587:     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);
    588:     */

    590:     Ni[0] = 1.0 - xi_p[0] - xi_p[1];
    591:     Ni[1] = xi_p[0];
    592:     Ni[2] = xi_p[1];

    594:     point_located = PETSC_TRUE;
    595:     for (k=0; k<3; k++) {
    596:       if (PetscRealPart(Ni[k]) < -PLEX_C_EPS) point_located = PETSC_FALSE;
    597:       if (PetscRealPart(Ni[k]) > (1.0+PLEX_C_EPS)) point_located = PETSC_FALSE;
    598:     }
    599:     if (!point_located){
    600:       PetscPrintf(PETSC_COMM_SELF,"[Error] xi,eta = %+1.8e, %+1.8e\n",(double)xi_p[0],(double)xi_p[1]);
    601:       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]));
    602:     }
    603:     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);

    605:     for (k=0; k<3; k++) {
    606:       Ni[k] = Ni[k] * dJ;
    607:       elfield[k] = Ni[k] * swarm_field[p];
    608:     }

    610:     DMPlexVecRestoreClosure(dm,coordSection,coor_l,e,NULL,&elcoor);

    612:     DMPlexVecSetClosure(dm, NULL,v_field_l, e, elfield, ADD_VALUES);
    613:     DMPlexVecSetClosure(dm, NULL,denom_l, e, Ni, ADD_VALUES);
    614:   }

    616:   DMSwarmRestoreField(swarm,DMSwarmPICField_cellid,NULL,NULL,(void**)&mpfield_cell);
    617:   DMSwarmRestoreField(swarm,DMSwarmPICField_coor,NULL,NULL,(void**)&mpfield_coor);

    619:   DMLocalToGlobalBegin(dm,v_field_l,ADD_VALUES,v_field);
    620:   DMLocalToGlobalEnd(dm,v_field_l,ADD_VALUES,v_field);
    621:   DMLocalToGlobalBegin(dm,denom_l,ADD_VALUES,denom);
    622:   DMLocalToGlobalEnd(dm,denom_l,ADD_VALUES,denom);

    624:   VecPointwiseDivide(v_field,v_field,denom);

    626:   DMRestoreLocalVector(dm,&v_field_l);
    627:   DMRestoreLocalVector(dm,&denom_l);
    628:   DMRestoreGlobalVector(dm,&denom);

    630:   return(0);
    631: }

    633: PetscErrorCode private_DMSwarmProjectFields_PLEX(DM swarm,DM celldm,PetscInt project_type,PetscInt nfields,DMSwarmDataField dfield[],Vec vecs[])
    634: {
    636:   PetscInt f,dim;

    639:   DMGetDimension(swarm,&dim);
    640:   switch (dim) {
    641:     case 2:
    642:       for (f=0; f<nfields; f++) {
    643:         PetscReal *swarm_field;

    645:         DMSwarmDataFieldGetEntries(dfield[f],(void**)&swarm_field);
    646:         DMSwarmProjectField_ApproxP1_PLEX_2D(swarm,swarm_field,celldm,vecs[f]);
    647:       }
    648:       break;
    649:     case 3:
    650:       SETERRQ(PetscObjectComm((PetscObject)swarm),PETSC_ERR_SUP,"No support for 3D");
    651:     default:
    652:       break;
    653:   }

    655:   return(0);
    656: }

    658: PetscErrorCode private_DMSwarmSetPointCoordinatesCellwise_PLEX(DM dm,DM dmc,PetscInt npoints,PetscReal xi[])
    659: {
    660:   PetscBool is_simplex,is_tensorcell;
    662:   PetscInt dim,nfaces,ps,pe,p,d,nbasis,pcnt,e,k,nel;
    663:   PetscFE fe;
    664:   PetscQuadrature quadrature;
    665:   PetscTabulation T;
    666:   PetscReal *xiq;
    667:   Vec coorlocal;
    668:   PetscSection coordSection;
    669:   PetscScalar *elcoor = NULL;
    670:   PetscReal *swarm_coor;
    671:   PetscInt *swarm_cellid;

    674:   DMGetDimension(dmc,&dim);

    676:   is_simplex = PETSC_FALSE;
    677:   is_tensorcell = PETSC_FALSE;
    678:   DMPlexGetHeightStratum(dmc,0,&ps,&pe);
    679:   DMPlexGetConeSize(dmc, ps, &nfaces);

    681:   if (nfaces == (dim+1)) { is_simplex = PETSC_TRUE; }

    683:   switch (dim) {
    684:     case 2:
    685:       if (nfaces == 4) { is_tensorcell = PETSC_TRUE; }
    686:       break;
    687:     case 3:
    688:       if (nfaces == 6) { is_tensorcell = PETSC_TRUE; }
    689:       break;
    690:     default:
    691:       SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only support for 2D, 3D");
    692:   }

    694:   /* check points provided fail inside the reference cell */
    695:   if (is_simplex) {
    696:     for (p=0; p<npoints; p++) {
    697:       PetscReal sum;
    698:       for (d=0; d<dim; d++) {
    699:         if (xi[dim*p+d] < -1.0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Points do not fail inside the simplex domain");
    700:       }
    701:       sum = 0.0;
    702:       for (d=0; d<dim; d++) {
    703:         sum += xi[dim*p+d];
    704:       }
    705:       if (sum > 0.0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Points do not fail inside the simplex domain");
    706:     }
    707:   } else if (is_tensorcell) {
    708:     for (p=0; p<npoints; p++) {
    709:       for (d=0; d<dim; d++) {
    710:         if (PetscAbsReal(xi[dim*p+d]) > 1.0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Points do not fail inside the tensor domain [-1,1]^d");
    711:       }
    712:     }
    713:   } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only support for d-simplex and d-tensorcell");

    715:   PetscQuadratureCreate(PetscObjectComm((PetscObject)dm),&quadrature);
    716:   PetscMalloc1(npoints*dim,&xiq);
    717:   PetscArraycpy(xiq,xi,npoints*dim);
    718:   PetscQuadratureSetData(quadrature,dim,1,npoints,(const PetscReal*)xiq,NULL);
    719:   private_PetscFECreateDefault_scalar_pk1(dmc, dim, is_simplex, 0, &fe);
    720:   PetscFESetQuadrature(fe,quadrature);
    721:   PetscFEGetDimension(fe,&nbasis);
    722:   PetscFEGetCellTabulation(fe, 1, &T);

    724:   /* for each cell, interpolate coordaintes and insert the interpolated points coordinates into swarm */
    725:   /* 0->cell, 1->edge, 2->vert */
    726:   DMPlexGetHeightStratum(dmc,0,&ps,&pe);
    727:   nel = pe - ps;

    729:   DMSwarmSetLocalSizes(dm,npoints*nel,-1);
    730:   DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
    731:   DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);

    733:   DMGetCoordinatesLocal(dmc,&coorlocal);
    734:   DMGetCoordinateSection(dmc,&coordSection);

    736:   pcnt = 0;
    737:   for (e=0; e<nel; e++) {
    738:     DMPlexVecGetClosure(dmc,coordSection,coorlocal,ps+e,NULL,&elcoor);

    740:     for (p=0; p<npoints; p++) {
    741:       for (d=0; d<dim; d++) {
    742:         swarm_coor[dim*pcnt+d] = 0.0;
    743:         for (k=0; k<nbasis; k++) {
    744:           swarm_coor[dim*pcnt+d] += T->T[0][p*nbasis + k] * PetscRealPart(elcoor[dim*k+d]);
    745:         }
    746:       }
    747:       swarm_cellid[pcnt] = e;
    748:       pcnt++;
    749:     }
    750:     DMPlexVecRestoreClosure(dmc,coordSection,coorlocal,ps+e,NULL,&elcoor);
    751:   }
    752:   DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
    753:   DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);

    755:   PetscQuadratureDestroy(&quadrature);
    756:   PetscFEDestroy(&fe);

    758:   return(0);
    759: }