Actual source code: dminterpolatesnes.c

  1: #include <petsc/private/dmpleximpl.h>
  2: #include <petsc/private/snesimpl.h>
  3: #include <petscds.h>
  4: #include <petsc/private/petscimpl.h>
  5: #include <petsc/private/petscfeimpl.h>

  7: /*@C
  8:   DMInterpolationCreate - Creates a `DMInterpolationInfo` context

 10:   Collective

 12:   Input Parameter:
 13: . comm - the communicator

 15:   Output Parameter:
 16: . ctx - the context

 18:   Level: beginner

 20:   Developer Note:
 21:   The naming is incorrect, either the object should be named `DMInterpolation` or all the routines should begin with `DMInterpolationInfo`

 23: .seealso: [](ch_dmbase), `DM`, `DMInterpolationInfo`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationDestroy()`
 24: @*/
 25: PetscErrorCode DMInterpolationCreate(MPI_Comm comm, DMInterpolationInfo *ctx)
 26: {
 27:   PetscFunctionBegin;
 28:   PetscAssertPointer(ctx, 2);
 29:   PetscCall(PetscNew(ctx));

 31:   (*ctx)->comm   = comm;
 32:   (*ctx)->dim    = -1;
 33:   (*ctx)->nInput = 0;
 34:   (*ctx)->points = NULL;
 35:   (*ctx)->cells  = NULL;
 36:   (*ctx)->n      = -1;
 37:   (*ctx)->coords = NULL;
 38:   PetscFunctionReturn(PETSC_SUCCESS);
 39: }

 41: /*@C
 42:   DMInterpolationSetDim - Sets the spatial dimension for the interpolation context

 44:   Not Collective

 46:   Input Parameters:
 47: + ctx - the context
 48: - dim - the spatial dimension

 50:   Level: intermediate

 52: .seealso: [](ch_dmbase), `DM`, `DMInterpolationInfo`, `DMInterpolationGetDim()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`
 53: @*/
 54: PetscErrorCode DMInterpolationSetDim(DMInterpolationInfo ctx, PetscInt dim)
 55: {
 56:   PetscFunctionBegin;
 57:   PetscCheck(!(dim < 1) && !(dim > 3), ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension for points: %" PetscInt_FMT, dim);
 58:   ctx->dim = dim;
 59:   PetscFunctionReturn(PETSC_SUCCESS);
 60: }

 62: /*@C
 63:   DMInterpolationGetDim - Gets the spatial dimension for the interpolation context

 65:   Not Collective

 67:   Input Parameter:
 68: . ctx - the context

 70:   Output Parameter:
 71: . dim - the spatial dimension

 73:   Level: intermediate

 75: .seealso: [](ch_dmbase), `DM`, `DMInterpolationInfo`, `DMInterpolationSetDim()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`
 76: @*/
 77: PetscErrorCode DMInterpolationGetDim(DMInterpolationInfo ctx, PetscInt *dim)
 78: {
 79:   PetscFunctionBegin;
 80:   PetscAssertPointer(dim, 2);
 81:   *dim = ctx->dim;
 82:   PetscFunctionReturn(PETSC_SUCCESS);
 83: }

 85: /*@C
 86:   DMInterpolationSetDof - Sets the number of fields interpolated at a point for the interpolation context

 88:   Not Collective

 90:   Input Parameters:
 91: + ctx - the context
 92: - dof - the number of fields

 94:   Level: intermediate

 96: .seealso: [](ch_dmbase), `DM`, `DMInterpolationInfo`, `DMInterpolationGetDof()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`
 97: @*/
 98: PetscErrorCode DMInterpolationSetDof(DMInterpolationInfo ctx, PetscInt dof)
 99: {
100:   PetscFunctionBegin;
101:   PetscCheck(dof >= 1, ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid number of components: %" PetscInt_FMT, dof);
102:   ctx->dof = dof;
103:   PetscFunctionReturn(PETSC_SUCCESS);
104: }

106: /*@C
107:   DMInterpolationGetDof - Gets the number of fields interpolated at a point for the interpolation context

109:   Not Collective

111:   Input Parameter:
112: . ctx - the context

114:   Output Parameter:
115: . dof - the number of fields

117:   Level: intermediate

119: .seealso: [](ch_dmbase), `DM`, `DMInterpolationInfo`, `DMInterpolationSetDof()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`
120: @*/
121: PetscErrorCode DMInterpolationGetDof(DMInterpolationInfo ctx, PetscInt *dof)
122: {
123:   PetscFunctionBegin;
124:   PetscAssertPointer(dof, 2);
125:   *dof = ctx->dof;
126:   PetscFunctionReturn(PETSC_SUCCESS);
127: }

129: /*@C
130:   DMInterpolationAddPoints - Add points at which we will interpolate the fields

132:   Not Collective

134:   Input Parameters:
135: + ctx    - the context
136: . n      - the number of points
137: - points - the coordinates for each point, an array of size `n` * dim

139:   Level: intermediate

141:   Note:
142:   The input coordinate information is copied into the object.

144: .seealso: [](ch_dmbase), `DM`, `DMInterpolationInfo`, `DMInterpolationSetDim()`, `DMInterpolationEvaluate()`, `DMInterpolationCreate()`
145: @*/
146: PetscErrorCode DMInterpolationAddPoints(DMInterpolationInfo ctx, PetscInt n, PetscReal points[])
147: {
148:   PetscFunctionBegin;
149:   PetscCheck(ctx->dim >= 0, ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The spatial dimension has not been set");
150:   PetscCheck(!ctx->points, ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "Cannot add points multiple times yet");
151:   ctx->nInput = n;

153:   PetscCall(PetscMalloc1(n * ctx->dim, &ctx->points));
154:   PetscCall(PetscArraycpy(ctx->points, points, n * ctx->dim));
155:   PetscFunctionReturn(PETSC_SUCCESS);
156: }

158: /*@C
159:   DMInterpolationSetUp - Compute spatial indices for point location during interpolation

161:   Collective

163:   Input Parameters:
164: + ctx                 - the context
165: . dm                  - the `DM` for the function space used for interpolation
166: . redundantPoints     - If `PETSC_TRUE`, all processes are passing in the same array of points. Otherwise, points need to be communicated among processes.
167: - ignoreOutsideDomain - If `PETSC_TRUE`, ignore points outside the domain, otherwise return an error

169:   Level: intermediate

171: .seealso: [](ch_dmbase), `DM`, `DMInterpolationInfo`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`
172: @*/
173: PetscErrorCode DMInterpolationSetUp(DMInterpolationInfo ctx, DM dm, PetscBool redundantPoints, PetscBool ignoreOutsideDomain)
174: {
175:   MPI_Comm           comm = ctx->comm;
176:   PetscScalar       *a;
177:   PetscInt           p, q, i;
178:   PetscMPIInt        rank, size;
179:   Vec                pointVec;
180:   PetscSF            cellSF;
181:   PetscLayout        layout;
182:   PetscReal         *globalPoints;
183:   PetscScalar       *globalPointsScalar;
184:   const PetscInt    *ranges;
185:   PetscMPIInt       *counts, *displs;
186:   const PetscSFNode *foundCells;
187:   const PetscInt    *foundPoints;
188:   PetscMPIInt       *foundProcs, *globalProcs;
189:   PetscInt           n, N, numFound;

191:   PetscFunctionBegin;
193:   PetscCallMPI(MPI_Comm_size(comm, &size));
194:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
195:   PetscCheck(ctx->dim >= 0, comm, PETSC_ERR_ARG_WRONGSTATE, "The spatial dimension has not been set");
196:   /* Locate points */
197:   n = ctx->nInput;
198:   if (!redundantPoints) {
199:     PetscCall(PetscLayoutCreate(comm, &layout));
200:     PetscCall(PetscLayoutSetBlockSize(layout, 1));
201:     PetscCall(PetscLayoutSetLocalSize(layout, n));
202:     PetscCall(PetscLayoutSetUp(layout));
203:     PetscCall(PetscLayoutGetSize(layout, &N));
204:     /* Communicate all points to all processes */
205:     PetscCall(PetscMalloc3(N * ctx->dim, &globalPoints, size, &counts, size, &displs));
206:     PetscCall(PetscLayoutGetRanges(layout, &ranges));
207:     for (p = 0; p < size; ++p) {
208:       counts[p] = (ranges[p + 1] - ranges[p]) * ctx->dim;
209:       displs[p] = ranges[p] * ctx->dim;
210:     }
211:     PetscCallMPI(MPI_Allgatherv(ctx->points, n * ctx->dim, MPIU_REAL, globalPoints, counts, displs, MPIU_REAL, comm));
212:   } else {
213:     N            = n;
214:     globalPoints = ctx->points;
215:     counts = displs = NULL;
216:     layout          = NULL;
217:   }
218: #if 0
219:   PetscCall(PetscMalloc3(N,&foundCells,N,&foundProcs,N,&globalProcs));
220:   /* foundCells[p] = m->locatePoint(&globalPoints[p*ctx->dim]); */
221: #else
222:   #if defined(PETSC_USE_COMPLEX)
223:   PetscCall(PetscMalloc1(N * ctx->dim, &globalPointsScalar));
224:   for (i = 0; i < N * ctx->dim; i++) globalPointsScalar[i] = globalPoints[i];
225:   #else
226:   globalPointsScalar = globalPoints;
227:   #endif
228:   PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, ctx->dim, N * ctx->dim, globalPointsScalar, &pointVec));
229:   PetscCall(PetscMalloc2(N, &foundProcs, N, &globalProcs));
230:   for (p = 0; p < N; ++p) foundProcs[p] = size;
231:   cellSF = NULL;
232:   PetscCall(DMLocatePoints(dm, pointVec, DM_POINTLOCATION_REMOVE, &cellSF));
233:   PetscCall(PetscSFGetGraph(cellSF, NULL, &numFound, &foundPoints, &foundCells));
234: #endif
235:   for (p = 0; p < numFound; ++p) {
236:     if (foundCells[p].index >= 0) foundProcs[foundPoints ? foundPoints[p] : p] = rank;
237:   }
238:   /* Let the lowest rank process own each point */
239:   PetscCall(MPIU_Allreduce(foundProcs, globalProcs, N, MPI_INT, MPI_MIN, comm));
240:   ctx->n = 0;
241:   for (p = 0; p < N; ++p) {
242:     if (globalProcs[p] == size) {
243:       PetscCheck(ignoreOutsideDomain, comm, PETSC_ERR_PLIB, "Point %" PetscInt_FMT ": %g %g %g not located in mesh", p, (double)globalPoints[p * ctx->dim + 0], (double)(ctx->dim > 1 ? globalPoints[p * ctx->dim + 1] : 0.0),
244:                  (double)(ctx->dim > 2 ? globalPoints[p * ctx->dim + 2] : 0.0));
245:       if (rank == 0) ++ctx->n;
246:     } else if (globalProcs[p] == rank) ++ctx->n;
247:   }
248:   /* Create coordinates vector and array of owned cells */
249:   PetscCall(PetscMalloc1(ctx->n, &ctx->cells));
250:   PetscCall(VecCreate(comm, &ctx->coords));
251:   PetscCall(VecSetSizes(ctx->coords, ctx->n * ctx->dim, PETSC_DECIDE));
252:   PetscCall(VecSetBlockSize(ctx->coords, ctx->dim));
253:   PetscCall(VecSetType(ctx->coords, VECSTANDARD));
254:   PetscCall(VecGetArray(ctx->coords, &a));
255:   for (p = 0, q = 0, i = 0; p < N; ++p) {
256:     if (globalProcs[p] == rank) {
257:       PetscInt d;

259:       for (d = 0; d < ctx->dim; ++d, ++i) a[i] = globalPoints[p * ctx->dim + d];
260:       ctx->cells[q] = foundCells[q].index;
261:       ++q;
262:     }
263:     if (globalProcs[p] == size && rank == 0) {
264:       PetscInt d;

266:       for (d = 0; d < ctx->dim; ++d, ++i) a[i] = 0.;
267:       ctx->cells[q] = -1;
268:       ++q;
269:     }
270:   }
271:   PetscCall(VecRestoreArray(ctx->coords, &a));
272: #if 0
273:   PetscCall(PetscFree3(foundCells,foundProcs,globalProcs));
274: #else
275:   PetscCall(PetscFree2(foundProcs, globalProcs));
276:   PetscCall(PetscSFDestroy(&cellSF));
277:   PetscCall(VecDestroy(&pointVec));
278: #endif
279:   if ((void *)globalPointsScalar != (void *)globalPoints) PetscCall(PetscFree(globalPointsScalar));
280:   if (!redundantPoints) PetscCall(PetscFree3(globalPoints, counts, displs));
281:   PetscCall(PetscLayoutDestroy(&layout));
282:   PetscFunctionReturn(PETSC_SUCCESS);
283: }

285: /*@C
286:   DMInterpolationGetCoordinates - Gets a `Vec` with the coordinates of each interpolation point

288:   Collective

290:   Input Parameter:
291: . ctx - the context

293:   Output Parameter:
294: . coordinates - the coordinates of interpolation points

296:   Level: intermediate

298:   Note:
299:   The local vector entries correspond to interpolation points lying on this process, according to the associated `DM`.
300:   This is a borrowed vector that the user should not destroy.

302: .seealso: [](ch_dmbase), `DM`, `DMInterpolationInfo`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`
303: @*/
304: PetscErrorCode DMInterpolationGetCoordinates(DMInterpolationInfo ctx, Vec *coordinates)
305: {
306:   PetscFunctionBegin;
307:   PetscAssertPointer(coordinates, 2);
308:   PetscCheck(ctx->coords, ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup.");
309:   *coordinates = ctx->coords;
310:   PetscFunctionReturn(PETSC_SUCCESS);
311: }

313: /*@C
314:   DMInterpolationGetVector - Gets a `Vec` which can hold all the interpolated field values

316:   Collective

318:   Input Parameter:
319: . ctx - the context

321:   Output Parameter:
322: . v - a vector capable of holding the interpolated field values

324:   Level: intermediate

326:   Note:
327:   This vector should be returned using `DMInterpolationRestoreVector()`.

329: .seealso: [](ch_dmbase), `DM`, `DMInterpolationInfo`, `DMInterpolationRestoreVector()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`
330: @*/
331: PetscErrorCode DMInterpolationGetVector(DMInterpolationInfo ctx, Vec *v)
332: {
333:   PetscFunctionBegin;
334:   PetscAssertPointer(v, 2);
335:   PetscCheck(ctx->coords, ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup.");
336:   PetscCall(VecCreate(ctx->comm, v));
337:   PetscCall(VecSetSizes(*v, ctx->n * ctx->dof, PETSC_DECIDE));
338:   PetscCall(VecSetBlockSize(*v, ctx->dof));
339:   PetscCall(VecSetType(*v, VECSTANDARD));
340:   PetscFunctionReturn(PETSC_SUCCESS);
341: }

343: /*@C
344:   DMInterpolationRestoreVector - Returns a `Vec` which can hold all the interpolated field values

346:   Collective

348:   Input Parameters:
349: + ctx - the context
350: - v   - a vector capable of holding the interpolated field values

352:   Level: intermediate

354: .seealso: [](ch_dmbase), `DM`, `DMInterpolationInfo`, `DMInterpolationGetVector()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`
355: @*/
356: PetscErrorCode DMInterpolationRestoreVector(DMInterpolationInfo ctx, Vec *v)
357: {
358:   PetscFunctionBegin;
359:   PetscAssertPointer(v, 2);
360:   PetscCheck(ctx->coords, ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup.");
361:   PetscCall(VecDestroy(v));
362:   PetscFunctionReturn(PETSC_SUCCESS);
363: }

365: static inline PetscErrorCode DMInterpolate_Segment_Private(DMInterpolationInfo ctx, DM dm, PetscInt p, Vec xLocal, Vec v)
366: {
367:   const PetscInt     c   = ctx->cells[p];
368:   const PetscInt     dof = ctx->dof;
369:   PetscScalar       *x   = NULL;
370:   const PetscScalar *coords;
371:   PetscScalar       *a;
372:   PetscReal          v0, J, invJ, detJ, xir[1];
373:   PetscInt           xSize, comp;

375:   PetscFunctionBegin;
376:   PetscCall(VecGetArrayRead(ctx->coords, &coords));
377:   PetscCall(VecGetArray(v, &a));
378:   PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, &v0, &J, &invJ, &detJ));
379:   PetscCheck(detJ > 0.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %" PetscInt_FMT, (double)detJ, c);
380:   xir[0] = invJ * PetscRealPart(coords[p] - v0);
381:   PetscCall(DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x));
382:   if (2 * dof == xSize) {
383:     for (comp = 0; comp < dof; ++comp) a[p * dof + comp] = x[0 * dof + comp] * (1 - xir[0]) + x[1 * dof + comp] * xir[0];
384:   } else if (dof == xSize) {
385:     for (comp = 0; comp < dof; ++comp) a[p * dof + comp] = x[0 * dof + comp];
386:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Input closure size %" PetscInt_FMT " must be either %" PetscInt_FMT " or %" PetscInt_FMT, xSize, 2 * dof, dof);
387:   PetscCall(DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x));
388:   PetscCall(VecRestoreArray(v, &a));
389:   PetscCall(VecRestoreArrayRead(ctx->coords, &coords));
390:   PetscFunctionReturn(PETSC_SUCCESS);
391: }

393: static inline PetscErrorCode DMInterpolate_Triangle_Private(DMInterpolationInfo ctx, DM dm, PetscInt p, Vec xLocal, Vec v)
394: {
395:   const PetscInt     c = ctx->cells[p];
396:   PetscScalar       *x = NULL;
397:   const PetscScalar *coords;
398:   PetscScalar       *a;
399:   PetscReal         *v0, *J, *invJ, detJ;
400:   PetscReal          xi[4];

402:   PetscFunctionBegin;
403:   PetscCall(PetscMalloc3(ctx->dim, &v0, ctx->dim * ctx->dim, &J, ctx->dim * ctx->dim, &invJ));
404:   PetscCall(VecGetArrayRead(ctx->coords, &coords));
405:   PetscCall(VecGetArray(v, &a));
406:   PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ));
407:   PetscCheck(detJ > 0.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %" PetscInt_FMT, (double)detJ, c);
408:   PetscCall(DMPlexVecGetClosure(dm, NULL, xLocal, c, NULL, &x));
409:   for (PetscInt comp = 0; comp < ctx->dof; ++comp) a[p * ctx->dof + comp] = x[0 * ctx->dof + comp];
410:   for (PetscInt d = 0; d < ctx->dim; ++d) {
411:     xi[d] = 0.0;
412:     for (PetscInt f = 0; f < ctx->dim; ++f) xi[d] += invJ[d * ctx->dim + f] * 0.5 * PetscRealPart(coords[p * ctx->dim + f] - v0[f]);
413:     for (PetscInt comp = 0; comp < ctx->dof; ++comp) a[p * ctx->dof + comp] += PetscRealPart(x[(d + 1) * ctx->dof + comp] - x[0 * ctx->dof + comp]) * xi[d];
414:   }
415:   PetscCall(DMPlexVecRestoreClosure(dm, NULL, xLocal, c, NULL, &x));
416:   PetscCall(VecRestoreArray(v, &a));
417:   PetscCall(VecRestoreArrayRead(ctx->coords, &coords));
418:   PetscCall(PetscFree3(v0, J, invJ));
419:   PetscFunctionReturn(PETSC_SUCCESS);
420: }

422: static inline PetscErrorCode DMInterpolate_Tetrahedron_Private(DMInterpolationInfo ctx, DM dm, PetscInt p, Vec xLocal, Vec v)
423: {
424:   const PetscInt     c        = ctx->cells[p];
425:   const PetscInt     order[3] = {2, 1, 3};
426:   PetscScalar       *x        = NULL;
427:   PetscReal         *v0, *J, *invJ, detJ;
428:   const PetscScalar *coords;
429:   PetscScalar       *a;
430:   PetscReal          xi[4];

432:   PetscFunctionBegin;
433:   PetscCall(PetscMalloc3(ctx->dim, &v0, ctx->dim * ctx->dim, &J, ctx->dim * ctx->dim, &invJ));
434:   PetscCall(VecGetArrayRead(ctx->coords, &coords));
435:   PetscCall(VecGetArray(v, &a));
436:   PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ));
437:   PetscCheck(detJ > 0.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %" PetscInt_FMT, (double)detJ, c);
438:   PetscCall(DMPlexVecGetClosure(dm, NULL, xLocal, c, NULL, &x));
439:   for (PetscInt comp = 0; comp < ctx->dof; ++comp) a[p * ctx->dof + comp] = x[0 * ctx->dof + comp];
440:   for (PetscInt d = 0; d < ctx->dim; ++d) {
441:     xi[d] = 0.0;
442:     for (PetscInt f = 0; f < ctx->dim; ++f) xi[d] += invJ[d * ctx->dim + f] * 0.5 * PetscRealPart(coords[p * ctx->dim + f] - v0[f]);
443:     for (PetscInt comp = 0; comp < ctx->dof; ++comp) a[p * ctx->dof + comp] += PetscRealPart(x[order[d] * ctx->dof + comp] - x[0 * ctx->dof + comp]) * xi[d];
444:   }
445:   PetscCall(DMPlexVecRestoreClosure(dm, NULL, xLocal, c, NULL, &x));
446:   PetscCall(VecRestoreArray(v, &a));
447:   PetscCall(VecRestoreArrayRead(ctx->coords, &coords));
448:   PetscCall(PetscFree3(v0, J, invJ));
449:   PetscFunctionReturn(PETSC_SUCCESS);
450: }

452: static inline PetscErrorCode QuadMap_Private(SNES snes, Vec Xref, Vec Xreal, void *ctx)
453: {
454:   const PetscScalar *vertices = (const PetscScalar *)ctx;
455:   const PetscScalar  x0       = vertices[0];
456:   const PetscScalar  y0       = vertices[1];
457:   const PetscScalar  x1       = vertices[2];
458:   const PetscScalar  y1       = vertices[3];
459:   const PetscScalar  x2       = vertices[4];
460:   const PetscScalar  y2       = vertices[5];
461:   const PetscScalar  x3       = vertices[6];
462:   const PetscScalar  y3       = vertices[7];
463:   const PetscScalar  f_1      = x1 - x0;
464:   const PetscScalar  g_1      = y1 - y0;
465:   const PetscScalar  f_3      = x3 - x0;
466:   const PetscScalar  g_3      = y3 - y0;
467:   const PetscScalar  f_01     = x2 - x1 - x3 + x0;
468:   const PetscScalar  g_01     = y2 - y1 - y3 + y0;
469:   const PetscScalar *ref;
470:   PetscScalar       *real;

472:   PetscFunctionBegin;
473:   PetscCall(VecGetArrayRead(Xref, &ref));
474:   PetscCall(VecGetArray(Xreal, &real));
475:   {
476:     const PetscScalar p0 = ref[0];
477:     const PetscScalar p1 = ref[1];

479:     real[0] = x0 + f_1 * p0 + f_3 * p1 + f_01 * p0 * p1;
480:     real[1] = y0 + g_1 * p0 + g_3 * p1 + g_01 * p0 * p1;
481:   }
482:   PetscCall(PetscLogFlops(28));
483:   PetscCall(VecRestoreArrayRead(Xref, &ref));
484:   PetscCall(VecRestoreArray(Xreal, &real));
485:   PetscFunctionReturn(PETSC_SUCCESS);
486: }

488: #include <petsc/private/dmimpl.h>
489: static inline PetscErrorCode QuadJacobian_Private(SNES snes, Vec Xref, Mat J, Mat M, void *ctx)
490: {
491:   const PetscScalar *vertices = (const PetscScalar *)ctx;
492:   const PetscScalar  x0       = vertices[0];
493:   const PetscScalar  y0       = vertices[1];
494:   const PetscScalar  x1       = vertices[2];
495:   const PetscScalar  y1       = vertices[3];
496:   const PetscScalar  x2       = vertices[4];
497:   const PetscScalar  y2       = vertices[5];
498:   const PetscScalar  x3       = vertices[6];
499:   const PetscScalar  y3       = vertices[7];
500:   const PetscScalar  f_01     = x2 - x1 - x3 + x0;
501:   const PetscScalar  g_01     = y2 - y1 - y3 + y0;
502:   const PetscScalar *ref;

504:   PetscFunctionBegin;
505:   PetscCall(VecGetArrayRead(Xref, &ref));
506:   {
507:     const PetscScalar x       = ref[0];
508:     const PetscScalar y       = ref[1];
509:     const PetscInt    rows[2] = {0, 1};
510:     PetscScalar       values[4];

512:     values[0] = (x1 - x0 + f_01 * y) * 0.5;
513:     values[1] = (x3 - x0 + f_01 * x) * 0.5;
514:     values[2] = (y1 - y0 + g_01 * y) * 0.5;
515:     values[3] = (y3 - y0 + g_01 * x) * 0.5;
516:     PetscCall(MatSetValues(J, 2, rows, 2, rows, values, INSERT_VALUES));
517:   }
518:   PetscCall(PetscLogFlops(30));
519:   PetscCall(VecRestoreArrayRead(Xref, &ref));
520:   PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
521:   PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
522:   PetscFunctionReturn(PETSC_SUCCESS);
523: }

525: static inline PetscErrorCode DMInterpolate_Quad_Private(DMInterpolationInfo ctx, DM dm, PetscInt p, Vec xLocal, Vec v)
526: {
527:   const PetscInt     c   = ctx->cells[p];
528:   PetscFE            fem = NULL;
529:   DM                 dmCoord;
530:   SNES               snes;
531:   KSP                ksp;
532:   PC                 pc;
533:   Vec                coordsLocal, r, ref, real;
534:   Mat                J;
535:   PetscTabulation    T = NULL;
536:   const PetscScalar *coords;
537:   PetscScalar       *a;
538:   PetscReal          xir[2] = {0., 0.};
539:   PetscInt           Nf;
540:   const PetscInt     dof = ctx->dof;
541:   PetscScalar       *x = NULL, *vertices = NULL;
542:   PetscScalar       *xi;
543:   PetscInt           coordSize, xSize;

545:   PetscFunctionBegin;
546:   PetscCall(DMGetNumFields(dm, &Nf));
547:   if (Nf) {
548:     PetscObject  obj;
549:     PetscClassId id;

551:     PetscCall(DMGetField(dm, 0, NULL, &obj));
552:     PetscCall(PetscObjectGetClassId(obj, &id));
553:     if (id == PETSCFE_CLASSID) {
554:       fem = (PetscFE)obj;
555:       PetscCall(PetscFECreateTabulation(fem, 1, 1, xir, 0, &T));
556:     }
557:   }
558:   PetscCall(DMGetCoordinatesLocal(dm, &coordsLocal));
559:   PetscCall(DMGetCoordinateDM(dm, &dmCoord));
560:   PetscCall(SNESCreate(PETSC_COMM_SELF, &snes));
561:   PetscCall(SNESSetOptionsPrefix(snes, "quad_interp_"));
562:   PetscCall(VecCreate(PETSC_COMM_SELF, &r));
563:   PetscCall(VecSetSizes(r, 2, 2));
564:   PetscCall(VecSetType(r, dm->vectype));
565:   PetscCall(VecDuplicate(r, &ref));
566:   PetscCall(VecDuplicate(r, &real));
567:   PetscCall(MatCreate(PETSC_COMM_SELF, &J));
568:   PetscCall(MatSetSizes(J, 2, 2, 2, 2));
569:   PetscCall(MatSetType(J, MATSEQDENSE));
570:   PetscCall(MatSetUp(J));
571:   PetscCall(SNESSetFunction(snes, r, QuadMap_Private, NULL));
572:   PetscCall(SNESSetJacobian(snes, J, J, QuadJacobian_Private, NULL));
573:   PetscCall(SNESGetKSP(snes, &ksp));
574:   PetscCall(KSPGetPC(ksp, &pc));
575:   PetscCall(PCSetType(pc, PCLU));
576:   PetscCall(SNESSetFromOptions(snes));

578:   PetscCall(VecGetArrayRead(ctx->coords, &coords));
579:   PetscCall(VecGetArray(v, &a));
580:   PetscCall(DMPlexVecGetClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices));
581:   PetscCheck(4 * 2 == coordSize, ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid closure size %" PetscInt_FMT " should be %d", coordSize, 4 * 2);
582:   PetscCall(DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x));
583:   PetscCall(SNESSetFunction(snes, NULL, NULL, vertices));
584:   PetscCall(SNESSetJacobian(snes, NULL, NULL, NULL, vertices));
585:   PetscCall(VecGetArray(real, &xi));
586:   xi[0] = coords[p * ctx->dim + 0];
587:   xi[1] = coords[p * ctx->dim + 1];
588:   PetscCall(VecRestoreArray(real, &xi));
589:   PetscCall(SNESSolve(snes, real, ref));
590:   PetscCall(VecGetArray(ref, &xi));
591:   xir[0] = PetscRealPart(xi[0]);
592:   xir[1] = PetscRealPart(xi[1]);
593:   if (4 * dof == xSize) {
594:     for (PetscInt comp = 0; comp < dof; ++comp) a[p * dof + comp] = x[0 * dof + comp] * (1 - xir[0]) * (1 - xir[1]) + x[1 * dof + comp] * xir[0] * (1 - xir[1]) + x[2 * dof + comp] * xir[0] * xir[1] + x[3 * dof + comp] * (1 - xir[0]) * xir[1];
595:   } else if (dof == xSize) {
596:     for (PetscInt comp = 0; comp < dof; ++comp) a[p * dof + comp] = x[0 * dof + comp];
597:   } else {
598:     PetscCheck(fem, ctx->comm, PETSC_ERR_ARG_WRONG, "Cannot have a higher order interpolant if the discretization is not PetscFE");
599:     xir[0] = 2.0 * xir[0] - 1.0;
600:     xir[1] = 2.0 * xir[1] - 1.0;
601:     PetscCall(PetscFEComputeTabulation(fem, 1, xir, 0, T));
602:     for (PetscInt comp = 0; comp < dof; ++comp) {
603:       a[p * dof + comp] = 0.0;
604:       for (PetscInt d = 0; d < xSize / dof; ++d) a[p * dof + comp] += x[d * dof + comp] * T->T[0][d * dof + comp];
605:     }
606:   }
607:   PetscCall(VecRestoreArray(ref, &xi));
608:   PetscCall(DMPlexVecRestoreClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices));
609:   PetscCall(DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x));
610:   PetscCall(PetscTabulationDestroy(&T));
611:   PetscCall(VecRestoreArray(v, &a));
612:   PetscCall(VecRestoreArrayRead(ctx->coords, &coords));

614:   PetscCall(SNESDestroy(&snes));
615:   PetscCall(VecDestroy(&r));
616:   PetscCall(VecDestroy(&ref));
617:   PetscCall(VecDestroy(&real));
618:   PetscCall(MatDestroy(&J));
619:   PetscFunctionReturn(PETSC_SUCCESS);
620: }

622: static inline PetscErrorCode HexMap_Private(SNES snes, Vec Xref, Vec Xreal, void *ctx)
623: {
624:   const PetscScalar *vertices = (const PetscScalar *)ctx;
625:   const PetscScalar  x0       = vertices[0];
626:   const PetscScalar  y0       = vertices[1];
627:   const PetscScalar  z0       = vertices[2];
628:   const PetscScalar  x1       = vertices[9];
629:   const PetscScalar  y1       = vertices[10];
630:   const PetscScalar  z1       = vertices[11];
631:   const PetscScalar  x2       = vertices[6];
632:   const PetscScalar  y2       = vertices[7];
633:   const PetscScalar  z2       = vertices[8];
634:   const PetscScalar  x3       = vertices[3];
635:   const PetscScalar  y3       = vertices[4];
636:   const PetscScalar  z3       = vertices[5];
637:   const PetscScalar  x4       = vertices[12];
638:   const PetscScalar  y4       = vertices[13];
639:   const PetscScalar  z4       = vertices[14];
640:   const PetscScalar  x5       = vertices[15];
641:   const PetscScalar  y5       = vertices[16];
642:   const PetscScalar  z5       = vertices[17];
643:   const PetscScalar  x6       = vertices[18];
644:   const PetscScalar  y6       = vertices[19];
645:   const PetscScalar  z6       = vertices[20];
646:   const PetscScalar  x7       = vertices[21];
647:   const PetscScalar  y7       = vertices[22];
648:   const PetscScalar  z7       = vertices[23];
649:   const PetscScalar  f_1      = x1 - x0;
650:   const PetscScalar  g_1      = y1 - y0;
651:   const PetscScalar  h_1      = z1 - z0;
652:   const PetscScalar  f_3      = x3 - x0;
653:   const PetscScalar  g_3      = y3 - y0;
654:   const PetscScalar  h_3      = z3 - z0;
655:   const PetscScalar  f_4      = x4 - x0;
656:   const PetscScalar  g_4      = y4 - y0;
657:   const PetscScalar  h_4      = z4 - z0;
658:   const PetscScalar  f_01     = x2 - x1 - x3 + x0;
659:   const PetscScalar  g_01     = y2 - y1 - y3 + y0;
660:   const PetscScalar  h_01     = z2 - z1 - z3 + z0;
661:   const PetscScalar  f_12     = x7 - x3 - x4 + x0;
662:   const PetscScalar  g_12     = y7 - y3 - y4 + y0;
663:   const PetscScalar  h_12     = z7 - z3 - z4 + z0;
664:   const PetscScalar  f_02     = x5 - x1 - x4 + x0;
665:   const PetscScalar  g_02     = y5 - y1 - y4 + y0;
666:   const PetscScalar  h_02     = z5 - z1 - z4 + z0;
667:   const PetscScalar  f_012    = x6 - x0 + x1 - x2 + x3 + x4 - x5 - x7;
668:   const PetscScalar  g_012    = y6 - y0 + y1 - y2 + y3 + y4 - y5 - y7;
669:   const PetscScalar  h_012    = z6 - z0 + z1 - z2 + z3 + z4 - z5 - z7;
670:   const PetscScalar *ref;
671:   PetscScalar       *real;

673:   PetscFunctionBegin;
674:   PetscCall(VecGetArrayRead(Xref, &ref));
675:   PetscCall(VecGetArray(Xreal, &real));
676:   {
677:     const PetscScalar p0 = ref[0];
678:     const PetscScalar p1 = ref[1];
679:     const PetscScalar p2 = ref[2];

681:     real[0] = x0 + f_1 * p0 + f_3 * p1 + f_4 * p2 + f_01 * p0 * p1 + f_12 * p1 * p2 + f_02 * p0 * p2 + f_012 * p0 * p1 * p2;
682:     real[1] = y0 + g_1 * p0 + g_3 * p1 + g_4 * p2 + g_01 * p0 * p1 + g_01 * p0 * p1 + g_12 * p1 * p2 + g_02 * p0 * p2 + g_012 * p0 * p1 * p2;
683:     real[2] = z0 + h_1 * p0 + h_3 * p1 + h_4 * p2 + h_01 * p0 * p1 + h_01 * p0 * p1 + h_12 * p1 * p2 + h_02 * p0 * p2 + h_012 * p0 * p1 * p2;
684:   }
685:   PetscCall(PetscLogFlops(114));
686:   PetscCall(VecRestoreArrayRead(Xref, &ref));
687:   PetscCall(VecRestoreArray(Xreal, &real));
688:   PetscFunctionReturn(PETSC_SUCCESS);
689: }

691: static inline PetscErrorCode HexJacobian_Private(SNES snes, Vec Xref, Mat J, Mat M, void *ctx)
692: {
693:   const PetscScalar *vertices = (const PetscScalar *)ctx;
694:   const PetscScalar  x0       = vertices[0];
695:   const PetscScalar  y0       = vertices[1];
696:   const PetscScalar  z0       = vertices[2];
697:   const PetscScalar  x1       = vertices[9];
698:   const PetscScalar  y1       = vertices[10];
699:   const PetscScalar  z1       = vertices[11];
700:   const PetscScalar  x2       = vertices[6];
701:   const PetscScalar  y2       = vertices[7];
702:   const PetscScalar  z2       = vertices[8];
703:   const PetscScalar  x3       = vertices[3];
704:   const PetscScalar  y3       = vertices[4];
705:   const PetscScalar  z3       = vertices[5];
706:   const PetscScalar  x4       = vertices[12];
707:   const PetscScalar  y4       = vertices[13];
708:   const PetscScalar  z4       = vertices[14];
709:   const PetscScalar  x5       = vertices[15];
710:   const PetscScalar  y5       = vertices[16];
711:   const PetscScalar  z5       = vertices[17];
712:   const PetscScalar  x6       = vertices[18];
713:   const PetscScalar  y6       = vertices[19];
714:   const PetscScalar  z6       = vertices[20];
715:   const PetscScalar  x7       = vertices[21];
716:   const PetscScalar  y7       = vertices[22];
717:   const PetscScalar  z7       = vertices[23];
718:   const PetscScalar  f_xy     = x2 - x1 - x3 + x0;
719:   const PetscScalar  g_xy     = y2 - y1 - y3 + y0;
720:   const PetscScalar  h_xy     = z2 - z1 - z3 + z0;
721:   const PetscScalar  f_yz     = x7 - x3 - x4 + x0;
722:   const PetscScalar  g_yz     = y7 - y3 - y4 + y0;
723:   const PetscScalar  h_yz     = z7 - z3 - z4 + z0;
724:   const PetscScalar  f_xz     = x5 - x1 - x4 + x0;
725:   const PetscScalar  g_xz     = y5 - y1 - y4 + y0;
726:   const PetscScalar  h_xz     = z5 - z1 - z4 + z0;
727:   const PetscScalar  f_xyz    = x6 - x0 + x1 - x2 + x3 + x4 - x5 - x7;
728:   const PetscScalar  g_xyz    = y6 - y0 + y1 - y2 + y3 + y4 - y5 - y7;
729:   const PetscScalar  h_xyz    = z6 - z0 + z1 - z2 + z3 + z4 - z5 - z7;
730:   const PetscScalar *ref;

732:   PetscFunctionBegin;
733:   PetscCall(VecGetArrayRead(Xref, &ref));
734:   {
735:     const PetscScalar x       = ref[0];
736:     const PetscScalar y       = ref[1];
737:     const PetscScalar z       = ref[2];
738:     const PetscInt    rows[3] = {0, 1, 2};
739:     PetscScalar       values[9];

741:     values[0] = (x1 - x0 + f_xy * y + f_xz * z + f_xyz * y * z) / 2.0;
742:     values[1] = (x3 - x0 + f_xy * x + f_yz * z + f_xyz * x * z) / 2.0;
743:     values[2] = (x4 - x0 + f_yz * y + f_xz * x + f_xyz * x * y) / 2.0;
744:     values[3] = (y1 - y0 + g_xy * y + g_xz * z + g_xyz * y * z) / 2.0;
745:     values[4] = (y3 - y0 + g_xy * x + g_yz * z + g_xyz * x * z) / 2.0;
746:     values[5] = (y4 - y0 + g_yz * y + g_xz * x + g_xyz * x * y) / 2.0;
747:     values[6] = (z1 - z0 + h_xy * y + h_xz * z + h_xyz * y * z) / 2.0;
748:     values[7] = (z3 - z0 + h_xy * x + h_yz * z + h_xyz * x * z) / 2.0;
749:     values[8] = (z4 - z0 + h_yz * y + h_xz * x + h_xyz * x * y) / 2.0;

751:     PetscCall(MatSetValues(J, 3, rows, 3, rows, values, INSERT_VALUES));
752:   }
753:   PetscCall(PetscLogFlops(152));
754:   PetscCall(VecRestoreArrayRead(Xref, &ref));
755:   PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
756:   PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
757:   PetscFunctionReturn(PETSC_SUCCESS);
758: }

760: static inline PetscErrorCode DMInterpolate_Hex_Private(DMInterpolationInfo ctx, DM dm, PetscInt p, Vec xLocal, Vec v)
761: {
762:   const PetscInt     c = ctx->cells[p];
763:   DM                 dmCoord;
764:   SNES               snes;
765:   KSP                ksp;
766:   PC                 pc;
767:   Vec                coordsLocal, r, ref, real;
768:   Mat                J;
769:   const PetscScalar *coords;
770:   PetscScalar       *a;
771:   PetscScalar       *x = NULL, *vertices = NULL;
772:   PetscScalar       *xi;
773:   PetscReal          xir[3];
774:   PetscInt           coordSize, xSize;

776:   PetscFunctionBegin;
777:   PetscCall(DMGetCoordinatesLocal(dm, &coordsLocal));
778:   PetscCall(DMGetCoordinateDM(dm, &dmCoord));
779:   PetscCall(SNESCreate(PETSC_COMM_SELF, &snes));
780:   PetscCall(SNESSetOptionsPrefix(snes, "hex_interp_"));
781:   PetscCall(VecCreate(PETSC_COMM_SELF, &r));
782:   PetscCall(VecSetSizes(r, 3, 3));
783:   PetscCall(VecSetType(r, dm->vectype));
784:   PetscCall(VecDuplicate(r, &ref));
785:   PetscCall(VecDuplicate(r, &real));
786:   PetscCall(MatCreate(PETSC_COMM_SELF, &J));
787:   PetscCall(MatSetSizes(J, 3, 3, 3, 3));
788:   PetscCall(MatSetType(J, MATSEQDENSE));
789:   PetscCall(MatSetUp(J));
790:   PetscCall(SNESSetFunction(snes, r, HexMap_Private, NULL));
791:   PetscCall(SNESSetJacobian(snes, J, J, HexJacobian_Private, NULL));
792:   PetscCall(SNESGetKSP(snes, &ksp));
793:   PetscCall(KSPGetPC(ksp, &pc));
794:   PetscCall(PCSetType(pc, PCLU));
795:   PetscCall(SNESSetFromOptions(snes));

797:   PetscCall(VecGetArrayRead(ctx->coords, &coords));
798:   PetscCall(VecGetArray(v, &a));
799:   PetscCall(DMPlexVecGetClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices));
800:   PetscCheck(8 * 3 == coordSize, ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid coordinate closure size %" PetscInt_FMT " should be %d", coordSize, 8 * 3);
801:   PetscCall(DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x));
802:   PetscCheck((8 * ctx->dof == xSize) || (ctx->dof == xSize), ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid input closure size %" PetscInt_FMT " should be %" PetscInt_FMT " or %" PetscInt_FMT, xSize, 8 * ctx->dof, ctx->dof);
803:   PetscCall(SNESSetFunction(snes, NULL, NULL, vertices));
804:   PetscCall(SNESSetJacobian(snes, NULL, NULL, NULL, vertices));
805:   PetscCall(VecGetArray(real, &xi));
806:   xi[0] = coords[p * ctx->dim + 0];
807:   xi[1] = coords[p * ctx->dim + 1];
808:   xi[2] = coords[p * ctx->dim + 2];
809:   PetscCall(VecRestoreArray(real, &xi));
810:   PetscCall(SNESSolve(snes, real, ref));
811:   PetscCall(VecGetArray(ref, &xi));
812:   xir[0] = PetscRealPart(xi[0]);
813:   xir[1] = PetscRealPart(xi[1]);
814:   xir[2] = PetscRealPart(xi[2]);
815:   if (8 * ctx->dof == xSize) {
816:     for (PetscInt comp = 0; comp < ctx->dof; ++comp) {
817:       a[p * ctx->dof + comp] = x[0 * ctx->dof + comp] * (1 - xir[0]) * (1 - xir[1]) * (1 - xir[2]) + x[3 * ctx->dof + comp] * xir[0] * (1 - xir[1]) * (1 - xir[2]) + x[2 * ctx->dof + comp] * xir[0] * xir[1] * (1 - xir[2]) + x[1 * ctx->dof + comp] * (1 - xir[0]) * xir[1] * (1 - xir[2]) +
818:                                x[4 * ctx->dof + comp] * (1 - xir[0]) * (1 - xir[1]) * xir[2] + x[5 * ctx->dof + comp] * xir[0] * (1 - xir[1]) * xir[2] + x[6 * ctx->dof + comp] * xir[0] * xir[1] * xir[2] + x[7 * ctx->dof + comp] * (1 - xir[0]) * xir[1] * xir[2];
819:     }
820:   } else {
821:     for (PetscInt comp = 0; comp < ctx->dof; ++comp) a[p * ctx->dof + comp] = x[0 * ctx->dof + comp];
822:   }
823:   PetscCall(VecRestoreArray(ref, &xi));
824:   PetscCall(DMPlexVecRestoreClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices));
825:   PetscCall(DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x));
826:   PetscCall(VecRestoreArray(v, &a));
827:   PetscCall(VecRestoreArrayRead(ctx->coords, &coords));

829:   PetscCall(SNESDestroy(&snes));
830:   PetscCall(VecDestroy(&r));
831:   PetscCall(VecDestroy(&ref));
832:   PetscCall(VecDestroy(&real));
833:   PetscCall(MatDestroy(&J));
834:   PetscFunctionReturn(PETSC_SUCCESS);
835: }

837: /*@C
838:   DMInterpolationEvaluate - Using the input from `dm` and `x`, calculates interpolated field values at the interpolation points.

840:   Input Parameters:
841: + ctx - The `DMInterpolationInfo` context obtained with `DMInterpolationCreate()`
842: . dm  - The `DM`
843: - x   - The local vector containing the field to be interpolated, can be created with `DMCreateGlobalVector()`

845:   Output Parameter:
846: . v - The vector containing the interpolated values, obtained with `DMInterpolationGetVector()`

848:   Level: beginner

850: .seealso: [](ch_dmbase), `DM`, `DMInterpolationInfo`, `DMInterpolationGetVector()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`, `DMInterpolationGetCoordinates()`
851: @*/
852: PetscErrorCode DMInterpolationEvaluate(DMInterpolationInfo ctx, DM dm, Vec x, Vec v)
853: {
854:   PetscDS   ds;
855:   PetscInt  n, p, Nf, field;
856:   PetscBool useDS = PETSC_FALSE;

858:   PetscFunctionBegin;
862:   PetscCall(VecGetLocalSize(v, &n));
863:   PetscCheck(n == ctx->n * ctx->dof, ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid input vector size %" PetscInt_FMT " should be %" PetscInt_FMT, n, ctx->n * ctx->dof);
864:   if (!n) PetscFunctionReturn(PETSC_SUCCESS);
865:   PetscCall(DMGetDS(dm, &ds));
866:   if (ds) {
867:     useDS = PETSC_TRUE;
868:     PetscCall(PetscDSGetNumFields(ds, &Nf));
869:     for (field = 0; field < Nf; ++field) {
870:       PetscObject  obj;
871:       PetscClassId id;

873:       PetscCall(PetscDSGetDiscretization(ds, field, &obj));
874:       PetscCall(PetscObjectGetClassId(obj, &id));
875:       if (id != PETSCFE_CLASSID && id != PETSCFV_CLASSID) {
876:         useDS = PETSC_FALSE;
877:         break;
878:       }
879:     }
880:   }
881:   if (useDS) {
882:     const PetscScalar *coords;
883:     PetscScalar       *interpolant;
884:     PetscInt           cdim, d;

886:     PetscCall(DMGetCoordinateDim(dm, &cdim));
887:     PetscCall(VecGetArrayRead(ctx->coords, &coords));
888:     PetscCall(VecGetArrayWrite(v, &interpolant));
889:     for (p = 0; p < ctx->n; ++p) {
890:       PetscReal    pcoords[3], xi[3];
891:       PetscScalar *xa   = NULL;
892:       PetscInt     coff = 0, foff = 0, clSize;

894:       if (ctx->cells[p] < 0) continue;
895:       for (d = 0; d < cdim; ++d) pcoords[d] = PetscRealPart(coords[p * cdim + d]);
896:       PetscCall(DMPlexCoordinatesToReference(dm, ctx->cells[p], 1, pcoords, xi));
897:       PetscCall(DMPlexVecGetClosure(dm, NULL, x, ctx->cells[p], &clSize, &xa));
898:       for (field = 0; field < Nf; ++field) {
899:         PetscTabulation T;
900:         PetscObject     obj;
901:         PetscClassId    id;

903:         PetscCall(PetscDSGetDiscretization(ds, field, &obj));
904:         PetscCall(PetscObjectGetClassId(obj, &id));
905:         if (id == PETSCFE_CLASSID) {
906:           PetscFE fe = (PetscFE)obj;

908:           PetscCall(PetscFECreateTabulation(fe, 1, 1, xi, 0, &T));
909:           {
910:             const PetscReal *basis = T->T[0];
911:             const PetscInt   Nb    = T->Nb;
912:             const PetscInt   Nc    = T->Nc;

914:             for (PetscInt fc = 0; fc < Nc; ++fc) {
915:               interpolant[p * ctx->dof + coff + fc] = 0.0;
916:               for (PetscInt f = 0; f < Nb; ++f) interpolant[p * ctx->dof + coff + fc] += xa[foff + f] * basis[(0 * Nb + f) * Nc + fc];
917:             }
918:             coff += Nc;
919:             foff += Nb;
920:           }
921:           PetscCall(PetscTabulationDestroy(&T));
922:         } else if (id == PETSCFV_CLASSID) {
923:           PetscFV  fv = (PetscFV)obj;
924:           PetscInt Nc;

926:           // TODO Could use reconstruction if available
927:           PetscCall(PetscFVGetNumComponents(fv, &Nc));
928:           for (PetscInt fc = 0; fc < Nc; ++fc) interpolant[p * ctx->dof + coff + fc] = xa[foff + fc];
929:           coff += Nc;
930:           foff += Nc;
931:         }
932:       }
933:       PetscCall(DMPlexVecRestoreClosure(dm, NULL, x, ctx->cells[p], &clSize, &xa));
934:       PetscCheck(coff == ctx->dof, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Total components %" PetscInt_FMT " != %" PetscInt_FMT " dof specified for interpolation", coff, ctx->dof);
935:       PetscCheck(foff == clSize, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Total FE/FV space size %" PetscInt_FMT " != %" PetscInt_FMT " closure size", foff, clSize);
936:     }
937:     PetscCall(VecRestoreArrayRead(ctx->coords, &coords));
938:     PetscCall(VecRestoreArrayWrite(v, &interpolant));
939:   } else {
940:     for (PetscInt p = 0; p < ctx->n; ++p) {
941:       const PetscInt cell = ctx->cells[p];
942:       DMPolytopeType ct;

944:       PetscCall(DMPlexGetCellType(dm, cell, &ct));
945:       switch (ct) {
946:       case DM_POLYTOPE_SEGMENT:
947:         PetscCall(DMInterpolate_Segment_Private(ctx, dm, p, x, v));
948:         break;
949:       case DM_POLYTOPE_TRIANGLE:
950:         PetscCall(DMInterpolate_Triangle_Private(ctx, dm, p, x, v));
951:         break;
952:       case DM_POLYTOPE_QUADRILATERAL:
953:         PetscCall(DMInterpolate_Quad_Private(ctx, dm, p, x, v));
954:         break;
955:       case DM_POLYTOPE_TETRAHEDRON:
956:         PetscCall(DMInterpolate_Tetrahedron_Private(ctx, dm, p, x, v));
957:         break;
958:       case DM_POLYTOPE_HEXAHEDRON:
959:         PetscCall(DMInterpolate_Hex_Private(ctx, dm, cell, x, v));
960:         break;
961:       default:
962:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for cell type %s", DMPolytopeTypes[PetscMax(0, PetscMin(ct, DM_NUM_POLYTOPES))]);
963:       }
964:     }
965:   }
966:   PetscFunctionReturn(PETSC_SUCCESS);
967: }

969: /*@C
970:   DMInterpolationDestroy - Destroys a `DMInterpolationInfo` context

972:   Collective

974:   Input Parameter:
975: . ctx - the context

977:   Level: beginner

979: .seealso: [](ch_dmbase), `DM`, `DMInterpolationInfo`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`
980: @*/
981: PetscErrorCode DMInterpolationDestroy(DMInterpolationInfo *ctx)
982: {
983:   PetscFunctionBegin;
984:   PetscAssertPointer(ctx, 1);
985:   PetscCall(VecDestroy(&(*ctx)->coords));
986:   PetscCall(PetscFree((*ctx)->points));
987:   PetscCall(PetscFree((*ctx)->cells));
988:   PetscCall(PetscFree(*ctx));
989:   *ctx = NULL;
990:   PetscFunctionReturn(PETSC_SUCCESS);
991: }