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: /*@
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: /*@
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: /*@
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: /*@
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: /*@
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: /*@
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: /*@
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, in;
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: PetscCall(PetscMPIIntCast((ranges[p + 1] - ranges[p]) * ctx->dim, &counts[p]));
209: PetscCall(PetscMPIIntCast(ranges[p] * ctx->dim, &displs[p]));
210: }
211: PetscCall(PetscMPIIntCast(n * ctx->dim, &in));
212: PetscCallMPI(MPI_Allgatherv(ctx->points, in, MPIU_REAL, globalPoints, counts, displs, MPIU_REAL, comm));
213: } else {
214: N = n;
215: globalPoints = ctx->points;
216: counts = displs = NULL;
217: layout = NULL;
218: }
219: #if 0
220: PetscCall(PetscMalloc3(N,&foundCells,N,&foundProcs,N,&globalProcs));
221: /* foundCells[p] = m->locatePoint(&globalPoints[p*ctx->dim]); */
222: #else
223: #if defined(PETSC_USE_COMPLEX)
224: PetscCall(PetscMalloc1(N * ctx->dim, &globalPointsScalar));
225: for (i = 0; i < N * ctx->dim; i++) globalPointsScalar[i] = globalPoints[i];
226: #else
227: globalPointsScalar = globalPoints;
228: #endif
229: PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, ctx->dim, N * ctx->dim, globalPointsScalar, &pointVec));
230: PetscCall(PetscMalloc2(N, &foundProcs, N, &globalProcs));
231: for (p = 0; p < N; ++p) foundProcs[p] = size;
232: cellSF = NULL;
233: PetscCall(DMLocatePoints(dm, pointVec, DM_POINTLOCATION_REMOVE, &cellSF));
234: PetscCall(PetscSFGetGraph(cellSF, NULL, &numFound, &foundPoints, &foundCells));
235: #endif
236: for (p = 0; p < numFound; ++p) {
237: if (foundCells[p].index >= 0) foundProcs[foundPoints ? foundPoints[p] : p] = rank;
238: }
239: /* Let the lowest rank process own each point */
240: PetscCallMPI(MPIU_Allreduce(foundProcs, globalProcs, N, MPI_INT, MPI_MIN, comm));
241: ctx->n = 0;
242: for (p = 0; p < N; ++p) {
243: if (globalProcs[p] == size) {
244: 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),
245: (double)(ctx->dim > 2 ? globalPoints[p * ctx->dim + 2] : 0));
246: if (rank == 0) ++ctx->n;
247: } else if (globalProcs[p] == rank) ++ctx->n;
248: }
249: /* Create coordinates vector and array of owned cells */
250: PetscCall(PetscMalloc1(ctx->n, &ctx->cells));
251: PetscCall(VecCreate(comm, &ctx->coords));
252: PetscCall(VecSetSizes(ctx->coords, ctx->n * ctx->dim, PETSC_DECIDE));
253: PetscCall(VecSetBlockSize(ctx->coords, ctx->dim));
254: PetscCall(VecSetType(ctx->coords, VECSTANDARD));
255: PetscCall(VecGetArray(ctx->coords, &a));
256: for (p = 0, q = 0, i = 0; p < N; ++p) {
257: if (globalProcs[p] == rank) {
258: for (PetscInt d = 0; d < ctx->dim; ++d, ++i) a[i] = globalPoints[p * ctx->dim + d];
259: ctx->cells[q] = foundCells[q].index;
260: ++q;
261: }
262: if (globalProcs[p] == size && rank == 0) {
263: for (PetscInt 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: }
282: /*@
283: DMInterpolationGetCoordinates - Gets a `Vec` with the coordinates of each interpolation point
285: Collective
287: Input Parameter:
288: . ctx - the context
290: Output Parameter:
291: . coordinates - the coordinates of interpolation points
293: Level: intermediate
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.
299: .seealso: [](ch_dmbase), `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: }
310: /*@
311: DMInterpolationGetVector - Gets a `Vec` which can hold all the interpolated field values
313: Collective
315: Input Parameter:
316: . ctx - the context
318: Output Parameter:
319: . v - a vector capable of holding the interpolated field values
321: Level: intermediate
323: Note:
324: This vector should be returned using `DMInterpolationRestoreVector()`.
326: .seealso: [](ch_dmbase), `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: }
340: /*@
341: DMInterpolationRestoreVector - Returns a `Vec` which can hold all the interpolated field values
343: Collective
345: Input Parameters:
346: + ctx - the context
347: - v - a vector capable of holding the interpolated field values
349: Level: intermediate
351: .seealso: [](ch_dmbase), `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: }
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;
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 (PetscInt 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 (PetscInt 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: }
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];
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: }
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];
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: }
449: static inline PetscErrorCode QuadMap_Private(SNES snes, Vec Xref, Vec Xreal, PetscCtx 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;
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];
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: }
485: #include <petsc/private/dmimpl.h>
486: static inline PetscErrorCode QuadJacobian_Private(SNES snes, Vec Xref, Mat J, Mat M, PetscCtx 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;
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];
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: }
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;
542: PetscFunctionBegin;
543: PetscCall(DMGetNumFields(dm, &Nf));
544: if (Nf) {
545: PetscObject obj;
546: PetscClassId id;
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));
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));
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: }
619: static inline PetscErrorCode HexMap_Private(SNES snes, Vec Xref, Vec Xreal, PetscCtx 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;
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];
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: }
688: static inline PetscErrorCode HexJacobian_Private(SNES snes, Vec Xref, Mat J, Mat M, PetscCtx 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;
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];
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;
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: }
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;
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));
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));
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: }
834: /*@
835: DMInterpolationEvaluate - Using the input from `dm` and `x`, calculates interpolated field values at the interpolation points.
837: Input Parameters:
838: + ctx - The `DMInterpolationInfo` context obtained with `DMInterpolationCreate()`
839: . dm - The `DM`
840: - x - The local vector containing the field to be interpolated, can be created with `DMCreateGlobalVector()`
842: Output Parameter:
843: . v - The vector containing the interpolated values, obtained with `DMInterpolationGetVector()`
845: Level: beginner
847: .seealso: [](ch_dmbase), `DM`, `DMInterpolationInfo`, `DMInterpolationGetVector()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`, `DMInterpolationGetCoordinates()`
848: @*/
849: PetscErrorCode DMInterpolationEvaluate(DMInterpolationInfo ctx, DM dm, Vec x, Vec v)
850: {
851: PetscDS ds;
852: PetscInt n, p, Nf, field;
853: PetscBool useDS = PETSC_FALSE;
855: PetscFunctionBegin;
859: PetscCall(VecGetLocalSize(v, &n));
860: 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);
861: if (!n) PetscFunctionReturn(PETSC_SUCCESS);
862: PetscCall(DMGetDS(dm, &ds));
863: if (ds) {
864: useDS = PETSC_TRUE;
865: PetscCall(PetscDSGetNumFields(ds, &Nf));
866: for (field = 0; field < Nf; ++field) {
867: PetscObject obj;
868: PetscClassId id;
870: PetscCall(PetscDSGetDiscretization(ds, field, &obj));
871: PetscCall(PetscObjectGetClassId(obj, &id));
872: if (id != PETSCFE_CLASSID && id != PETSCFV_CLASSID) {
873: useDS = PETSC_FALSE;
874: break;
875: }
876: }
877: }
878: if (useDS) {
879: const PetscScalar *coords;
880: PetscScalar *interpolant;
881: PetscInt cdim;
883: PetscCall(DMGetCoordinateDim(dm, &cdim));
884: PetscCall(VecGetArrayRead(ctx->coords, &coords));
885: PetscCall(VecGetArrayWrite(v, &interpolant));
886: for (p = 0; p < ctx->n; ++p) {
887: PetscReal pcoords[3], xi[3];
888: PetscScalar *xa = NULL;
889: PetscInt coff = 0, foff = 0, clSize;
891: if (ctx->cells[p] < 0) continue;
892: for (PetscInt d = 0; d < cdim; ++d) pcoords[d] = PetscRealPart(coords[p * cdim + d]);
893: PetscCall(DMPlexCoordinatesToReference(dm, ctx->cells[p], 1, pcoords, xi));
894: PetscCall(DMPlexVecGetClosure(dm, NULL, x, ctx->cells[p], &clSize, &xa));
895: for (field = 0; field < Nf; ++field) {
896: PetscTabulation T;
897: PetscObject obj;
898: PetscClassId id;
900: PetscCall(PetscDSGetDiscretization(ds, field, &obj));
901: PetscCall(PetscObjectGetClassId(obj, &id));
902: if (id == PETSCFE_CLASSID) {
903: PetscFE fe = (PetscFE)obj;
905: PetscCall(PetscFECreateTabulation(fe, 1, 1, xi, 0, &T));
906: {
907: const PetscReal *basis = T->T[0];
908: const PetscInt Nb = T->Nb;
909: const PetscInt Nc = T->Nc;
911: for (PetscInt fc = 0; fc < Nc; ++fc) {
912: interpolant[p * ctx->dof + coff + fc] = 0.0;
913: for (PetscInt f = 0; f < Nb; ++f) interpolant[p * ctx->dof + coff + fc] += xa[foff + f] * basis[(0 * Nb + f) * Nc + fc];
914: }
915: coff += Nc;
916: foff += Nb;
917: }
918: PetscCall(PetscTabulationDestroy(&T));
919: } else if (id == PETSCFV_CLASSID) {
920: PetscFV fv = (PetscFV)obj;
921: PetscInt Nc;
923: // TODO Could use reconstruction if available
924: PetscCall(PetscFVGetNumComponents(fv, &Nc));
925: for (PetscInt fc = 0; fc < Nc; ++fc) interpolant[p * ctx->dof + coff + fc] = xa[foff + fc];
926: coff += Nc;
927: foff += Nc;
928: }
929: }
930: PetscCall(DMPlexVecRestoreClosure(dm, NULL, x, ctx->cells[p], &clSize, &xa));
931: PetscCheck(coff == ctx->dof, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Total components %" PetscInt_FMT " != %" PetscInt_FMT " dof specified for interpolation", coff, ctx->dof);
932: PetscCheck(foff == clSize, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Total FE/FV space size %" PetscInt_FMT " != %" PetscInt_FMT " closure size", foff, clSize);
933: }
934: PetscCall(VecRestoreArrayRead(ctx->coords, &coords));
935: PetscCall(VecRestoreArrayWrite(v, &interpolant));
936: } else {
937: for (PetscInt p = 0; p < ctx->n; ++p) {
938: const PetscInt cell = ctx->cells[p];
939: DMPolytopeType ct;
941: PetscCall(DMPlexGetCellType(dm, cell, &ct));
942: switch (ct) {
943: case DM_POLYTOPE_SEGMENT:
944: PetscCall(DMInterpolate_Segment_Private(ctx, dm, p, x, v));
945: break;
946: case DM_POLYTOPE_TRIANGLE:
947: PetscCall(DMInterpolate_Triangle_Private(ctx, dm, p, x, v));
948: break;
949: case DM_POLYTOPE_QUADRILATERAL:
950: PetscCall(DMInterpolate_Quad_Private(ctx, dm, p, x, v));
951: break;
952: case DM_POLYTOPE_TETRAHEDRON:
953: PetscCall(DMInterpolate_Tetrahedron_Private(ctx, dm, p, x, v));
954: break;
955: case DM_POLYTOPE_HEXAHEDRON:
956: PetscCall(DMInterpolate_Hex_Private(ctx, dm, cell, x, v));
957: break;
958: default:
959: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for cell type %s", DMPolytopeTypes[PetscMax(0, PetscMin(ct, DM_NUM_POLYTOPES))]);
960: }
961: }
962: }
963: PetscFunctionReturn(PETSC_SUCCESS);
964: }
966: /*@
967: DMInterpolationDestroy - Destroys a `DMInterpolationInfo` context
969: Collective
971: Input Parameter:
972: . ctx - the context
974: Level: beginner
976: .seealso: [](ch_dmbase), `DM`, `DMInterpolationInfo`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`
977: @*/
978: PetscErrorCode DMInterpolationDestroy(DMInterpolationInfo *ctx)
979: {
980: PetscFunctionBegin;
981: PetscAssertPointer(ctx, 1);
982: PetscCall(VecDestroy(&(*ctx)->coords));
983: PetscCall(PetscFree((*ctx)->points));
984: PetscCall(PetscFree((*ctx)->cells));
985: PetscCall(PetscFree(*ctx));
986: *ctx = NULL;
987: PetscFunctionReturn(PETSC_SUCCESS);
988: }