xref: /petsc/src/snes/utils/dm/dminterpolatesnes.c (revision c6f61ee217dbd23bd0792f1ef4cfacda5212558b) 
1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/
2 #include <petsc/private/snesimpl.h>   /*I "petscsnes.h"   I*/
3 #include <petscds.h>
4 #include <petsc/private/petscimpl.h>
5 #include <petsc/private/petscfeimpl.h>
6 
7 /*@C
8   DMInterpolationCreate - Creates a `DMInterpolationInfo` context
9 
10   Collective
11 
12   Input Parameter:
13 . comm - the communicator
14 
15   Output Parameter:
16 . ctx - the context
17 
18   Level: beginner
19 
20 .seealso: [](ch_snes), `DM`, `DMInterpolationInfo`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationDestroy()`
21 @*/
22 PetscErrorCode DMInterpolationCreate(MPI_Comm comm, DMInterpolationInfo *ctx)
23 {
24   PetscFunctionBegin;
25   PetscAssertPointer(ctx, 2);
26   PetscCall(PetscNew(ctx));
27 
28   (*ctx)->comm   = comm;
29   (*ctx)->dim    = -1;
30   (*ctx)->nInput = 0;
31   (*ctx)->points = NULL;
32   (*ctx)->cells  = NULL;
33   (*ctx)->n      = -1;
34   (*ctx)->coords = NULL;
35   PetscFunctionReturn(PETSC_SUCCESS);
36 }
37 
38 /*@C
39   DMInterpolationSetDim - Sets the spatial dimension for the interpolation context
40 
41   Not Collective
42 
43   Input Parameters:
44 + ctx - the context
45 - dim - the spatial dimension
46 
47   Level: intermediate
48 
49 .seealso: [](ch_snes), `DM`, `DMInterpolationInfo`, `DMInterpolationGetDim()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`
50 @*/
51 PetscErrorCode DMInterpolationSetDim(DMInterpolationInfo ctx, PetscInt dim)
52 {
53   PetscFunctionBegin;
54   PetscCheck(!(dim < 1) && !(dim > 3), ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension for points: %" PetscInt_FMT, dim);
55   ctx->dim = dim;
56   PetscFunctionReturn(PETSC_SUCCESS);
57 }
58 
59 /*@C
60   DMInterpolationGetDim - Gets the spatial dimension for the interpolation context
61 
62   Not Collective
63 
64   Input Parameter:
65 . ctx - the context
66 
67   Output Parameter:
68 . dim - the spatial dimension
69 
70   Level: intermediate
71 
72 .seealso: [](ch_snes), `DM`, `DMInterpolationInfo`, `DMInterpolationSetDim()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`
73 @*/
74 PetscErrorCode DMInterpolationGetDim(DMInterpolationInfo ctx, PetscInt *dim)
75 {
76   PetscFunctionBegin;
77   PetscAssertPointer(dim, 2);
78   *dim = ctx->dim;
79   PetscFunctionReturn(PETSC_SUCCESS);
80 }
81 
82 /*@C
83   DMInterpolationSetDof - Sets the number of fields interpolated at a point for the interpolation context
84 
85   Not Collective
86 
87   Input Parameters:
88 + ctx - the context
89 - dof - the number of fields
90 
91   Level: intermediate
92 
93 .seealso: [](ch_snes), `DM`, `DMInterpolationInfo`, `DMInterpolationGetDof()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`
94 @*/
95 PetscErrorCode DMInterpolationSetDof(DMInterpolationInfo ctx, PetscInt dof)
96 {
97   PetscFunctionBegin;
98   PetscCheck(dof >= 1, ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid number of components: %" PetscInt_FMT, dof);
99   ctx->dof = dof;
100   PetscFunctionReturn(PETSC_SUCCESS);
101 }
102 
103 /*@C
104   DMInterpolationGetDof - Gets the number of fields interpolated at a point for the interpolation context
105 
106   Not Collective
107 
108   Input Parameter:
109 . ctx - the context
110 
111   Output Parameter:
112 . dof - the number of fields
113 
114   Level: intermediate
115 
116 .seealso: [](ch_snes), `DM`, `DMInterpolationInfo`, `DMInterpolationSetDof()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`
117 @*/
118 PetscErrorCode DMInterpolationGetDof(DMInterpolationInfo ctx, PetscInt *dof)
119 {
120   PetscFunctionBegin;
121   PetscAssertPointer(dof, 2);
122   *dof = ctx->dof;
123   PetscFunctionReturn(PETSC_SUCCESS);
124 }
125 
126 /*@C
127   DMInterpolationAddPoints - Add points at which we will interpolate the fields
128 
129   Not Collective
130 
131   Input Parameters:
132 + ctx    - the context
133 . n      - the number of points
134 - points - the coordinates for each point, an array of size n * dim
135 
136   Level: intermediate
137 
138   Note:
139   The coordinate information is copied.
140 
141 .seealso: [](ch_snes), `DM`, `DMInterpolationInfo`, `DMInterpolationSetDim()`, `DMInterpolationEvaluate()`, `DMInterpolationCreate()`
142 @*/
143 PetscErrorCode DMInterpolationAddPoints(DMInterpolationInfo ctx, PetscInt n, PetscReal points[])
144 {
145   PetscFunctionBegin;
146   PetscCheck(ctx->dim >= 0, ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The spatial dimension has not been set");
147   PetscCheck(!ctx->points, ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "Cannot add points multiple times yet");
148   ctx->nInput = n;
149 
150   PetscCall(PetscMalloc1(n * ctx->dim, &ctx->points));
151   PetscCall(PetscArraycpy(ctx->points, points, n * ctx->dim));
152   PetscFunctionReturn(PETSC_SUCCESS);
153 }
154 
155 /*@C
156   DMInterpolationSetUp - Compute spatial indices for point location during interpolation
157 
158   Collective
159 
160   Input Parameters:
161 + ctx                 - the context
162 . dm                  - the `DM` for the function space used for interpolation
163 . redundantPoints     - If `PETSC_TRUE`, all processes are passing in the same array of points. Otherwise, points need to be communicated among processes.
164 - ignoreOutsideDomain - If `PETSC_TRUE`, ignore points outside the domain, otherwise return an error
165 
166   Level: intermediate
167 
168 .seealso: [](ch_snes), `DM`, `DMInterpolationInfo`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`
169 @*/
170 PetscErrorCode DMInterpolationSetUp(DMInterpolationInfo ctx, DM dm, PetscBool redundantPoints, PetscBool ignoreOutsideDomain)
171 {
172   MPI_Comm           comm = ctx->comm;
173   PetscScalar       *a;
174   PetscInt           p, q, i;
175   PetscMPIInt        rank, size;
176   Vec                pointVec;
177   PetscSF            cellSF;
178   PetscLayout        layout;
179   PetscReal         *globalPoints;
180   PetscScalar       *globalPointsScalar;
181   const PetscInt    *ranges;
182   PetscMPIInt       *counts, *displs;
183   const PetscSFNode *foundCells;
184   const PetscInt    *foundPoints;
185   PetscMPIInt       *foundProcs, *globalProcs;
186   PetscInt           n, N, numFound;
187 
188   PetscFunctionBegin;
189   PetscValidHeaderSpecific(dm, DM_CLASSID, 2);
190   PetscCallMPI(MPI_Comm_size(comm, &size));
191   PetscCallMPI(MPI_Comm_rank(comm, &rank));
192   PetscCheck(ctx->dim >= 0, comm, PETSC_ERR_ARG_WRONGSTATE, "The spatial dimension has not been set");
193   /* Locate points */
194   n = ctx->nInput;
195   if (!redundantPoints) {
196     PetscCall(PetscLayoutCreate(comm, &layout));
197     PetscCall(PetscLayoutSetBlockSize(layout, 1));
198     PetscCall(PetscLayoutSetLocalSize(layout, n));
199     PetscCall(PetscLayoutSetUp(layout));
200     PetscCall(PetscLayoutGetSize(layout, &N));
201     /* Communicate all points to all processes */
202     PetscCall(PetscMalloc3(N * ctx->dim, &globalPoints, size, &counts, size, &displs));
203     PetscCall(PetscLayoutGetRanges(layout, &ranges));
204     for (p = 0; p < size; ++p) {
205       counts[p] = (ranges[p + 1] - ranges[p]) * ctx->dim;
206       displs[p] = ranges[p] * ctx->dim;
207     }
208     PetscCallMPI(MPI_Allgatherv(ctx->points, n * ctx->dim, MPIU_REAL, globalPoints, counts, displs, MPIU_REAL, comm));
209   } else {
210     N            = n;
211     globalPoints = ctx->points;
212     counts = displs = NULL;
213     layout          = NULL;
214   }
215 #if 0
216   PetscCall(PetscMalloc3(N,&foundCells,N,&foundProcs,N,&globalProcs));
217   /* foundCells[p] = m->locatePoint(&globalPoints[p*ctx->dim]); */
218 #else
219   #if defined(PETSC_USE_COMPLEX)
220   PetscCall(PetscMalloc1(N * ctx->dim, &globalPointsScalar));
221   for (i = 0; i < N * ctx->dim; i++) globalPointsScalar[i] = globalPoints[i];
222   #else
223   globalPointsScalar = globalPoints;
224   #endif
225   PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, ctx->dim, N * ctx->dim, globalPointsScalar, &pointVec));
226   PetscCall(PetscMalloc2(N, &foundProcs, N, &globalProcs));
227   for (p = 0; p < N; ++p) foundProcs[p] = size;
228   cellSF = NULL;
229   PetscCall(DMLocatePoints(dm, pointVec, DM_POINTLOCATION_REMOVE, &cellSF));
230   PetscCall(PetscSFGetGraph(cellSF, NULL, &numFound, &foundPoints, &foundCells));
231 #endif
232   for (p = 0; p < numFound; ++p) {
233     if (foundCells[p].index >= 0) foundProcs[foundPoints ? foundPoints[p] : p] = rank;
234   }
235   /* Let the lowest rank process own each point */
236   PetscCall(MPIU_Allreduce(foundProcs, globalProcs, N, MPI_INT, MPI_MIN, comm));
237   ctx->n = 0;
238   for (p = 0; p < N; ++p) {
239     if (globalProcs[p] == size) {
240       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),
241                  (double)(ctx->dim > 2 ? globalPoints[p * ctx->dim + 2] : 0.0));
242       if (rank == 0) ++ctx->n;
243     } else if (globalProcs[p] == rank) ++ctx->n;
244   }
245   /* Create coordinates vector and array of owned cells */
246   PetscCall(PetscMalloc1(ctx->n, &ctx->cells));
247   PetscCall(VecCreate(comm, &ctx->coords));
248   PetscCall(VecSetSizes(ctx->coords, ctx->n * ctx->dim, PETSC_DECIDE));
249   PetscCall(VecSetBlockSize(ctx->coords, ctx->dim));
250   PetscCall(VecSetType(ctx->coords, VECSTANDARD));
251   PetscCall(VecGetArray(ctx->coords, &a));
252   for (p = 0, q = 0, i = 0; p < N; ++p) {
253     if (globalProcs[p] == rank) {
254       PetscInt d;
255 
256       for (d = 0; d < ctx->dim; ++d, ++i) a[i] = globalPoints[p * ctx->dim + d];
257       ctx->cells[q] = foundCells[q].index;
258       ++q;
259     }
260     if (globalProcs[p] == size && rank == 0) {
261       PetscInt d;
262 
263       for (d = 0; d < ctx->dim; ++d, ++i) a[i] = 0.;
264       ctx->cells[q] = -1;
265       ++q;
266     }
267   }
268   PetscCall(VecRestoreArray(ctx->coords, &a));
269 #if 0
270   PetscCall(PetscFree3(foundCells,foundProcs,globalProcs));
271 #else
272   PetscCall(PetscFree2(foundProcs, globalProcs));
273   PetscCall(PetscSFDestroy(&cellSF));
274   PetscCall(VecDestroy(&pointVec));
275 #endif
276   if ((void *)globalPointsScalar != (void *)globalPoints) PetscCall(PetscFree(globalPointsScalar));
277   if (!redundantPoints) PetscCall(PetscFree3(globalPoints, counts, displs));
278   PetscCall(PetscLayoutDestroy(&layout));
279   PetscFunctionReturn(PETSC_SUCCESS);
280 }
281 
282 /*@C
283   DMInterpolationGetCoordinates - Gets a `Vec` with the coordinates of each interpolation point
284 
285   Collective
286 
287   Input Parameter:
288 . ctx - the context
289 
290   Output Parameter:
291 . coordinates - the coordinates of interpolation points
292 
293   Level: intermediate
294 
295   Note:
296   The local vector entries correspond to interpolation points lying on this process, according to the associated `DM`.
297   This is a borrowed vector that the user should not destroy.
298 
299 .seealso: [](ch_snes), `DM`, `DMInterpolationInfo`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`
300 @*/
301 PetscErrorCode DMInterpolationGetCoordinates(DMInterpolationInfo ctx, Vec *coordinates)
302 {
303   PetscFunctionBegin;
304   PetscAssertPointer(coordinates, 2);
305   PetscCheck(ctx->coords, ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup.");
306   *coordinates = ctx->coords;
307   PetscFunctionReturn(PETSC_SUCCESS);
308 }
309 
310 /*@C
311   DMInterpolationGetVector - Gets a `Vec` which can hold all the interpolated field values
312 
313   Collective
314 
315   Input Parameter:
316 . ctx - the context
317 
318   Output Parameter:
319 . v - a vector capable of holding the interpolated field values
320 
321   Level: intermediate
322 
323   Note:
324   This vector should be returned using `DMInterpolationRestoreVector()`.
325 
326 .seealso: [](ch_snes), `DM`, `DMInterpolationInfo`, `DMInterpolationRestoreVector()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`
327 @*/
328 PetscErrorCode DMInterpolationGetVector(DMInterpolationInfo ctx, Vec *v)
329 {
330   PetscFunctionBegin;
331   PetscAssertPointer(v, 2);
332   PetscCheck(ctx->coords, ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup.");
333   PetscCall(VecCreate(ctx->comm, v));
334   PetscCall(VecSetSizes(*v, ctx->n * ctx->dof, PETSC_DECIDE));
335   PetscCall(VecSetBlockSize(*v, ctx->dof));
336   PetscCall(VecSetType(*v, VECSTANDARD));
337   PetscFunctionReturn(PETSC_SUCCESS);
338 }
339 
340 /*@C
341   DMInterpolationRestoreVector - Returns a `Vec` which can hold all the interpolated field values
342 
343   Collective
344 
345   Input Parameters:
346 + ctx - the context
347 - v   - a vector capable of holding the interpolated field values
348 
349   Level: intermediate
350 
351 .seealso: [](ch_snes), `DM`, `DMInterpolationInfo`, `DMInterpolationGetVector()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`
352 @*/
353 PetscErrorCode DMInterpolationRestoreVector(DMInterpolationInfo ctx, Vec *v)
354 {
355   PetscFunctionBegin;
356   PetscAssertPointer(v, 2);
357   PetscCheck(ctx->coords, ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup.");
358   PetscCall(VecDestroy(v));
359   PetscFunctionReturn(PETSC_SUCCESS);
360 }
361 
362 static inline PetscErrorCode DMInterpolate_Segment_Private(DMInterpolationInfo ctx, DM dm, PetscInt p, Vec xLocal, Vec v)
363 {
364   const PetscInt     c   = ctx->cells[p];
365   const PetscInt     dof = ctx->dof;
366   PetscScalar       *x   = NULL;
367   const PetscScalar *coords;
368   PetscScalar       *a;
369   PetscReal          v0, J, invJ, detJ, xir[1];
370   PetscInt           xSize, comp;
371 
372   PetscFunctionBegin;
373   PetscCall(VecGetArrayRead(ctx->coords, &coords));
374   PetscCall(VecGetArray(v, &a));
375   PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, &v0, &J, &invJ, &detJ));
376   PetscCheck(detJ > 0.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %" PetscInt_FMT, (double)detJ, c);
377   xir[0] = invJ * PetscRealPart(coords[p] - v0);
378   PetscCall(DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x));
379   if (2 * dof == xSize) {
380     for (comp = 0; comp < dof; ++comp) a[p * dof + comp] = x[0 * dof + comp] * (1 - xir[0]) + x[1 * dof + comp] * xir[0];
381   } else if (dof == xSize) {
382     for (comp = 0; comp < dof; ++comp) a[p * dof + comp] = x[0 * dof + comp];
383   } 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);
384   PetscCall(DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x));
385   PetscCall(VecRestoreArray(v, &a));
386   PetscCall(VecRestoreArrayRead(ctx->coords, &coords));
387   PetscFunctionReturn(PETSC_SUCCESS);
388 }
389 
390 static inline PetscErrorCode DMInterpolate_Triangle_Private(DMInterpolationInfo ctx, DM dm, PetscInt p, Vec xLocal, Vec v)
391 {
392   const PetscInt     c = ctx->cells[p];
393   PetscScalar       *x = NULL;
394   const PetscScalar *coords;
395   PetscScalar       *a;
396   PetscReal         *v0, *J, *invJ, detJ;
397   PetscReal          xi[4];
398 
399   PetscFunctionBegin;
400   PetscCall(PetscMalloc3(ctx->dim, &v0, ctx->dim * ctx->dim, &J, ctx->dim * ctx->dim, &invJ));
401   PetscCall(VecGetArrayRead(ctx->coords, &coords));
402   PetscCall(VecGetArray(v, &a));
403   PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ));
404   PetscCheck(detJ > 0.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %" PetscInt_FMT, (double)detJ, c);
405   PetscCall(DMPlexVecGetClosure(dm, NULL, xLocal, c, NULL, &x));
406   for (PetscInt comp = 0; comp < ctx->dof; ++comp) a[p * ctx->dof + comp] = x[0 * ctx->dof + comp];
407   for (PetscInt d = 0; d < ctx->dim; ++d) {
408     xi[d] = 0.0;
409     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]);
410     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];
411   }
412   PetscCall(DMPlexVecRestoreClosure(dm, NULL, xLocal, c, NULL, &x));
413   PetscCall(VecRestoreArray(v, &a));
414   PetscCall(VecRestoreArrayRead(ctx->coords, &coords));
415   PetscCall(PetscFree3(v0, J, invJ));
416   PetscFunctionReturn(PETSC_SUCCESS);
417 }
418 
419 static inline PetscErrorCode DMInterpolate_Tetrahedron_Private(DMInterpolationInfo ctx, DM dm, PetscInt p, Vec xLocal, Vec v)
420 {
421   const PetscInt     c        = ctx->cells[p];
422   const PetscInt     order[3] = {2, 1, 3};
423   PetscScalar       *x        = NULL;
424   PetscReal         *v0, *J, *invJ, detJ;
425   const PetscScalar *coords;
426   PetscScalar       *a;
427   PetscReal          xi[4];
428 
429   PetscFunctionBegin;
430   PetscCall(PetscMalloc3(ctx->dim, &v0, ctx->dim * ctx->dim, &J, ctx->dim * ctx->dim, &invJ));
431   PetscCall(VecGetArrayRead(ctx->coords, &coords));
432   PetscCall(VecGetArray(v, &a));
433   PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ));
434   PetscCheck(detJ > 0.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %" PetscInt_FMT, (double)detJ, c);
435   PetscCall(DMPlexVecGetClosure(dm, NULL, xLocal, c, NULL, &x));
436   for (PetscInt comp = 0; comp < ctx->dof; ++comp) a[p * ctx->dof + comp] = x[0 * ctx->dof + comp];
437   for (PetscInt d = 0; d < ctx->dim; ++d) {
438     xi[d] = 0.0;
439     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]);
440     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];
441   }
442   PetscCall(DMPlexVecRestoreClosure(dm, NULL, xLocal, c, NULL, &x));
443   PetscCall(VecRestoreArray(v, &a));
444   PetscCall(VecRestoreArrayRead(ctx->coords, &coords));
445   PetscCall(PetscFree3(v0, J, invJ));
446   PetscFunctionReturn(PETSC_SUCCESS);
447 }
448 
449 static inline PetscErrorCode QuadMap_Private(SNES snes, Vec Xref, Vec Xreal, void *ctx)
450 {
451   const PetscScalar *vertices = (const PetscScalar *)ctx;
452   const PetscScalar  x0       = vertices[0];
453   const PetscScalar  y0       = vertices[1];
454   const PetscScalar  x1       = vertices[2];
455   const PetscScalar  y1       = vertices[3];
456   const PetscScalar  x2       = vertices[4];
457   const PetscScalar  y2       = vertices[5];
458   const PetscScalar  x3       = vertices[6];
459   const PetscScalar  y3       = vertices[7];
460   const PetscScalar  f_1      = x1 - x0;
461   const PetscScalar  g_1      = y1 - y0;
462   const PetscScalar  f_3      = x3 - x0;
463   const PetscScalar  g_3      = y3 - y0;
464   const PetscScalar  f_01     = x2 - x1 - x3 + x0;
465   const PetscScalar  g_01     = y2 - y1 - y3 + y0;
466   const PetscScalar *ref;
467   PetscScalar       *real;
468 
469   PetscFunctionBegin;
470   PetscCall(VecGetArrayRead(Xref, &ref));
471   PetscCall(VecGetArray(Xreal, &real));
472   {
473     const PetscScalar p0 = ref[0];
474     const PetscScalar p1 = ref[1];
475 
476     real[0] = x0 + f_1 * p0 + f_3 * p1 + f_01 * p0 * p1;
477     real[1] = y0 + g_1 * p0 + g_3 * p1 + g_01 * p0 * p1;
478   }
479   PetscCall(PetscLogFlops(28));
480   PetscCall(VecRestoreArrayRead(Xref, &ref));
481   PetscCall(VecRestoreArray(Xreal, &real));
482   PetscFunctionReturn(PETSC_SUCCESS);
483 }
484 
485 #include <petsc/private/dmimpl.h>
486 static inline PetscErrorCode QuadJacobian_Private(SNES snes, Vec Xref, Mat J, Mat M, void *ctx)
487 {
488   const PetscScalar *vertices = (const PetscScalar *)ctx;
489   const PetscScalar  x0       = vertices[0];
490   const PetscScalar  y0       = vertices[1];
491   const PetscScalar  x1       = vertices[2];
492   const PetscScalar  y1       = vertices[3];
493   const PetscScalar  x2       = vertices[4];
494   const PetscScalar  y2       = vertices[5];
495   const PetscScalar  x3       = vertices[6];
496   const PetscScalar  y3       = vertices[7];
497   const PetscScalar  f_01     = x2 - x1 - x3 + x0;
498   const PetscScalar  g_01     = y2 - y1 - y3 + y0;
499   const PetscScalar *ref;
500 
501   PetscFunctionBegin;
502   PetscCall(VecGetArrayRead(Xref, &ref));
503   {
504     const PetscScalar x       = ref[0];
505     const PetscScalar y       = ref[1];
506     const PetscInt    rows[2] = {0, 1};
507     PetscScalar       values[4];
508 
509     values[0] = (x1 - x0 + f_01 * y) * 0.5;
510     values[1] = (x3 - x0 + f_01 * x) * 0.5;
511     values[2] = (y1 - y0 + g_01 * y) * 0.5;
512     values[3] = (y3 - y0 + g_01 * x) * 0.5;
513     PetscCall(MatSetValues(J, 2, rows, 2, rows, values, INSERT_VALUES));
514   }
515   PetscCall(PetscLogFlops(30));
516   PetscCall(VecRestoreArrayRead(Xref, &ref));
517   PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
518   PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
519   PetscFunctionReturn(PETSC_SUCCESS);
520 }
521 
522 static inline PetscErrorCode DMInterpolate_Quad_Private(DMInterpolationInfo ctx, DM dm, PetscInt p, Vec xLocal, Vec v)
523 {
524   const PetscInt     c   = ctx->cells[p];
525   PetscFE            fem = NULL;
526   DM                 dmCoord;
527   SNES               snes;
528   KSP                ksp;
529   PC                 pc;
530   Vec                coordsLocal, r, ref, real;
531   Mat                J;
532   PetscTabulation    T = NULL;
533   const PetscScalar *coords;
534   PetscScalar       *a;
535   PetscReal          xir[2] = {0., 0.};
536   PetscInt           Nf;
537   const PetscInt     dof = ctx->dof;
538   PetscScalar       *x = NULL, *vertices = NULL;
539   PetscScalar       *xi;
540   PetscInt           coordSize, xSize;
541 
542   PetscFunctionBegin;
543   PetscCall(DMGetNumFields(dm, &Nf));
544   if (Nf) {
545     PetscObject  obj;
546     PetscClassId id;
547 
548     PetscCall(DMGetField(dm, 0, NULL, &obj));
549     PetscCall(PetscObjectGetClassId(obj, &id));
550     if (id == PETSCFE_CLASSID) {
551       fem = (PetscFE)obj;
552       PetscCall(PetscFECreateTabulation(fem, 1, 1, xir, 0, &T));
553     }
554   }
555   PetscCall(DMGetCoordinatesLocal(dm, &coordsLocal));
556   PetscCall(DMGetCoordinateDM(dm, &dmCoord));
557   PetscCall(SNESCreate(PETSC_COMM_SELF, &snes));
558   PetscCall(SNESSetOptionsPrefix(snes, "quad_interp_"));
559   PetscCall(VecCreate(PETSC_COMM_SELF, &r));
560   PetscCall(VecSetSizes(r, 2, 2));
561   PetscCall(VecSetType(r, dm->vectype));
562   PetscCall(VecDuplicate(r, &ref));
563   PetscCall(VecDuplicate(r, &real));
564   PetscCall(MatCreate(PETSC_COMM_SELF, &J));
565   PetscCall(MatSetSizes(J, 2, 2, 2, 2));
566   PetscCall(MatSetType(J, MATSEQDENSE));
567   PetscCall(MatSetUp(J));
568   PetscCall(SNESSetFunction(snes, r, QuadMap_Private, NULL));
569   PetscCall(SNESSetJacobian(snes, J, J, QuadJacobian_Private, NULL));
570   PetscCall(SNESGetKSP(snes, &ksp));
571   PetscCall(KSPGetPC(ksp, &pc));
572   PetscCall(PCSetType(pc, PCLU));
573   PetscCall(SNESSetFromOptions(snes));
574 
575   PetscCall(VecGetArrayRead(ctx->coords, &coords));
576   PetscCall(VecGetArray(v, &a));
577   PetscCall(DMPlexVecGetClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices));
578   PetscCheck(4 * 2 == coordSize, ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid closure size %" PetscInt_FMT " should be %d", coordSize, 4 * 2);
579   PetscCall(DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x));
580   PetscCall(SNESSetFunction(snes, NULL, NULL, vertices));
581   PetscCall(SNESSetJacobian(snes, NULL, NULL, NULL, vertices));
582   PetscCall(VecGetArray(real, &xi));
583   xi[0] = coords[p * ctx->dim + 0];
584   xi[1] = coords[p * ctx->dim + 1];
585   PetscCall(VecRestoreArray(real, &xi));
586   PetscCall(SNESSolve(snes, real, ref));
587   PetscCall(VecGetArray(ref, &xi));
588   xir[0] = PetscRealPart(xi[0]);
589   xir[1] = PetscRealPart(xi[1]);
590   if (4 * dof == xSize) {
591     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];
592   } else if (dof == xSize) {
593     for (PetscInt comp = 0; comp < dof; ++comp) a[p * dof + comp] = x[0 * dof + comp];
594   } else {
595     PetscCheck(fem, ctx->comm, PETSC_ERR_ARG_WRONG, "Cannot have a higher order interpolant if the discretization is not PetscFE");
596     xir[0] = 2.0 * xir[0] - 1.0;
597     xir[1] = 2.0 * xir[1] - 1.0;
598     PetscCall(PetscFEComputeTabulation(fem, 1, xir, 0, T));
599     for (PetscInt comp = 0; comp < dof; ++comp) {
600       a[p * dof + comp] = 0.0;
601       for (PetscInt d = 0; d < xSize / dof; ++d) a[p * dof + comp] += x[d * dof + comp] * T->T[0][d * dof + comp];
602     }
603   }
604   PetscCall(VecRestoreArray(ref, &xi));
605   PetscCall(DMPlexVecRestoreClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices));
606   PetscCall(DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x));
607   PetscCall(PetscTabulationDestroy(&T));
608   PetscCall(VecRestoreArray(v, &a));
609   PetscCall(VecRestoreArrayRead(ctx->coords, &coords));
610 
611   PetscCall(SNESDestroy(&snes));
612   PetscCall(VecDestroy(&r));
613   PetscCall(VecDestroy(&ref));
614   PetscCall(VecDestroy(&real));
615   PetscCall(MatDestroy(&J));
616   PetscFunctionReturn(PETSC_SUCCESS);
617 }
618 
619 static inline PetscErrorCode HexMap_Private(SNES snes, Vec Xref, Vec Xreal, void *ctx)
620 {
621   const PetscScalar *vertices = (const PetscScalar *)ctx;
622   const PetscScalar  x0       = vertices[0];
623   const PetscScalar  y0       = vertices[1];
624   const PetscScalar  z0       = vertices[2];
625   const PetscScalar  x1       = vertices[9];
626   const PetscScalar  y1       = vertices[10];
627   const PetscScalar  z1       = vertices[11];
628   const PetscScalar  x2       = vertices[6];
629   const PetscScalar  y2       = vertices[7];
630   const PetscScalar  z2       = vertices[8];
631   const PetscScalar  x3       = vertices[3];
632   const PetscScalar  y3       = vertices[4];
633   const PetscScalar  z3       = vertices[5];
634   const PetscScalar  x4       = vertices[12];
635   const PetscScalar  y4       = vertices[13];
636   const PetscScalar  z4       = vertices[14];
637   const PetscScalar  x5       = vertices[15];
638   const PetscScalar  y5       = vertices[16];
639   const PetscScalar  z5       = vertices[17];
640   const PetscScalar  x6       = vertices[18];
641   const PetscScalar  y6       = vertices[19];
642   const PetscScalar  z6       = vertices[20];
643   const PetscScalar  x7       = vertices[21];
644   const PetscScalar  y7       = vertices[22];
645   const PetscScalar  z7       = vertices[23];
646   const PetscScalar  f_1      = x1 - x0;
647   const PetscScalar  g_1      = y1 - y0;
648   const PetscScalar  h_1      = z1 - z0;
649   const PetscScalar  f_3      = x3 - x0;
650   const PetscScalar  g_3      = y3 - y0;
651   const PetscScalar  h_3      = z3 - z0;
652   const PetscScalar  f_4      = x4 - x0;
653   const PetscScalar  g_4      = y4 - y0;
654   const PetscScalar  h_4      = z4 - z0;
655   const PetscScalar  f_01     = x2 - x1 - x3 + x0;
656   const PetscScalar  g_01     = y2 - y1 - y3 + y0;
657   const PetscScalar  h_01     = z2 - z1 - z3 + z0;
658   const PetscScalar  f_12     = x7 - x3 - x4 + x0;
659   const PetscScalar  g_12     = y7 - y3 - y4 + y0;
660   const PetscScalar  h_12     = z7 - z3 - z4 + z0;
661   const PetscScalar  f_02     = x5 - x1 - x4 + x0;
662   const PetscScalar  g_02     = y5 - y1 - y4 + y0;
663   const PetscScalar  h_02     = z5 - z1 - z4 + z0;
664   const PetscScalar  f_012    = x6 - x0 + x1 - x2 + x3 + x4 - x5 - x7;
665   const PetscScalar  g_012    = y6 - y0 + y1 - y2 + y3 + y4 - y5 - y7;
666   const PetscScalar  h_012    = z6 - z0 + z1 - z2 + z3 + z4 - z5 - z7;
667   const PetscScalar *ref;
668   PetscScalar       *real;
669 
670   PetscFunctionBegin;
671   PetscCall(VecGetArrayRead(Xref, &ref));
672   PetscCall(VecGetArray(Xreal, &real));
673   {
674     const PetscScalar p0 = ref[0];
675     const PetscScalar p1 = ref[1];
676     const PetscScalar p2 = ref[2];
677 
678     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;
679     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;
680     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;
681   }
682   PetscCall(PetscLogFlops(114));
683   PetscCall(VecRestoreArrayRead(Xref, &ref));
684   PetscCall(VecRestoreArray(Xreal, &real));
685   PetscFunctionReturn(PETSC_SUCCESS);
686 }
687 
688 static inline PetscErrorCode HexJacobian_Private(SNES snes, Vec Xref, Mat J, Mat M, void *ctx)
689 {
690   const PetscScalar *vertices = (const PetscScalar *)ctx;
691   const PetscScalar  x0       = vertices[0];
692   const PetscScalar  y0       = vertices[1];
693   const PetscScalar  z0       = vertices[2];
694   const PetscScalar  x1       = vertices[9];
695   const PetscScalar  y1       = vertices[10];
696   const PetscScalar  z1       = vertices[11];
697   const PetscScalar  x2       = vertices[6];
698   const PetscScalar  y2       = vertices[7];
699   const PetscScalar  z2       = vertices[8];
700   const PetscScalar  x3       = vertices[3];
701   const PetscScalar  y3       = vertices[4];
702   const PetscScalar  z3       = vertices[5];
703   const PetscScalar  x4       = vertices[12];
704   const PetscScalar  y4       = vertices[13];
705   const PetscScalar  z4       = vertices[14];
706   const PetscScalar  x5       = vertices[15];
707   const PetscScalar  y5       = vertices[16];
708   const PetscScalar  z5       = vertices[17];
709   const PetscScalar  x6       = vertices[18];
710   const PetscScalar  y6       = vertices[19];
711   const PetscScalar  z6       = vertices[20];
712   const PetscScalar  x7       = vertices[21];
713   const PetscScalar  y7       = vertices[22];
714   const PetscScalar  z7       = vertices[23];
715   const PetscScalar  f_xy     = x2 - x1 - x3 + x0;
716   const PetscScalar  g_xy     = y2 - y1 - y3 + y0;
717   const PetscScalar  h_xy     = z2 - z1 - z3 + z0;
718   const PetscScalar  f_yz     = x7 - x3 - x4 + x0;
719   const PetscScalar  g_yz     = y7 - y3 - y4 + y0;
720   const PetscScalar  h_yz     = z7 - z3 - z4 + z0;
721   const PetscScalar  f_xz     = x5 - x1 - x4 + x0;
722   const PetscScalar  g_xz     = y5 - y1 - y4 + y0;
723   const PetscScalar  h_xz     = z5 - z1 - z4 + z0;
724   const PetscScalar  f_xyz    = x6 - x0 + x1 - x2 + x3 + x4 - x5 - x7;
725   const PetscScalar  g_xyz    = y6 - y0 + y1 - y2 + y3 + y4 - y5 - y7;
726   const PetscScalar  h_xyz    = z6 - z0 + z1 - z2 + z3 + z4 - z5 - z7;
727   const PetscScalar *ref;
728 
729   PetscFunctionBegin;
730   PetscCall(VecGetArrayRead(Xref, &ref));
731   {
732     const PetscScalar x       = ref[0];
733     const PetscScalar y       = ref[1];
734     const PetscScalar z       = ref[2];
735     const PetscInt    rows[3] = {0, 1, 2};
736     PetscScalar       values[9];
737 
738     values[0] = (x1 - x0 + f_xy * y + f_xz * z + f_xyz * y * z) / 2.0;
739     values[1] = (x3 - x0 + f_xy * x + f_yz * z + f_xyz * x * z) / 2.0;
740     values[2] = (x4 - x0 + f_yz * y + f_xz * x + f_xyz * x * y) / 2.0;
741     values[3] = (y1 - y0 + g_xy * y + g_xz * z + g_xyz * y * z) / 2.0;
742     values[4] = (y3 - y0 + g_xy * x + g_yz * z + g_xyz * x * z) / 2.0;
743     values[5] = (y4 - y0 + g_yz * y + g_xz * x + g_xyz * x * y) / 2.0;
744     values[6] = (z1 - z0 + h_xy * y + h_xz * z + h_xyz * y * z) / 2.0;
745     values[7] = (z3 - z0 + h_xy * x + h_yz * z + h_xyz * x * z) / 2.0;
746     values[8] = (z4 - z0 + h_yz * y + h_xz * x + h_xyz * x * y) / 2.0;
747 
748     PetscCall(MatSetValues(J, 3, rows, 3, rows, values, INSERT_VALUES));
749   }
750   PetscCall(PetscLogFlops(152));
751   PetscCall(VecRestoreArrayRead(Xref, &ref));
752   PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
753   PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
754   PetscFunctionReturn(PETSC_SUCCESS);
755 }
756 
757 static inline PetscErrorCode DMInterpolate_Hex_Private(DMInterpolationInfo ctx, DM dm, PetscInt p, Vec xLocal, Vec v)
758 {
759   const PetscInt     c = ctx->cells[p];
760   DM                 dmCoord;
761   SNES               snes;
762   KSP                ksp;
763   PC                 pc;
764   Vec                coordsLocal, r, ref, real;
765   Mat                J;
766   const PetscScalar *coords;
767   PetscScalar       *a;
768   PetscScalar       *x = NULL, *vertices = NULL;
769   PetscScalar       *xi;
770   PetscReal          xir[3];
771   PetscInt           coordSize, xSize;
772 
773   PetscFunctionBegin;
774   PetscCall(DMGetCoordinatesLocal(dm, &coordsLocal));
775   PetscCall(DMGetCoordinateDM(dm, &dmCoord));
776   PetscCall(SNESCreate(PETSC_COMM_SELF, &snes));
777   PetscCall(SNESSetOptionsPrefix(snes, "hex_interp_"));
778   PetscCall(VecCreate(PETSC_COMM_SELF, &r));
779   PetscCall(VecSetSizes(r, 3, 3));
780   PetscCall(VecSetType(r, dm->vectype));
781   PetscCall(VecDuplicate(r, &ref));
782   PetscCall(VecDuplicate(r, &real));
783   PetscCall(MatCreate(PETSC_COMM_SELF, &J));
784   PetscCall(MatSetSizes(J, 3, 3, 3, 3));
785   PetscCall(MatSetType(J, MATSEQDENSE));
786   PetscCall(MatSetUp(J));
787   PetscCall(SNESSetFunction(snes, r, HexMap_Private, NULL));
788   PetscCall(SNESSetJacobian(snes, J, J, HexJacobian_Private, NULL));
789   PetscCall(SNESGetKSP(snes, &ksp));
790   PetscCall(KSPGetPC(ksp, &pc));
791   PetscCall(PCSetType(pc, PCLU));
792   PetscCall(SNESSetFromOptions(snes));
793 
794   PetscCall(VecGetArrayRead(ctx->coords, &coords));
795   PetscCall(VecGetArray(v, &a));
796   PetscCall(DMPlexVecGetClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices));
797   PetscCheck(8 * 3 == coordSize, ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid coordinate closure size %" PetscInt_FMT " should be %d", coordSize, 8 * 3);
798   PetscCall(DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x));
799   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);
800   PetscCall(SNESSetFunction(snes, NULL, NULL, vertices));
801   PetscCall(SNESSetJacobian(snes, NULL, NULL, NULL, vertices));
802   PetscCall(VecGetArray(real, &xi));
803   xi[0] = coords[p * ctx->dim + 0];
804   xi[1] = coords[p * ctx->dim + 1];
805   xi[2] = coords[p * ctx->dim + 2];
806   PetscCall(VecRestoreArray(real, &xi));
807   PetscCall(SNESSolve(snes, real, ref));
808   PetscCall(VecGetArray(ref, &xi));
809   xir[0] = PetscRealPart(xi[0]);
810   xir[1] = PetscRealPart(xi[1]);
811   xir[2] = PetscRealPart(xi[2]);
812   if (8 * ctx->dof == xSize) {
813     for (PetscInt comp = 0; comp < ctx->dof; ++comp) {
814       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]) +
815                                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];
816     }
817   } else {
818     for (PetscInt comp = 0; comp < ctx->dof; ++comp) a[p * ctx->dof + comp] = x[0 * ctx->dof + comp];
819   }
820   PetscCall(VecRestoreArray(ref, &xi));
821   PetscCall(DMPlexVecRestoreClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices));
822   PetscCall(DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x));
823   PetscCall(VecRestoreArray(v, &a));
824   PetscCall(VecRestoreArrayRead(ctx->coords, &coords));
825 
826   PetscCall(SNESDestroy(&snes));
827   PetscCall(VecDestroy(&r));
828   PetscCall(VecDestroy(&ref));
829   PetscCall(VecDestroy(&real));
830   PetscCall(MatDestroy(&J));
831   PetscFunctionReturn(PETSC_SUCCESS);
832 }
833 
834 /*@C
835   DMInterpolationEvaluate - Using the input from dm and x, calculates interpolated field values at the interpolation points.
836 
837   Input Parameters:
838 + ctx - The `DMInterpolationInfo` context
839 . dm  - The `DM`
840 - x   - The local vector containing the field to be interpolated
841 
842   Output Parameter:
843 . v - The vector containing the interpolated values
844 
845   Level: beginner
846 
847   Note:
848   A suitable `v` can be obtained using `DMInterpolationGetVector()`.
849 
850 .seealso: [](ch_snes), `DM`, `DMInterpolationInfo`, `DMInterpolationGetVector()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`
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;
857 
858   PetscFunctionBegin;
859   PetscValidHeaderSpecific(dm, DM_CLASSID, 2);
860   PetscValidHeaderSpecific(x, VEC_CLASSID, 3);
861   PetscValidHeaderSpecific(v, VEC_CLASSID, 4);
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;
872 
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;
885 
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;
893 
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;
902 
903         PetscCall(PetscDSGetDiscretization(ds, field, &obj));
904         PetscCall(PetscObjectGetClassId(obj, &id));
905         if (id == PETSCFE_CLASSID) {
906           PetscFE fe = (PetscFE)obj;
907 
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;
913 
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;
925 
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;
943 
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 }
968 
969 /*@C
970   DMInterpolationDestroy - Destroys a `DMInterpolationInfo` context
971 
972   Collective
973 
974   Input Parameter:
975 . ctx - the context
976 
977   Level: beginner
978 
979 .seealso: [](ch_snes), `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 }
992