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:   } else {
 61:     PetscDTGaussTensorQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);
 62:     PetscDTGaussTensorQuadrature(dim - 1, 1, quadPointsPerEdge, -1.0, 1.0, &fq);
 63:   }
 64:   PetscFESetQuadrature(*fem, q);
 65:   PetscFESetFaceQuadrature(*fem, fq);
 66:   PetscQuadratureDestroy(&q);
 67:   PetscQuadratureDestroy(&fq);
 68:   return 0;
 69: }

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

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

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

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

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

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

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

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

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

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

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

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

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

148:   DMGetCoordinatesLocal(dmc, &coorlocal);
149:   DMGetCoordinateSection(dmc, &coordSection);

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

155:     for (q = 0; q < npoints_q; q++) {
156:       for (d = 0; d < dim; d++) {
157:         swarm_coor[dim * pcnt + d] = 0.0;
158:         for (k = 0; k < npe; k++) swarm_coor[dim * pcnt + d] += basis[q][k] * PetscRealPart(elcoor[dim * k + d]);
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++) PetscFree(basis[q]);
170:   PetscFree(basis);
171:   return 0;
172: }

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

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

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

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

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

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

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

216:   DMGetCoordinatesLocal(dmc, &coorlocal);
217:   DMGetCoordinateSection(dmc, &coordSection);

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

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

236:   PetscFEDestroy(&fe);
237:   return 0;
238: }

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

253:   DMGetDimension(dmc, &dim);
255:   is_simplex = PETSC_FALSE;
256:   DMPlexGetHeightStratum(dmc, 0, &ps, &pe);
257:   DMPlexGetConeSize(dmc, ps, &nfaces);
258:   if (nfaces == (dim + 1)) is_simplex = PETSC_TRUE;

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

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

273:       xi[dim * pcnt + 0] += 0.35 * ds2;
274:       xi[dim * pcnt + 1] += 0.35 * ds2;
275:       pcnt++;
276:     }
277:   }
278:   npoints_q = pcnt;

280:   npe = 3; /* nodes per element (triangle) */
281:   PetscMalloc1(npoints_q, &basis);
282:   for (q = 0; q < npoints_q; q++) {
283:     PetscMalloc1(npe, &basis[q]);

285:     basis[q][0] = 1.0 - xi[dim * q + 0] - xi[dim * q + 1];
286:     basis[q][1] = xi[dim * q + 0];
287:     basis[q][2] = xi[dim * q + 1];
288:   }

290:   /* 0->cell, 1->edge, 2->vert */
291:   DMPlexGetHeightStratum(dmc, 0, &ps, &pe);
292:   nel = pe - ps;

294:   DMSwarmSetLocalSizes(dm, npoints_q * nel, -1);
295:   DMSwarmGetField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor);
296:   DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid);

298:   DMGetCoordinatesLocal(dmc, &coorlocal);
299:   DMGetCoordinateSection(dmc, &coordSection);

301:   pcnt = 0;
302:   for (e = 0; e < nel; e++) {
303:     DMPlexVecGetClosure(dmc, coordSection, coorlocal, e, NULL, &elcoor);

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

318:   PetscFree(xi);
319:   for (q = 0; q < npoints_q; q++) PetscFree(basis[q]);
320:   PetscFree(basis);
321:   return 0;
322: }

324: PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX(DM dm, DM celldm, DMSwarmPICLayoutType layout, PetscInt layout_param)
325: {
326:   PetscInt dim;

328:   DMGetDimension(celldm, &dim);
329:   switch (layout) {
330:   case DMSWARMPIC_LAYOUT_REGULAR:
332:     private_DMSwarmInsertPointsUsingCellDM_PLEX2D_Regular(dm, celldm, layout_param);
333:     break;
334:   case DMSWARMPIC_LAYOUT_GAUSS: {
335:     PetscInt         npoints, npoints1, ps, pe, nfaces;
336:     const PetscReal *xi;
337:     PetscBool        is_simplex;
338:     PetscQuadrature  quadrature;

340:     is_simplex = PETSC_FALSE;
341:     DMPlexGetHeightStratum(celldm, 0, &ps, &pe);
342:     DMPlexGetConeSize(celldm, ps, &nfaces);
343:     if (nfaces == (dim + 1)) is_simplex = PETSC_TRUE;

345:     npoints1 = layout_param;
346:     if (is_simplex) {
347:       PetscDTStroudConicalQuadrature(dim, 1, npoints1, -1.0, 1.0, &quadrature);
348:     } else {
349:       PetscDTGaussTensorQuadrature(dim, 1, npoints1, -1.0, 1.0, &quadrature);
350:     }
351:     PetscQuadratureGetData(quadrature, NULL, NULL, &npoints, &xi, NULL);
352:     private_DMSwarmSetPointCoordinatesCellwise_PLEX(dm, celldm, npoints, (PetscReal *)xi);
353:     PetscQuadratureDestroy(&quadrature);
354:   } break;
355:   case DMSWARMPIC_LAYOUT_SUBDIVISION:
356:     private_DMSwarmInsertPointsUsingCellDM_PLEX_SubDivide(dm, celldm, layout_param);
357:     break;
358:   }
359:   return 0;
360: }

362: /*
363: typedef struct {
364:   PetscReal x,y;
365: } Point2d;

367: static PetscErrorCode signp2d(Point2d p1, Point2d p2, Point2d p3,PetscReal *s)
368: {
369:   *s = (p1.x - p3.x) * (p2.y - p3.y) - (p2.x - p3.x) * (p1.y - p3.y);
370:   return 0;
371: }
372: */
373: /*
374: static PetscErrorCode PointInTriangle(Point2d pt, Point2d v1, Point2d v2, Point2d v3,PetscBool *v)
375: {
376:   PetscReal s1,s2,s3;
377:   PetscBool b1, b2, b3;

379:   signp2d(pt, v1, v2,&s1); b1 = s1 < 0.0f;
380:   signp2d(pt, v2, v3,&s2); b2 = s2 < 0.0f;
381:   signp2d(pt, v3, v1,&s3); b3 = s3 < 0.0f;

383:   *v = ((b1 == b2) && (b2 == b3));
384:   return 0;
385: }
386: */
387: /*
388: static PetscErrorCode _ComputeLocalCoordinateAffine2d(PetscReal xp[],PetscReal coords[],PetscReal xip[],PetscReal *dJ)
389: {
390:   PetscReal x1,y1,x2,y2,x3,y3;
391:   PetscReal c,b[2],A[2][2],inv[2][2],detJ,od;

393:   x1 = coords[2*0+0];
394:   x2 = coords[2*1+0];
395:   x3 = coords[2*2+0];

397:   y1 = coords[2*0+1];
398:   y2 = coords[2*1+1];
399:   y3 = coords[2*2+1];

401:   c = x1 - 0.5*x1 - 0.5*x1 + 0.5*x2 + 0.5*x3;
402:   b[0] = xp[0] - c;
403:   c = y1 - 0.5*y1 - 0.5*y1 + 0.5*y2 + 0.5*y3;
404:   b[1] = xp[1] - c;

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

409:   detJ = A[0][0]*A[1][1] - A[0][1]*A[1][0];
410:   *dJ = PetscAbsReal(detJ);
411:   od = 1.0/detJ;

413:   inv[0][0] =  A[1][1] * od;
414:   inv[0][1] = -A[0][1] * od;
415:   inv[1][0] = -A[1][0] * od;
416:   inv[1][1] =  A[0][0] * od;

418:   xip[0] = inv[0][0]*b[0] + inv[0][1]*b[1];
419:   xip[1] = inv[1][0]*b[0] + inv[1][1]*b[1];
420:   return 0;
421: }
422: */

424: static PetscErrorCode ComputeLocalCoordinateAffine2d(PetscReal xp[], PetscScalar coords[], PetscReal xip[], PetscReal *dJ)
425: {
426:   PetscReal x1, y1, x2, y2, x3, y3;
427:   PetscReal b[2], A[2][2], inv[2][2], detJ, od;

429:   x1 = PetscRealPart(coords[2 * 0 + 0]);
430:   x2 = PetscRealPart(coords[2 * 1 + 0]);
431:   x3 = PetscRealPart(coords[2 * 2 + 0]);

433:   y1 = PetscRealPart(coords[2 * 0 + 1]);
434:   y2 = PetscRealPart(coords[2 * 1 + 1]);
435:   y3 = PetscRealPart(coords[2 * 2 + 1]);

437:   b[0] = xp[0] - x1;
438:   b[1] = xp[1] - y1;

440:   A[0][0] = x2 - x1;
441:   A[0][1] = x3 - x1;
442:   A[1][0] = y2 - y1;
443:   A[1][1] = y3 - y1;

445:   detJ = A[0][0] * A[1][1] - A[0][1] * A[1][0];
446:   *dJ  = PetscAbsReal(detJ);
447:   od   = 1.0 / detJ;

449:   inv[0][0] = A[1][1] * od;
450:   inv[0][1] = -A[0][1] * od;
451:   inv[1][0] = -A[1][0] * od;
452:   inv[1][1] = A[0][0] * od;

454:   xip[0] = inv[0][0] * b[0] + inv[0][1] * b[1];
455:   xip[1] = inv[1][0] * b[0] + inv[1][1] * b[1];
456:   return 0;
457: }

459: PetscErrorCode DMSwarmProjectField_ApproxP1_PLEX_2D(DM swarm, PetscReal *swarm_field, DM dm, Vec v_field)
460: {
461:   const PetscReal PLEX_C_EPS = 1.0e-8;
462:   Vec             v_field_l, denom_l, coor_l, denom;
463:   PetscInt        k, p, e, npoints;
464:   PetscInt       *mpfield_cell;
465:   PetscReal      *mpfield_coor;
466:   PetscReal       xi_p[2];
467:   PetscScalar     Ni[3];
468:   PetscSection    coordSection;
469:   PetscScalar    *elcoor = NULL;

471:   VecZeroEntries(v_field);

473:   DMGetLocalVector(dm, &v_field_l);
474:   DMGetGlobalVector(dm, &denom);
475:   DMGetLocalVector(dm, &denom_l);
476:   VecZeroEntries(v_field_l);
477:   VecZeroEntries(denom);
478:   VecZeroEntries(denom_l);

480:   DMGetCoordinatesLocal(dm, &coor_l);
481:   DMGetCoordinateSection(dm, &coordSection);

483:   DMSwarmGetLocalSize(swarm, &npoints);
484:   DMSwarmGetField(swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&mpfield_coor);
485:   DMSwarmGetField(swarm, DMSwarmPICField_cellid, NULL, NULL, (void **)&mpfield_cell);

487:   for (p = 0; p < npoints; p++) {
488:     PetscReal  *coor_p, dJ;
489:     PetscScalar elfield[3];
490:     PetscBool   point_located;

492:     e      = mpfield_cell[p];
493:     coor_p = &mpfield_coor[2 * p];

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

497:     /*
498:     while (!point_located && (failed_counter < 25)) {
499:       PointInTriangle(point, coords[0], coords[1], coords[2], &point_located);
500:       point.x = coor_p[0];
501:       point.y = coor_p[1];
502:       point.x += 1.0e-10 * (2.0 * rand()/((double)RAND_MAX)-1.0);
503:       point.y += 1.0e-10 * (2.0 * rand()/((double)RAND_MAX)-1.0);
504:       failed_counter++;
505:     }

507:     if (!point_located) {
508:         PetscPrintf(PETSC_COMM_SELF,"Failed to locate point (%1.8e,%1.8e) in local mesh (cell %" PetscInt_FMT ") with triangle coords (%1.8e,%1.8e) : (%1.8e,%1.8e) : (%1.8e,%1.8e) in %" PetscInt_FMT " 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);
509:     }

512:     else {
513:       _ComputeLocalCoordinateAffine2d(coor_p,elcoor,xi_p,&dJ);
514:       xi_p[0] = 0.5*(xi_p[0] + 1.0);
515:       xi_p[1] = 0.5*(xi_p[1] + 1.0);

517:       PetscPrintf(PETSC_COMM_SELF,"[p=%" PetscInt_FMT "] x(%+1.4e,%+1.4e) -> mapped to element %" PetscInt_FMT " xi(%+1.4e,%+1.4e)\n",p,point.x,point.y,e,xi_p[0],xi_p[1]);

519:     }
520: */

522:     ComputeLocalCoordinateAffine2d(coor_p, elcoor, xi_p, &dJ);
523:     /*
524:     PetscPrintf(PETSC_COMM_SELF,"[p=%" PetscInt_FMT "] x(%+1.4e,%+1.4e) -> mapped to element %" PetscInt_FMT " xi(%+1.4e,%+1.4e)\n",p,point.x,point.y,e,xi_p[0],xi_p[1]);
525:     */
526:     /*
527:      point_located = PETSC_TRUE;
528:     if (xi_p[0] < 0.0) {
529:       if (xi_p[0] > -PLEX_C_EPS) {
530:         xi_p[0] = 0.0;
531:       } else {
532:         point_located = PETSC_FALSE;
533:       }
534:     }
535:     if (xi_p[1] < 0.0) {
536:       if (xi_p[1] > -PLEX_C_EPS) {
537:         xi_p[1] = 0.0;
538:       } else {
539:         point_located = PETSC_FALSE;
540:       }
541:     }
542:     if (xi_p[1] > (1.0-xi_p[0])) {
543:       if ((xi_p[1] - 1.0 + xi_p[0]) < PLEX_C_EPS) {
544:         xi_p[1] = 1.0 - xi_p[0];
545:       } else {
546:         point_located = PETSC_FALSE;
547:       }
548:     }
549:     if (!point_located) {
550:       PetscPrintf(PETSC_COMM_SELF,"[Error] xi,eta = %+1.8e, %+1.8e\n",xi_p[0],xi_p[1]);
551:       PetscPrintf(PETSC_COMM_SELF,"[Error] Failed to locate point (%1.8e,%1.8e) in local mesh (cell %" PetscInt_FMT ") 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]);
552:     }
554:     */

556:     Ni[0] = 1.0 - xi_p[0] - xi_p[1];
557:     Ni[1] = xi_p[0];
558:     Ni[2] = xi_p[1];

560:     point_located = PETSC_TRUE;
561:     for (k = 0; k < 3; k++) {
562:       if (PetscRealPart(Ni[k]) < -PLEX_C_EPS) point_located = PETSC_FALSE;
563:       if (PetscRealPart(Ni[k]) > (1.0 + PLEX_C_EPS)) point_located = PETSC_FALSE;
564:     }
565:     if (!point_located) {
566:       PetscPrintf(PETSC_COMM_SELF, "[Error] xi,eta = %+1.8e, %+1.8e\n", (double)xi_p[0], (double)xi_p[1]);
567:       PetscPrintf(PETSC_COMM_SELF, "[Error] Failed to locate point (%1.8e,%1.8e) in local mesh (cell %" PetscInt_FMT ") 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]));
568:     }

571:     for (k = 0; k < 3; k++) {
572:       Ni[k]      = Ni[k] * dJ;
573:       elfield[k] = Ni[k] * swarm_field[p];
574:     }
575:     DMPlexVecRestoreClosure(dm, coordSection, coor_l, e, NULL, &elcoor);

577:     DMPlexVecSetClosure(dm, NULL, v_field_l, e, elfield, ADD_VALUES);
578:     DMPlexVecSetClosure(dm, NULL, denom_l, e, Ni, ADD_VALUES);
579:   }

581:   DMSwarmRestoreField(swarm, DMSwarmPICField_cellid, NULL, NULL, (void **)&mpfield_cell);
582:   DMSwarmRestoreField(swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&mpfield_coor);

584:   DMLocalToGlobalBegin(dm, v_field_l, ADD_VALUES, v_field);
585:   DMLocalToGlobalEnd(dm, v_field_l, ADD_VALUES, v_field);
586:   DMLocalToGlobalBegin(dm, denom_l, ADD_VALUES, denom);
587:   DMLocalToGlobalEnd(dm, denom_l, ADD_VALUES, denom);

589:   VecPointwiseDivide(v_field, v_field, denom);

591:   DMRestoreLocalVector(dm, &v_field_l);
592:   DMRestoreLocalVector(dm, &denom_l);
593:   DMRestoreGlobalVector(dm, &denom);
594:   return 0;
595: }

597: PetscErrorCode private_DMSwarmProjectFields_PLEX(DM swarm, DM celldm, PetscInt project_type, PetscInt nfields, DMSwarmDataField dfield[], Vec vecs[])
598: {
599:   PetscInt f, dim;

601:   DMGetDimension(swarm, &dim);
602:   switch (dim) {
603:   case 2:
604:     for (f = 0; f < nfields; f++) {
605:       PetscReal *swarm_field;

607:       DMSwarmDataFieldGetEntries(dfield[f], (void **)&swarm_field);
608:       DMSwarmProjectField_ApproxP1_PLEX_2D(swarm, swarm_field, celldm, vecs[f]);
609:     }
610:     break;
611:   case 3:
612:     SETERRQ(PetscObjectComm((PetscObject)swarm), PETSC_ERR_SUP, "No support for 3D");
613:   default:
614:     break;
615:   }
616:   return 0;
617: }

619: PetscErrorCode private_DMSwarmSetPointCoordinatesCellwise_PLEX(DM dm, DM dmc, PetscInt npoints, PetscReal xi[])
620: {
621:   PetscBool       is_simplex, is_tensorcell;
622:   PetscInt        dim, nfaces, ps, pe, p, d, nbasis, pcnt, e, k, nel;
623:   PetscFE         fe;
624:   PetscQuadrature quadrature;
625:   PetscTabulation T;
626:   PetscReal      *xiq;
627:   Vec             coorlocal;
628:   PetscSection    coordSection;
629:   PetscScalar    *elcoor = NULL;
630:   PetscReal      *swarm_coor;
631:   PetscInt       *swarm_cellid;

633:   DMGetDimension(dmc, &dim);

635:   is_simplex    = PETSC_FALSE;
636:   is_tensorcell = PETSC_FALSE;
637:   DMPlexGetHeightStratum(dmc, 0, &ps, &pe);
638:   DMPlexGetConeSize(dmc, ps, &nfaces);

640:   if (nfaces == (dim + 1)) is_simplex = PETSC_TRUE;

642:   switch (dim) {
643:   case 2:
644:     if (nfaces == 4) is_tensorcell = PETSC_TRUE;
645:     break;
646:   case 3:
647:     if (nfaces == 6) is_tensorcell = PETSC_TRUE;
648:     break;
649:   default:
650:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only support for 2D, 3D");
651:   }

653:   /* check points provided fail inside the reference cell */
654:   if (is_simplex) {
655:     for (p = 0; p < npoints; p++) {
656:       PetscReal sum;
658:       sum = 0.0;
659:       for (d = 0; d < dim; d++) sum += xi[dim * p + d];
661:     }
662:   } else if (is_tensorcell) {
663:     for (p = 0; p < npoints; p++) {
665:     }
666:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only support for d-simplex and d-tensorcell");

668:   PetscQuadratureCreate(PetscObjectComm((PetscObject)dm), &quadrature);
669:   PetscMalloc1(npoints * dim, &xiq);
670:   PetscArraycpy(xiq, xi, npoints * dim);
671:   PetscQuadratureSetData(quadrature, dim, 1, npoints, (const PetscReal *)xiq, NULL);
672:   private_PetscFECreateDefault_scalar_pk1(dmc, dim, is_simplex, 0, &fe);
673:   PetscFESetQuadrature(fe, quadrature);
674:   PetscFEGetDimension(fe, &nbasis);
675:   PetscFEGetCellTabulation(fe, 1, &T);

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

682:   DMSwarmSetLocalSizes(dm, npoints * nel, -1);
683:   DMSwarmGetField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor);
684:   DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid);

686:   DMGetCoordinatesLocal(dmc, &coorlocal);
687:   DMGetCoordinateSection(dmc, &coordSection);

689:   pcnt = 0;
690:   for (e = 0; e < nel; e++) {
691:     DMPlexVecGetClosure(dmc, coordSection, coorlocal, ps + e, NULL, &elcoor);

693:     for (p = 0; p < npoints; p++) {
694:       for (d = 0; d < dim; d++) {
695:         swarm_coor[dim * pcnt + d] = 0.0;
696:         for (k = 0; k < nbasis; k++) swarm_coor[dim * pcnt + d] += T->T[0][p * nbasis + k] * PetscRealPart(elcoor[dim * k + d]);
697:       }
698:       swarm_cellid[pcnt] = e;
699:       pcnt++;
700:     }
701:     DMPlexVecRestoreClosure(dmc, coordSection, coorlocal, ps + e, NULL, &elcoor);
702:   }
703:   DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid);
704:   DMSwarmRestoreField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor);

706:   PetscQuadratureDestroy(&quadrature);
707:   PetscFEDestroy(&fe);
708:   return 0;
709: }