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

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

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

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

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

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

118:   PetscFunctionBegin;
119:   npoints_q = 1;
120:   for (d = 0; d < nsub; d++) npoints_q *= 4;
121:   PetscCall(PetscMalloc1(dim * npoints_q, &xi));

123:   v1[0] = 0.0;
124:   v1[1] = 0.0;
125:   v2[0] = 1.0;
126:   v2[1] = 0.0;
127:   v3[0] = 0.0;
128:   v3[1] = 1.0;
129:   depth = 0;
130:   pcnt  = 0;
131:   PetscCall(subdivide_triangle(v1, v2, v3, depth, nsub, xi, &pcnt));

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

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

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

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

151:   PetscCall(DMGetCoordinatesLocal(dmc, &coorlocal));
152:   PetscCall(DMGetCoordinateSection(dmc, &coordSection));

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

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

171:   PetscCall(PetscFree(xi));
172:   for (q = 0; q < npoints_q; q++) PetscCall(PetscFree(basis[q]));
173:   PetscCall(PetscFree(basis));
174:   PetscFunctionReturn(PETSC_SUCCESS);
175: }

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

192:   PetscFunctionBegin;
193:   PetscCall(DMGetDimension(dmc, &dim));
194:   is_simplex = PETSC_FALSE;
195:   PetscCall(DMPlexGetHeightStratum(dmc, 0, &ps, &pe));
196:   PetscCall(DMPlexGetConeSize(dmc, ps, &nfaces));
197:   if (nfaces == (dim + 1)) is_simplex = PETSC_TRUE;

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

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

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

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

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

220:   PetscCall(DMGetCoordinatesLocal(dmc, &coorlocal));
221:   PetscCall(DMGetCoordinateSection(dmc, &coordSection));

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

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

240:   PetscCall(PetscFEDestroy(&fe));
241:   PetscFunctionReturn(PETSC_SUCCESS);
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:   PetscFunctionBegin;
258:   PetscCall(DMGetDimension(dmc, &dim));
259:   PetscCheck(dim == 2, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only 2D is supported");
260:   is_simplex = PETSC_FALSE;
261:   PetscCall(DMPlexGetHeightStratum(dmc, 0, &ps, &pe));
262:   PetscCall(DMPlexGetConeSize(dmc, ps, &nfaces));
263:   if (nfaces == (dim + 1)) is_simplex = PETSC_TRUE;
264:   PetscCheck(is_simplex, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only the simplex is supported");

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

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

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

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

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

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

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

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

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

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

323:   PetscCall(PetscFree(xi));
324:   for (q = 0; q < npoints_q; q++) PetscCall(PetscFree(basis[q]));
325:   PetscCall(PetscFree(basis));
326:   PetscFunctionReturn(PETSC_SUCCESS);
327: }

329: PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX(DM dm, DM celldm, DMSwarmPICLayoutType layout, PetscInt layout_param)
330: {
331:   PetscInt dim;

333:   PetscFunctionBegin;
334:   PetscCall(DMGetDimension(celldm, &dim));
335:   switch (layout) {
336:   case DMSWARMPIC_LAYOUT_REGULAR:
337:     PetscCheck(dim != 3, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "No 3D support for REGULAR+PLEX");
338:     PetscCall(private_DMSwarmInsertPointsUsingCellDM_PLEX2D_Regular(dm, celldm, layout_param));
339:     break;
340:   case DMSWARMPIC_LAYOUT_GAUSS: {
341:     PetscInt         npoints, npoints1, ps, pe, nfaces;
342:     const PetscReal *xi;
343:     PetscBool        is_simplex;
344:     PetscQuadrature  quadrature;

346:     is_simplex = PETSC_FALSE;
347:     PetscCall(DMPlexGetHeightStratum(celldm, 0, &ps, &pe));
348:     PetscCall(DMPlexGetConeSize(celldm, ps, &nfaces));
349:     if (nfaces == (dim + 1)) is_simplex = PETSC_TRUE;

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

368: /*
369: typedef struct {
370:   PetscReal x,y;
371: } Point2d;

373: static PetscErrorCode signp2d(Point2d p1, Point2d p2, Point2d p3,PetscReal *s)
374: {
375:   PetscFunctionBegin;
376:   *s = (p1.x - p3.x) * (p2.y - p3.y) - (p2.x - p3.x) * (p1.y - p3.y);
377:   PetscFunctionReturn(PETSC_SUCCESS);
378: }
379: */
380: /*
381: static PetscErrorCode PointInTriangle(Point2d pt, Point2d v1, Point2d v2, Point2d v3,PetscBool *v)
382: {
383:   PetscReal s1,s2,s3;
384:   PetscBool b1, b2, b3;

386:   PetscFunctionBegin;
387:   signp2d(pt, v1, v2,&s1); b1 = s1 < 0.0f;
388:   signp2d(pt, v2, v3,&s2); b2 = s2 < 0.0f;
389:   signp2d(pt, v3, v1,&s3); b3 = s3 < 0.0f;

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

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

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

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

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

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

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

427:   xip[0] = inv[0][0]*b[0] + inv[0][1]*b[1];
428:   xip[1] = inv[1][0]*b[0] + inv[1][1]*b[1];
429:   PetscFunctionReturn(PETSC_SUCCESS);
430: }
431: */

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

438:   PetscFunctionBegin;
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;
451:   A[0][1] = x3 - x1;
452:   A[1][0] = y2 - y1;
453:   A[1][1] = y3 - y1;

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

459:   inv[0][0] = A[1][1] * od;
460:   inv[0][1] = -A[0][1] * od;
461:   inv[1][0] = -A[1][0] * od;
462:   inv[1][1] = A[0][0] * od;

464:   xip[0] = inv[0][0] * b[0] + inv[0][1] * b[1];
465:   xip[1] = inv[1][0] * b[0] + inv[1][1] * b[1];
466:   PetscFunctionReturn(PETSC_SUCCESS);
467: }

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

481:   PetscFunctionBegin;
482:   PetscCall(VecZeroEntries(v_field));

484:   PetscCall(DMGetLocalVector(dm, &v_field_l));
485:   PetscCall(DMGetGlobalVector(dm, &denom));
486:   PetscCall(DMGetLocalVector(dm, &denom_l));
487:   PetscCall(VecZeroEntries(v_field_l));
488:   PetscCall(VecZeroEntries(denom));
489:   PetscCall(VecZeroEntries(denom_l));

491:   PetscCall(DMGetCoordinatesLocal(dm, &coor_l));
492:   PetscCall(DMGetCoordinateSection(dm, &coordSection));

494:   PetscCall(DMSwarmGetLocalSize(swarm, &npoints));
495:   PetscCall(DMSwarmGetField(swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&mpfield_coor));
496:   PetscCall(DMSwarmGetField(swarm, DMSwarmPICField_cellid, NULL, NULL, (void **)&mpfield_cell));

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

503:     e      = mpfield_cell[p];
504:     coor_p = &mpfield_coor[2 * p];

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

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

518:     if (!point_located) {
519:         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);
520:     }

522:     PetscCheck(point_located,PETSC_COMM_SELF,PETSC_ERR_SUP,"Failed to locate point (%1.8e,%1.8e) in local mesh (cell %" PetscInt_FMT ")",point.x,point.y,e);
523:     else {
524:       PetscCall(_ComputeLocalCoordinateAffine2d(coor_p,elcoor,xi_p,&dJ));
525:       xi_p[0] = 0.5*(xi_p[0] + 1.0);
526:       xi_p[1] = 0.5*(xi_p[1] + 1.0);

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

530:     }
531: */

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

567:     Ni[0] = 1.0 - xi_p[0] - xi_p[1];
568:     Ni[1] = xi_p[0];
569:     Ni[2] = xi_p[1];

571:     point_located = PETSC_TRUE;
572:     for (k = 0; k < 3; k++) {
573:       if (PetscRealPart(Ni[k]) < -PLEX_C_EPS) point_located = PETSC_FALSE;
574:       if (PetscRealPart(Ni[k]) > (1.0 + PLEX_C_EPS)) point_located = PETSC_FALSE;
575:     }
576:     if (!point_located) {
577:       PetscCall(PetscPrintf(PETSC_COMM_SELF, "[Error] xi,eta = %+1.8e, %+1.8e\n", (double)xi_p[0], (double)xi_p[1]));
578:       PetscCall(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])));
579:     }
580:     PetscCheck(point_located, PETSC_COMM_SELF, PETSC_ERR_SUP, "Failed to locate point (%1.8e,%1.8e) in local mesh (cell %" PetscInt_FMT ")", (double)coor_p[0], (double)coor_p[1], e);

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

588:     PetscCall(DMPlexVecSetClosure(dm, NULL, v_field_l, e, elfield, ADD_VALUES));
589:     PetscCall(DMPlexVecSetClosure(dm, NULL, denom_l, e, Ni, ADD_VALUES));
590:   }

592:   PetscCall(DMSwarmRestoreField(swarm, DMSwarmPICField_cellid, NULL, NULL, (void **)&mpfield_cell));
593:   PetscCall(DMSwarmRestoreField(swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&mpfield_coor));

595:   PetscCall(DMLocalToGlobalBegin(dm, v_field_l, ADD_VALUES, v_field));
596:   PetscCall(DMLocalToGlobalEnd(dm, v_field_l, ADD_VALUES, v_field));
597:   PetscCall(DMLocalToGlobalBegin(dm, denom_l, ADD_VALUES, denom));
598:   PetscCall(DMLocalToGlobalEnd(dm, denom_l, ADD_VALUES, denom));

600:   PetscCall(VecPointwiseDivide(v_field, v_field, denom));

602:   PetscCall(DMRestoreLocalVector(dm, &v_field_l));
603:   PetscCall(DMRestoreLocalVector(dm, &denom_l));
604:   PetscCall(DMRestoreGlobalVector(dm, &denom));
605:   PetscFunctionReturn(PETSC_SUCCESS);
606: }

608: PetscErrorCode private_DMSwarmProjectFields_PLEX(DM swarm, DM celldm, PetscInt project_type, PetscInt nfields, DMSwarmDataField dfield[], Vec vecs[])
609: {
610:   PetscInt f, dim;

612:   PetscFunctionBegin;
613:   PetscCall(DMGetDimension(swarm, &dim));
614:   switch (dim) {
615:   case 2:
616:     for (f = 0; f < nfields; f++) {
617:       PetscReal *swarm_field;

619:       PetscCall(DMSwarmDataFieldGetEntries(dfield[f], (void **)&swarm_field));
620:       PetscCall(DMSwarmProjectField_ApproxP1_PLEX_2D(swarm, swarm_field, celldm, vecs[f]));
621:     }
622:     break;
623:   case 3:
624:     SETERRQ(PetscObjectComm((PetscObject)swarm), PETSC_ERR_SUP, "No support for 3D");
625:   default:
626:     break;
627:   }
628:   PetscFunctionReturn(PETSC_SUCCESS);
629: }

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

645:   PetscFunctionBegin;
646:   PetscCall(DMGetDimension(dmc, &dim));

648:   is_simplex    = PETSC_FALSE;
649:   is_tensorcell = PETSC_FALSE;
650:   PetscCall(DMPlexGetHeightStratum(dmc, 0, &ps, &pe));
651:   PetscCall(DMPlexGetConeSize(dmc, ps, &nfaces));

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

655:   switch (dim) {
656:   case 2:
657:     if (nfaces == 4) is_tensorcell = PETSC_TRUE;
658:     break;
659:   case 3:
660:     if (nfaces == 6) is_tensorcell = PETSC_TRUE;
661:     break;
662:   default:
663:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only support for 2D, 3D");
664:   }

666:   /* check points provided fail inside the reference cell */
667:   if (is_simplex) {
668:     for (p = 0; p < npoints; p++) {
669:       PetscReal sum;
670:       for (d = 0; d < dim; d++) PetscCheck(xi[dim * p + d] >= -1.0, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Points do not fail inside the simplex domain");
671:       sum = 0.0;
672:       for (d = 0; d < dim; d++) sum += xi[dim * p + d];
673:       PetscCheck(sum <= 0.0, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Points do not fail inside the simplex domain");
674:     }
675:   } else if (is_tensorcell) {
676:     for (p = 0; p < npoints; p++) {
677:       for (d = 0; d < dim; d++) PetscCheck(PetscAbsReal(xi[dim * p + d]) <= 1.0, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Points do not fail inside the tensor domain [-1,1]^d");
678:     }
679:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only support for d-simplex and d-tensorcell");

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

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

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

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

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

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

719:   PetscCall(PetscQuadratureDestroy(&quadrature));
720:   PetscCall(PetscFEDestroy(&fe));
721:   PetscFunctionReturn(PETSC_SUCCESS);
722: }