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, obtained with `DMInterpolationGetCoordinates()`
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: }