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"

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

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

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

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

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

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

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

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

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

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


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

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

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

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

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

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

153:   DMGetCoordinatesLocal(dmc,&coorlocal);
154:   DMGetCoordinateSection(dmc,&coordSection);

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

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

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

181:   return(0);
182: }

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


202:   DMGetDimension(dmc,&dim);

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

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

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

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

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

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

230:   DMGetCoordinatesLocal(dmc,&coorlocal);
231:   DMGetCoordinateSection(dmc,&coordSection);

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

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

252:   PetscFEDestroy(&fe);
253:   return(0);
254: }

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

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

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

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

291:       xi[dim*pcnt+0] += 0.35*ds2;
292:       xi[dim*pcnt+1] += 0.35*ds2;

294:       pcnt++;
295:     }
296:   }
297:   npoints_q = pcnt;

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

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

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

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

317:   DMGetCoordinatesLocal(dmc,&coorlocal);
318:   DMGetCoordinateSection(dmc,&coordSection);

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

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

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

345:   return(0);
346: }

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

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

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

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

390: /*
391: typedef struct {
392:   PetscReal x,y;
393: } Point2d;

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

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

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

424:   x1 = coords[2*0+0];
425:   x2 = coords[2*1+0];
426:   x3 = coords[2*2+0];

428:   y1 = coords[2*0+1];
429:   y2 = coords[2*1+1];
430:   y3 = coords[2*2+1];

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

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

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

444:   inv[0][0] =  A[1][1] * od;
445:   inv[0][1] = -A[0][1] * od;
446:   inv[1][0] = -A[1][0] * od;
447:   inv[1][1] =  A[0][0] * od;

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

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

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

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

469:   b[0] = xp[0] - x1;
470:   b[1] = xp[1] - y1;

472:   A[0][0] = x2-x1;   A[0][1] = x3-x1;
473:   A[1][0] = y2-y1;   A[1][1] = y3-y1;

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

479:   inv[0][0] =  A[1][1] * od;
480:   inv[0][1] = -A[0][1] * od;
481:   inv[1][0] = -A[1][0] * od;
482:   inv[1][1] =  A[0][0] * od;

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

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

503:   VecZeroEntries(v_field);

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

512:   DMGetCoordinatesLocal(dm,&coor_l);
513:   DMGetCoordinateSection(dm,&coordSection);

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

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

524:     e       = mpfield_cell[p];
525:     coor_p  = &mpfield_coor[2*p];

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

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

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

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

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

551:     }
552: */

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

588:     Ni[0] = 1.0 - xi_p[0] - xi_p[1];
589:     Ni[1] = xi_p[0];
590:     Ni[2] = xi_p[1];

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

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

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

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

614:   DMSwarmRestoreField(swarm,DMSwarmPICField_cellid,NULL,NULL,(void**)&mpfield_cell);
615:   DMSwarmRestoreField(swarm,DMSwarmPICField_coor,NULL,NULL,(void**)&mpfield_coor);

617:   DMLocalToGlobalBegin(dm,v_field_l,ADD_VALUES,v_field);
618:   DMLocalToGlobalEnd(dm,v_field_l,ADD_VALUES,v_field);
619:   DMLocalToGlobalBegin(dm,denom_l,ADD_VALUES,denom);
620:   DMLocalToGlobalEnd(dm,denom_l,ADD_VALUES,denom);

622:   VecPointwiseDivide(v_field,v_field,denom);

624:   DMRestoreLocalVector(dm,&v_field_l);
625:   DMRestoreLocalVector(dm,&denom_l);
626:   DMRestoreGlobalVector(dm,&denom);

628:   return(0);
629: }

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

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

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

653:   return(0);
654: }

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

672:   DMGetDimension(dmc,&dim);

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

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

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

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

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

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

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

731:   DMGetCoordinatesLocal(dmc,&coorlocal);
732:   DMGetCoordinateSection(dmc,&coordSection);

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

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

753:   PetscQuadratureDestroy(&quadrature);
754:   PetscFEDestroy(&fe);

756:   return(0);
757: }