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;

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

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

 77:   if (depth == max) {
 78:     PetscReal cx[2];

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

 83:     xi[2*(*np)+0] = cx[0];
 84:     xi[2*(*np)+1] = cx[1];
 85:     (*np)++;
 86:     return 0;
 87:   }

 89:   /* calculate midpoints of each side */
 90:   for (i = 0; i < 2; i++) {
 91:     v12[i] = (v1[i]+v2[i])/2.0;
 92:     v23[i] = (v2[i]+v3[i])/2.0;
 93:     v31[i] = (v3[i]+v1[i])/2.0;
 94:   }

 96:   /* recursively subdivide new triangles */
 97:   subdivide_triangle(v1,v12,v31,depth+1,max,xi,np);
 98:   subdivide_triangle(v2,v23,v12,depth+1,max,xi,np);
 99:   subdivide_triangle(v3,v31,v23,depth+1,max,xi,np);
100:   subdivide_triangle(v12,v23,v31,depth+1,max,xi,np);
101:   return 0;
102: }

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

117:   npoints_q = 1;
118:   for (d=0; d<nsub; d++) { npoints_q *= 4; }
119:   PetscMalloc1(dim*npoints_q,&xi);

121:   v1[0] = 0.0;  v1[1] = 0.0;
122:   v2[0] = 1.0;  v2[1] = 0.0;
123:   v3[0] = 0.0;  v3[1] = 1.0;
124:   depth = 0;
125:   pcnt = 0;
126:   subdivide_triangle(v1,v2,v3,depth,nsub,xi,&pcnt);

128:   npe = 3; /* nodes per element (triangle) */
129:   PetscMalloc1(npoints_q,&basis);
130:   for (q=0; q<npoints_q; q++) {
131:     PetscMalloc1(npe,&basis[q]);

133:     basis[q][0] = 1.0 - xi[dim*q+0] - xi[dim*q+1];
134:     basis[q][1] = xi[dim*q+0];
135:     basis[q][2] = xi[dim*q+1];
136:   }

138:   /* 0->cell, 1->edge, 2->vert */
139:   DMPlexGetHeightStratum(dmc,0,&ps,&pe);
140:   nel = pe - ps;

142:   DMSwarmSetLocalSizes(dm,npoints_q*nel,-1);
143:   DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
144:   DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);

146:   DMGetCoordinatesLocal(dmc,&coorlocal);
147:   DMGetCoordinateSection(dmc,&coordSection);

149:   pcnt = 0;
150:   for (e=0; e<nel; e++) {
151:     DMPlexVecGetClosure(dmc,coordSection,coorlocal,e,NULL,&elcoor);

153:     for (q=0; q<npoints_q; q++) {
154:       for (d=0; d<dim; d++) {
155:         swarm_coor[dim*pcnt+d] = 0.0;
156:         for (k=0; k<npe; k++) {
157:           swarm_coor[dim*pcnt+d] += basis[q][k] * PetscRealPart(elcoor[dim*k+d]);
158:         }
159:       }
160:       swarm_cellid[pcnt] = e;
161:       pcnt++;
162:     }
163:     DMPlexVecRestoreClosure(dmc,coordSection,coorlocal,e,NULL,&elcoor);
164:   }
165:   DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
166:   DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);

168:   PetscFree(xi);
169:   for (q=0; q<npoints_q; q++) {
170:     PetscFree(basis[q]);
171:   }
172:   PetscFree(basis);
173:   return 0;
174: }

176: PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX_SubDivide(DM dm,DM dmc,PetscInt nsub)
177: {
178:   PetscInt        dim,nfaces,nbasis;
179:   PetscInt        q,npoints_q,e,nel,pcnt,ps,pe,d,k,r;
180:   PetscTabulation T;
181:   Vec             coorlocal;
182:   PetscSection    coordSection;
183:   PetscScalar     *elcoor = NULL;
184:   PetscReal       *swarm_coor;
185:   PetscInt        *swarm_cellid;
186:   const PetscReal *xiq;
187:   PetscQuadrature quadrature;
188:   PetscFE         fe,feRef;
189:   PetscBool       is_simplex;

191:   DMGetDimension(dmc,&dim);
192:   is_simplex = PETSC_FALSE;
193:   DMPlexGetHeightStratum(dmc,0,&ps,&pe);
194:   DMPlexGetConeSize(dmc, ps, &nfaces);
195:   if (nfaces == (dim+1)) { is_simplex = PETSC_TRUE; }

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

199:   for (r=0; r<nsub; r++) {
200:     PetscFERefine(fe,&feRef);
201:     PetscFECopyQuadrature(feRef,fe);
202:     PetscFEDestroy(&feRef);
203:   }

205:   PetscFEGetQuadrature(fe,&quadrature);
206:   PetscQuadratureGetData(quadrature, NULL, NULL, &npoints_q, &xiq, NULL);
207:   PetscFEGetDimension(fe,&nbasis);
208:   PetscFEGetCellTabulation(fe, 1, &T);

210:   /* 0->cell, 1->edge, 2->vert */
211:   DMPlexGetHeightStratum(dmc,0,&ps,&pe);
212:   nel = pe - ps;

214:   DMSwarmSetLocalSizes(dm,npoints_q*nel,-1);
215:   DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
216:   DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);

218:   DMGetCoordinatesLocal(dmc,&coorlocal);
219:   DMGetCoordinateSection(dmc,&coordSection);

221:   pcnt = 0;
222:   for (e=0; e<nel; e++) {
223:     DMPlexVecGetClosure(dmc,coordSection,coorlocal,ps+e,NULL,&elcoor);

225:     for (q=0; q<npoints_q; q++) {
226:       for (d=0; d<dim; d++) {
227:         swarm_coor[dim*pcnt+d] = 0.0;
228:         for (k=0; k<nbasis; k++) {
229:           swarm_coor[dim*pcnt+d] += T->T[0][q*nbasis + k] * PetscRealPart(elcoor[dim*k+d]);
230:         }
231:       }
232:       swarm_cellid[pcnt] = e;
233:       pcnt++;
234:     }
235:     DMPlexVecRestoreClosure(dmc,coordSection,coorlocal,ps+e,NULL,&elcoor);
236:   }
237:   DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
238:   DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);

240:   PetscFEDestroy(&fe);
241:   return 0;
242: }

244: PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX2D_Regular(DM dm,DM dmc,PetscInt npoints)
245: {
246:   PetscInt       dim;
247:   PetscInt       ii,jj,q,npoints_q,e,nel,npe,pcnt,ps,pe,d,k,nfaces;
248:   PetscReal      *xi,ds,ds2;
249:   PetscReal      **basis;
250:   Vec            coorlocal;
251:   PetscSection   coordSection;
252:   PetscScalar    *elcoor = NULL;
253:   PetscReal      *swarm_coor;
254:   PetscInt       *swarm_cellid;
255:   PetscBool      is_simplex;

257:   DMGetDimension(dmc,&dim);
259:   is_simplex = PETSC_FALSE;
260:   DMPlexGetHeightStratum(dmc,0,&ps,&pe);
261:   DMPlexGetConeSize(dmc, ps, &nfaces);
262:   if (nfaces == (dim+1)) { is_simplex = PETSC_TRUE; }

265:   PetscMalloc1(dim*npoints*npoints,&xi);
266:   pcnt = 0;
267:   ds = 1.0/((PetscReal)(npoints-1));
268:   ds2 = 1.0/((PetscReal)(npoints));
269:   for (jj = 0; jj<npoints; jj++) {
270:     for (ii=0; ii<npoints-jj; ii++) {
271:       xi[dim*pcnt+0] = ii * ds;
272:       xi[dim*pcnt+1] = jj * ds;

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

277:       xi[dim*pcnt+0] += 0.35*ds2;
278:       xi[dim*pcnt+1] += 0.35*ds2;
279:       pcnt++;
280:     }
281:   }
282:   npoints_q = pcnt;

284:   npe = 3; /* nodes per element (triangle) */
285:   PetscMalloc1(npoints_q,&basis);
286:   for (q=0; q<npoints_q; q++) {
287:     PetscMalloc1(npe,&basis[q]);

289:     basis[q][0] = 1.0 - xi[dim*q+0] - xi[dim*q+1];
290:     basis[q][1] = xi[dim*q+0];
291:     basis[q][2] = xi[dim*q+1];
292:   }

294:   /* 0->cell, 1->edge, 2->vert */
295:   DMPlexGetHeightStratum(dmc,0,&ps,&pe);
296:   nel = pe - ps;

298:   DMSwarmSetLocalSizes(dm,npoints_q*nel,-1);
299:   DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
300:   DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);

302:   DMGetCoordinatesLocal(dmc,&coorlocal);
303:   DMGetCoordinateSection(dmc,&coordSection);

305:   pcnt = 0;
306:   for (e=0; e<nel; e++) {
307:     DMPlexVecGetClosure(dmc,coordSection,coorlocal,e,NULL,&elcoor);

309:     for (q=0; q<npoints_q; q++) {
310:       for (d=0; d<dim; d++) {
311:         swarm_coor[dim*pcnt+d] = 0.0;
312:         for (k=0; k<npe; k++) {
313:           swarm_coor[dim*pcnt+d] += basis[q][k] * PetscRealPart(elcoor[dim*k+d]);
314:         }
315:       }
316:       swarm_cellid[pcnt] = e;
317:       pcnt++;
318:     }
319:     DMPlexVecRestoreClosure(dmc,coordSection,coorlocal,e,NULL,&elcoor);
320:   }
321:   DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
322:   DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);

324:   PetscFree(xi);
325:   for (q=0; q<npoints_q; q++) {
326:     PetscFree(basis[q]);
327:   }
328:   PetscFree(basis);
329:   return 0;
330: }

332: PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX(DM dm,DM celldm,DMSwarmPICLayoutType layout,PetscInt layout_param)
333: {
334:   PetscInt       dim;

336:   DMGetDimension(celldm,&dim);
337:   switch (layout) {
338:     case DMSWARMPIC_LAYOUT_REGULAR:
340:       private_DMSwarmInsertPointsUsingCellDM_PLEX2D_Regular(dm,celldm,layout_param);
341:       break;
342:     case DMSWARMPIC_LAYOUT_GAUSS:
343:     {
344:       PetscInt npoints,npoints1,ps,pe,nfaces;
345:       const PetscReal *xi;
346:       PetscBool is_simplex;
347:       PetscQuadrature quadrature;

349:       is_simplex = PETSC_FALSE;
350:       DMPlexGetHeightStratum(celldm,0,&ps,&pe);
351:       DMPlexGetConeSize(celldm, ps, &nfaces);
352:       if (nfaces == (dim+1)) { is_simplex = PETSC_TRUE; }

354:       npoints1 = layout_param;
355:       if (is_simplex) {
356:         PetscDTStroudConicalQuadrature(dim,1,npoints1,-1.0,1.0,&quadrature);
357:       } else {
358:         PetscDTGaussTensorQuadrature(dim,1,npoints1,-1.0,1.0,&quadrature);
359:       }
360:       PetscQuadratureGetData(quadrature,NULL,NULL,&npoints,&xi,NULL);
361:       private_DMSwarmSetPointCoordinatesCellwise_PLEX(dm,celldm,npoints,(PetscReal*)xi);
362:       PetscQuadratureDestroy(&quadrature);
363:     }
364:       break;
365:     case DMSWARMPIC_LAYOUT_SUBDIVISION:
366:       private_DMSwarmInsertPointsUsingCellDM_PLEX_SubDivide(dm,celldm,layout_param);
367:       break;
368:   }
369:   return 0;
370: }

372: /*
373: typedef struct {
374:   PetscReal x,y;
375: } Point2d;

377: static PetscErrorCode signp2d(Point2d p1, Point2d p2, Point2d p3,PetscReal *s)
378: {
379:   *s = (p1.x - p3.x) * (p2.y - p3.y) - (p2.x - p3.x) * (p1.y - p3.y);
380:   return 0;
381: }
382: */
383: /*
384: static PetscErrorCode PointInTriangle(Point2d pt, Point2d v1, Point2d v2, Point2d v3,PetscBool *v)
385: {
386:   PetscReal s1,s2,s3;
387:   PetscBool b1, b2, b3;

389:   signp2d(pt, v1, v2,&s1); b1 = s1 < 0.0f;
390:   signp2d(pt, v2, v3,&s2); b2 = s2 < 0.0f;
391:   signp2d(pt, v3, v1,&s3); b3 = s3 < 0.0f;

393:   *v = ((b1 == b2) && (b2 == b3));
394:   return 0;
395: }
396: */
397: /*
398: static PetscErrorCode _ComputeLocalCoordinateAffine2d(PetscReal xp[],PetscReal coords[],PetscReal xip[],PetscReal *dJ)
399: {
400:   PetscReal x1,y1,x2,y2,x3,y3;
401:   PetscReal c,b[2],A[2][2],inv[2][2],detJ,od;

403:   x1 = coords[2*0+0];
404:   x2 = coords[2*1+0];
405:   x3 = coords[2*2+0];

407:   y1 = coords[2*0+1];
408:   y2 = coords[2*1+1];
409:   y3 = coords[2*2+1];

411:   c = x1 - 0.5*x1 - 0.5*x1 + 0.5*x2 + 0.5*x3;
412:   b[0] = xp[0] - c;
413:   c = y1 - 0.5*y1 - 0.5*y1 + 0.5*y2 + 0.5*y3;
414:   b[1] = xp[1] - c;

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

419:   detJ = A[0][0]*A[1][1] - A[0][1]*A[1][0];
420:   *dJ = PetscAbsReal(detJ);
421:   od = 1.0/detJ;

423:   inv[0][0] =  A[1][1] * od;
424:   inv[0][1] = -A[0][1] * od;
425:   inv[1][0] = -A[1][0] * od;
426:   inv[1][1] =  A[0][0] * od;

428:   xip[0] = inv[0][0]*b[0] + inv[0][1]*b[1];
429:   xip[1] = inv[1][0]*b[0] + inv[1][1]*b[1];
430:   return 0;
431: }
432: */

434: static PetscErrorCode ComputeLocalCoordinateAffine2d(PetscReal xp[],PetscScalar coords[],PetscReal xip[],PetscReal *dJ)
435: {
436:   PetscReal x1,y1,x2,y2,x3,y3;
437:   PetscReal b[2],A[2][2],inv[2][2],detJ,od;

439:   x1 = PetscRealPart(coords[2*0+0]);
440:   x2 = PetscRealPart(coords[2*1+0]);
441:   x3 = PetscRealPart(coords[2*2+0]);

443:   y1 = PetscRealPart(coords[2*0+1]);
444:   y2 = PetscRealPart(coords[2*1+1]);
445:   y3 = PetscRealPart(coords[2*2+1]);

447:   b[0] = xp[0] - x1;
448:   b[1] = xp[1] - y1;

450:   A[0][0] = x2-x1;   A[0][1] = x3-x1;
451:   A[1][0] = y2-y1;   A[1][1] = y3-y1;

453:   detJ = A[0][0]*A[1][1] - A[0][1]*A[1][0];
454:   *dJ = PetscAbsReal(detJ);
455:   od = 1.0/detJ;

457:   inv[0][0] =  A[1][1] * od;
458:   inv[0][1] = -A[0][1] * od;
459:   inv[1][0] = -A[1][0] * od;
460:   inv[1][1] =  A[0][0] * od;

462:   xip[0] = inv[0][0]*b[0] + inv[0][1]*b[1];
463:   xip[1] = inv[1][0]*b[0] + inv[1][1]*b[1];
464:   return 0;
465: }

467: PetscErrorCode DMSwarmProjectField_ApproxP1_PLEX_2D(DM swarm,PetscReal *swarm_field,DM dm,Vec v_field)
468: {
469:   const PetscReal PLEX_C_EPS = 1.0e-8;
470:   Vec             v_field_l,denom_l,coor_l,denom;
471:   PetscInt        k,p,e,npoints;
472:   PetscInt        *mpfield_cell;
473:   PetscReal       *mpfield_coor;
474:   PetscReal       xi_p[2];
475:   PetscScalar     Ni[3];
476:   PetscSection    coordSection;
477:   PetscScalar     *elcoor = NULL;

479:   VecZeroEntries(v_field);

481:   DMGetLocalVector(dm,&v_field_l);
482:   DMGetGlobalVector(dm,&denom);
483:   DMGetLocalVector(dm,&denom_l);
484:   VecZeroEntries(v_field_l);
485:   VecZeroEntries(denom);
486:   VecZeroEntries(denom_l);

488:   DMGetCoordinatesLocal(dm,&coor_l);
489:   DMGetCoordinateSection(dm,&coordSection);

491:   DMSwarmGetLocalSize(swarm,&npoints);
492:   DMSwarmGetField(swarm,DMSwarmPICField_coor,NULL,NULL,(void**)&mpfield_coor);
493:   DMSwarmGetField(swarm,DMSwarmPICField_cellid,NULL,NULL,(void**)&mpfield_cell);

495:   for (p=0; p<npoints; p++) {
496:     PetscReal   *coor_p,dJ;
497:     PetscScalar elfield[3];
498:     PetscBool   point_located;

500:     e       = mpfield_cell[p];
501:     coor_p  = &mpfield_coor[2*p];

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

505: /*
506:     while (!point_located && (failed_counter < 25)) {
507:       PointInTriangle(point, coords[0], coords[1], coords[2], &point_located);
508:       point.x = coor_p[0];
509:       point.y = coor_p[1];
510:       point.x += 1.0e-10 * (2.0 * rand()/((double)RAND_MAX)-1.0);
511:       point.y += 1.0e-10 * (2.0 * rand()/((double)RAND_MAX)-1.0);
512:       failed_counter++;
513:     }

515:     if (!point_located) {
516:         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);
517:     }

520:     else {
521:       _ComputeLocalCoordinateAffine2d(coor_p,elcoor,xi_p,&dJ);
522:       xi_p[0] = 0.5*(xi_p[0] + 1.0);
523:       xi_p[1] = 0.5*(xi_p[1] + 1.0);

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

527:     }
528: */

530:     ComputeLocalCoordinateAffine2d(coor_p,elcoor,xi_p,&dJ);
531:     /*
532:     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]);
533:     */
534:     /*
535:      point_located = PETSC_TRUE;
536:     if (xi_p[0] < 0.0) {
537:       if (xi_p[0] > -PLEX_C_EPS) {
538:         xi_p[0] = 0.0;
539:       } else {
540:         point_located = PETSC_FALSE;
541:       }
542:     }
543:     if (xi_p[1] < 0.0) {
544:       if (xi_p[1] > -PLEX_C_EPS) {
545:         xi_p[1] = 0.0;
546:       } else {
547:         point_located = PETSC_FALSE;
548:       }
549:     }
550:     if (xi_p[1] > (1.0-xi_p[0])) {
551:       if ((xi_p[1] - 1.0 + xi_p[0]) < PLEX_C_EPS) {
552:         xi_p[1] = 1.0 - xi_p[0];
553:       } else {
554:         point_located = PETSC_FALSE;
555:       }
556:     }
557:     if (!point_located) {
558:       PetscPrintf(PETSC_COMM_SELF,"[Error] xi,eta = %+1.8e, %+1.8e\n",xi_p[0],xi_p[1]);
559:       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]);
560:     }
562:     */

564:     Ni[0] = 1.0 - xi_p[0] - xi_p[1];
565:     Ni[1] = xi_p[0];
566:     Ni[2] = xi_p[1];

568:     point_located = PETSC_TRUE;
569:     for (k=0; k<3; k++) {
570:       if (PetscRealPart(Ni[k]) < -PLEX_C_EPS) point_located = PETSC_FALSE;
571:       if (PetscRealPart(Ni[k]) > (1.0+PLEX_C_EPS)) point_located = PETSC_FALSE;
572:     }
573:     if (!point_located) {
574:       PetscPrintf(PETSC_COMM_SELF,"[Error] xi,eta = %+1.8e, %+1.8e\n",(double)xi_p[0],(double)xi_p[1]);
575:       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]));
576:     }

579:     for (k=0; k<3; k++) {
580:       Ni[k] = Ni[k] * dJ;
581:       elfield[k] = Ni[k] * swarm_field[p];
582:     }
583:     DMPlexVecRestoreClosure(dm,coordSection,coor_l,e,NULL,&elcoor);

585:     DMPlexVecSetClosure(dm, NULL,v_field_l, e, elfield, ADD_VALUES);
586:     DMPlexVecSetClosure(dm, NULL,denom_l, e, Ni, ADD_VALUES);
587:   }

589:   DMSwarmRestoreField(swarm,DMSwarmPICField_cellid,NULL,NULL,(void**)&mpfield_cell);
590:   DMSwarmRestoreField(swarm,DMSwarmPICField_coor,NULL,NULL,(void**)&mpfield_coor);

592:   DMLocalToGlobalBegin(dm,v_field_l,ADD_VALUES,v_field);
593:   DMLocalToGlobalEnd(dm,v_field_l,ADD_VALUES,v_field);
594:   DMLocalToGlobalBegin(dm,denom_l,ADD_VALUES,denom);
595:   DMLocalToGlobalEnd(dm,denom_l,ADD_VALUES,denom);

597:   VecPointwiseDivide(v_field,v_field,denom);

599:   DMRestoreLocalVector(dm,&v_field_l);
600:   DMRestoreLocalVector(dm,&denom_l);
601:   DMRestoreGlobalVector(dm,&denom);
602:   return 0;
603: }

605: PetscErrorCode private_DMSwarmProjectFields_PLEX(DM swarm,DM celldm,PetscInt project_type,PetscInt nfields,DMSwarmDataField dfield[],Vec vecs[])
606: {
607:   PetscInt       f,dim;

609:   DMGetDimension(swarm,&dim);
610:   switch (dim) {
611:     case 2:
612:       for (f=0; f<nfields; f++) {
613:         PetscReal *swarm_field;

615:         DMSwarmDataFieldGetEntries(dfield[f],(void**)&swarm_field);
616:         DMSwarmProjectField_ApproxP1_PLEX_2D(swarm,swarm_field,celldm,vecs[f]);
617:       }
618:       break;
619:     case 3:
620:       SETERRQ(PetscObjectComm((PetscObject)swarm),PETSC_ERR_SUP,"No support for 3D");
621:     default:
622:       break;
623:   }
624:   return 0;
625: }

627: PetscErrorCode private_DMSwarmSetPointCoordinatesCellwise_PLEX(DM dm,DM dmc,PetscInt npoints,PetscReal xi[])
628: {
629:   PetscBool       is_simplex,is_tensorcell;
630:   PetscInt        dim,nfaces,ps,pe,p,d,nbasis,pcnt,e,k,nel;
631:   PetscFE         fe;
632:   PetscQuadrature quadrature;
633:   PetscTabulation T;
634:   PetscReal       *xiq;
635:   Vec             coorlocal;
636:   PetscSection    coordSection;
637:   PetscScalar     *elcoor = NULL;
638:   PetscReal       *swarm_coor;
639:   PetscInt        *swarm_cellid;

641:   DMGetDimension(dmc,&dim);

643:   is_simplex = PETSC_FALSE;
644:   is_tensorcell = PETSC_FALSE;
645:   DMPlexGetHeightStratum(dmc,0,&ps,&pe);
646:   DMPlexGetConeSize(dmc, ps, &nfaces);

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

650:   switch (dim) {
651:     case 2:
652:       if (nfaces == 4) { is_tensorcell = PETSC_TRUE; }
653:       break;
654:     case 3:
655:       if (nfaces == 6) { is_tensorcell = PETSC_TRUE; }
656:       break;
657:     default:
658:       SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only support for 2D, 3D");
659:   }

661:   /* check points provided fail inside the reference cell */
662:   if (is_simplex) {
663:     for (p=0; p<npoints; p++) {
664:       PetscReal sum;
665:       for (d=0; d<dim; d++) {
667:       }
668:       sum = 0.0;
669:       for (d=0; d<dim; d++) {
670:         sum += xi[dim*p+d];
671:       }
673:     }
674:   } else if (is_tensorcell) {
675:     for (p=0; p<npoints; p++) {
676:       for (d=0; d<dim; d++) {
678:       }
679:     }
680:   } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only support for d-simplex and d-tensorcell");

682:   PetscQuadratureCreate(PetscObjectComm((PetscObject)dm),&quadrature);
683:   PetscMalloc1(npoints*dim,&xiq);
684:   PetscArraycpy(xiq,xi,npoints*dim);
685:   PetscQuadratureSetData(quadrature,dim,1,npoints,(const PetscReal*)xiq,NULL);
686:   private_PetscFECreateDefault_scalar_pk1(dmc, dim, is_simplex, 0, &fe);
687:   PetscFESetQuadrature(fe,quadrature);
688:   PetscFEGetDimension(fe,&nbasis);
689:   PetscFEGetCellTabulation(fe, 1, &T);

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

696:   DMSwarmSetLocalSizes(dm,npoints*nel,-1);
697:   DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
698:   DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);

700:   DMGetCoordinatesLocal(dmc,&coorlocal);
701:   DMGetCoordinateSection(dmc,&coordSection);

703:   pcnt = 0;
704:   for (e=0; e<nel; e++) {
705:     DMPlexVecGetClosure(dmc,coordSection,coorlocal,ps+e,NULL,&elcoor);

707:     for (p=0; p<npoints; p++) {
708:       for (d=0; d<dim; d++) {
709:         swarm_coor[dim*pcnt+d] = 0.0;
710:         for (k=0; k<nbasis; k++) {
711:           swarm_coor[dim*pcnt+d] += T->T[0][p*nbasis + k] * PetscRealPart(elcoor[dim*k+d]);
712:         }
713:       }
714:       swarm_cellid[pcnt] = e;
715:       pcnt++;
716:     }
717:     DMPlexVecRestoreClosure(dmc,coordSection,coorlocal,ps+e,NULL,&elcoor);
718:   }
719:   DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
720:   DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);

722:   PetscQuadratureDestroy(&quadrature);
723:   PetscFEDestroy(&fe);
724:   return 0;
725: }