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, 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: PetscInt d;
260: for (d = 0; d < ctx->dim; ++d, ++i) a[i] = globalPoints[p * ctx->dim + d];
261: ctx->cells[q] = foundCells[q].index;
262: ++q;
263: }
264: if (globalProcs[p] == size && rank == 0) {
265: PetscInt d;
267: for (d = 0; d < ctx->dim; ++d, ++i) a[i] = 0.;
268: ctx->cells[q] = -1;
269: ++q;
270: }
271: }
272: PetscCall(VecRestoreArray(ctx->coords, &a));
273: #if 0
274: PetscCall(PetscFree3(foundCells,foundProcs,globalProcs));
275: #else
276: PetscCall(PetscFree2(foundProcs, globalProcs));
277: PetscCall(PetscSFDestroy(&cellSF));
278: PetscCall(VecDestroy(&pointVec));
279: #endif
280: if ((void *)globalPointsScalar != (void *)globalPoints) PetscCall(PetscFree(globalPointsScalar));
281: if (!redundantPoints) PetscCall(PetscFree3(globalPoints, counts, displs));
282: PetscCall(PetscLayoutDestroy(&layout));
283: PetscFunctionReturn(PETSC_SUCCESS);
284: }
286: /*@C
287: DMInterpolationGetCoordinates - Gets a `Vec` with the coordinates of each interpolation point
289: Collective
291: Input Parameter:
292: . ctx - the context
294: Output Parameter:
295: . coordinates - the coordinates of interpolation points
297: Level: intermediate
299: Note:
300: The local vector entries correspond to interpolation points lying on this process, according to the associated `DM`.
301: This is a borrowed vector that the user should not destroy.
303: .seealso: [](ch_dmbase), `DM`, `DMInterpolationInfo`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`
304: @*/
305: PetscErrorCode DMInterpolationGetCoordinates(DMInterpolationInfo ctx, Vec *coordinates)
306: {
307: PetscFunctionBegin;
308: PetscAssertPointer(coordinates, 2);
309: PetscCheck(ctx->coords, ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup.");
310: *coordinates = ctx->coords;
311: PetscFunctionReturn(PETSC_SUCCESS);
312: }
314: /*@C
315: DMInterpolationGetVector - Gets a `Vec` which can hold all the interpolated field values
317: Collective
319: Input Parameter:
320: . ctx - the context
322: Output Parameter:
323: . v - a vector capable of holding the interpolated field values
325: Level: intermediate
327: Note:
328: This vector should be returned using `DMInterpolationRestoreVector()`.
330: .seealso: [](ch_dmbase), `DM`, `DMInterpolationInfo`, `DMInterpolationRestoreVector()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`
331: @*/
332: PetscErrorCode DMInterpolationGetVector(DMInterpolationInfo ctx, Vec *v)
333: {
334: PetscFunctionBegin;
335: PetscAssertPointer(v, 2);
336: PetscCheck(ctx->coords, ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup.");
337: PetscCall(VecCreate(ctx->comm, v));
338: PetscCall(VecSetSizes(*v, ctx->n * ctx->dof, PETSC_DECIDE));
339: PetscCall(VecSetBlockSize(*v, ctx->dof));
340: PetscCall(VecSetType(*v, VECSTANDARD));
341: PetscFunctionReturn(PETSC_SUCCESS);
342: }
344: /*@C
345: DMInterpolationRestoreVector - Returns a `Vec` which can hold all the interpolated field values
347: Collective
349: Input Parameters:
350: + ctx - the context
351: - v - a vector capable of holding the interpolated field values
353: Level: intermediate
355: .seealso: [](ch_dmbase), `DM`, `DMInterpolationInfo`, `DMInterpolationGetVector()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`
356: @*/
357: PetscErrorCode DMInterpolationRestoreVector(DMInterpolationInfo ctx, Vec *v)
358: {
359: PetscFunctionBegin;
360: PetscAssertPointer(v, 2);
361: PetscCheck(ctx->coords, ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup.");
362: PetscCall(VecDestroy(v));
363: PetscFunctionReturn(PETSC_SUCCESS);
364: }
366: static inline PetscErrorCode DMInterpolate_Segment_Private(DMInterpolationInfo ctx, DM dm, PetscInt p, Vec xLocal, Vec v)
367: {
368: const PetscInt c = ctx->cells[p];
369: const PetscInt dof = ctx->dof;
370: PetscScalar *x = NULL;
371: const PetscScalar *coords;
372: PetscScalar *a;
373: PetscReal v0, J, invJ, detJ, xir[1];
374: PetscInt xSize, comp;
376: PetscFunctionBegin;
377: PetscCall(VecGetArrayRead(ctx->coords, &coords));
378: PetscCall(VecGetArray(v, &a));
379: PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, &v0, &J, &invJ, &detJ));
380: PetscCheck(detJ > 0.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %" PetscInt_FMT, (double)detJ, c);
381: xir[0] = invJ * PetscRealPart(coords[p] - v0);
382: PetscCall(DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x));
383: if (2 * dof == xSize) {
384: for (comp = 0; comp < dof; ++comp) a[p * dof + comp] = x[0 * dof + comp] * (1 - xir[0]) + x[1 * dof + comp] * xir[0];
385: } else if (dof == xSize) {
386: for (comp = 0; comp < dof; ++comp) a[p * dof + comp] = x[0 * dof + comp];
387: } 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);
388: PetscCall(DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x));
389: PetscCall(VecRestoreArray(v, &a));
390: PetscCall(VecRestoreArrayRead(ctx->coords, &coords));
391: PetscFunctionReturn(PETSC_SUCCESS);
392: }
394: static inline PetscErrorCode DMInterpolate_Triangle_Private(DMInterpolationInfo ctx, DM dm, PetscInt p, Vec xLocal, Vec v)
395: {
396: const PetscInt c = ctx->cells[p];
397: PetscScalar *x = NULL;
398: const PetscScalar *coords;
399: PetscScalar *a;
400: PetscReal *v0, *J, *invJ, detJ;
401: PetscReal xi[4];
403: PetscFunctionBegin;
404: PetscCall(PetscMalloc3(ctx->dim, &v0, ctx->dim * ctx->dim, &J, ctx->dim * ctx->dim, &invJ));
405: PetscCall(VecGetArrayRead(ctx->coords, &coords));
406: PetscCall(VecGetArray(v, &a));
407: PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ));
408: PetscCheck(detJ > 0.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %" PetscInt_FMT, (double)detJ, c);
409: PetscCall(DMPlexVecGetClosure(dm, NULL, xLocal, c, NULL, &x));
410: for (PetscInt comp = 0; comp < ctx->dof; ++comp) a[p * ctx->dof + comp] = x[0 * ctx->dof + comp];
411: for (PetscInt d = 0; d < ctx->dim; ++d) {
412: xi[d] = 0.0;
413: 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]);
414: 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];
415: }
416: PetscCall(DMPlexVecRestoreClosure(dm, NULL, xLocal, c, NULL, &x));
417: PetscCall(VecRestoreArray(v, &a));
418: PetscCall(VecRestoreArrayRead(ctx->coords, &coords));
419: PetscCall(PetscFree3(v0, J, invJ));
420: PetscFunctionReturn(PETSC_SUCCESS);
421: }
423: static inline PetscErrorCode DMInterpolate_Tetrahedron_Private(DMInterpolationInfo ctx, DM dm, PetscInt p, Vec xLocal, Vec v)
424: {
425: const PetscInt c = ctx->cells[p];
426: const PetscInt order[3] = {2, 1, 3};
427: PetscScalar *x = NULL;
428: PetscReal *v0, *J, *invJ, detJ;
429: const PetscScalar *coords;
430: PetscScalar *a;
431: PetscReal xi[4];
433: PetscFunctionBegin;
434: PetscCall(PetscMalloc3(ctx->dim, &v0, ctx->dim * ctx->dim, &J, ctx->dim * ctx->dim, &invJ));
435: PetscCall(VecGetArrayRead(ctx->coords, &coords));
436: PetscCall(VecGetArray(v, &a));
437: PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ));
438: PetscCheck(detJ > 0.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %" PetscInt_FMT, (double)detJ, c);
439: PetscCall(DMPlexVecGetClosure(dm, NULL, xLocal, c, NULL, &x));
440: for (PetscInt comp = 0; comp < ctx->dof; ++comp) a[p * ctx->dof + comp] = x[0 * ctx->dof + comp];
441: for (PetscInt d = 0; d < ctx->dim; ++d) {
442: xi[d] = 0.0;
443: 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]);
444: 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];
445: }
446: PetscCall(DMPlexVecRestoreClosure(dm, NULL, xLocal, c, NULL, &x));
447: PetscCall(VecRestoreArray(v, &a));
448: PetscCall(VecRestoreArrayRead(ctx->coords, &coords));
449: PetscCall(PetscFree3(v0, J, invJ));
450: PetscFunctionReturn(PETSC_SUCCESS);
451: }
453: static inline PetscErrorCode QuadMap_Private(SNES snes, Vec Xref, Vec Xreal, void *ctx)
454: {
455: const PetscScalar *vertices = (const PetscScalar *)ctx;
456: const PetscScalar x0 = vertices[0];
457: const PetscScalar y0 = vertices[1];
458: const PetscScalar x1 = vertices[2];
459: const PetscScalar y1 = vertices[3];
460: const PetscScalar x2 = vertices[4];
461: const PetscScalar y2 = vertices[5];
462: const PetscScalar x3 = vertices[6];
463: const PetscScalar y3 = vertices[7];
464: const PetscScalar f_1 = x1 - x0;
465: const PetscScalar g_1 = y1 - y0;
466: const PetscScalar f_3 = x3 - x0;
467: const PetscScalar g_3 = y3 - y0;
468: const PetscScalar f_01 = x2 - x1 - x3 + x0;
469: const PetscScalar g_01 = y2 - y1 - y3 + y0;
470: const PetscScalar *ref;
471: PetscScalar *real;
473: PetscFunctionBegin;
474: PetscCall(VecGetArrayRead(Xref, &ref));
475: PetscCall(VecGetArray(Xreal, &real));
476: {
477: const PetscScalar p0 = ref[0];
478: const PetscScalar p1 = ref[1];
480: real[0] = x0 + f_1 * p0 + f_3 * p1 + f_01 * p0 * p1;
481: real[1] = y0 + g_1 * p0 + g_3 * p1 + g_01 * p0 * p1;
482: }
483: PetscCall(PetscLogFlops(28));
484: PetscCall(VecRestoreArrayRead(Xref, &ref));
485: PetscCall(VecRestoreArray(Xreal, &real));
486: PetscFunctionReturn(PETSC_SUCCESS);
487: }
489: #include <petsc/private/dmimpl.h>
490: static inline PetscErrorCode QuadJacobian_Private(SNES snes, Vec Xref, Mat J, Mat M, void *ctx)
491: {
492: const PetscScalar *vertices = (const PetscScalar *)ctx;
493: const PetscScalar x0 = vertices[0];
494: const PetscScalar y0 = vertices[1];
495: const PetscScalar x1 = vertices[2];
496: const PetscScalar y1 = vertices[3];
497: const PetscScalar x2 = vertices[4];
498: const PetscScalar y2 = vertices[5];
499: const PetscScalar x3 = vertices[6];
500: const PetscScalar y3 = vertices[7];
501: const PetscScalar f_01 = x2 - x1 - x3 + x0;
502: const PetscScalar g_01 = y2 - y1 - y3 + y0;
503: const PetscScalar *ref;
505: PetscFunctionBegin;
506: PetscCall(VecGetArrayRead(Xref, &ref));
507: {
508: const PetscScalar x = ref[0];
509: const PetscScalar y = ref[1];
510: const PetscInt rows[2] = {0, 1};
511: PetscScalar values[4];
513: values[0] = (x1 - x0 + f_01 * y) * 0.5;
514: values[1] = (x3 - x0 + f_01 * x) * 0.5;
515: values[2] = (y1 - y0 + g_01 * y) * 0.5;
516: values[3] = (y3 - y0 + g_01 * x) * 0.5;
517: PetscCall(MatSetValues(J, 2, rows, 2, rows, values, INSERT_VALUES));
518: }
519: PetscCall(PetscLogFlops(30));
520: PetscCall(VecRestoreArrayRead(Xref, &ref));
521: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
522: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
523: PetscFunctionReturn(PETSC_SUCCESS);
524: }
526: static inline PetscErrorCode DMInterpolate_Quad_Private(DMInterpolationInfo ctx, DM dm, PetscInt p, Vec xLocal, Vec v)
527: {
528: const PetscInt c = ctx->cells[p];
529: PetscFE fem = NULL;
530: DM dmCoord;
531: SNES snes;
532: KSP ksp;
533: PC pc;
534: Vec coordsLocal, r, ref, real;
535: Mat J;
536: PetscTabulation T = NULL;
537: const PetscScalar *coords;
538: PetscScalar *a;
539: PetscReal xir[2] = {0., 0.};
540: PetscInt Nf;
541: const PetscInt dof = ctx->dof;
542: PetscScalar *x = NULL, *vertices = NULL;
543: PetscScalar *xi;
544: PetscInt coordSize, xSize;
546: PetscFunctionBegin;
547: PetscCall(DMGetNumFields(dm, &Nf));
548: if (Nf) {
549: PetscObject obj;
550: PetscClassId id;
552: PetscCall(DMGetField(dm, 0, NULL, &obj));
553: PetscCall(PetscObjectGetClassId(obj, &id));
554: if (id == PETSCFE_CLASSID) {
555: fem = (PetscFE)obj;
556: PetscCall(PetscFECreateTabulation(fem, 1, 1, xir, 0, &T));
557: }
558: }
559: PetscCall(DMGetCoordinatesLocal(dm, &coordsLocal));
560: PetscCall(DMGetCoordinateDM(dm, &dmCoord));
561: PetscCall(SNESCreate(PETSC_COMM_SELF, &snes));
562: PetscCall(SNESSetOptionsPrefix(snes, "quad_interp_"));
563: PetscCall(VecCreate(PETSC_COMM_SELF, &r));
564: PetscCall(VecSetSizes(r, 2, 2));
565: PetscCall(VecSetType(r, dm->vectype));
566: PetscCall(VecDuplicate(r, &ref));
567: PetscCall(VecDuplicate(r, &real));
568: PetscCall(MatCreate(PETSC_COMM_SELF, &J));
569: PetscCall(MatSetSizes(J, 2, 2, 2, 2));
570: PetscCall(MatSetType(J, MATSEQDENSE));
571: PetscCall(MatSetUp(J));
572: PetscCall(SNESSetFunction(snes, r, QuadMap_Private, NULL));
573: PetscCall(SNESSetJacobian(snes, J, J, QuadJacobian_Private, NULL));
574: PetscCall(SNESGetKSP(snes, &ksp));
575: PetscCall(KSPGetPC(ksp, &pc));
576: PetscCall(PCSetType(pc, PCLU));
577: PetscCall(SNESSetFromOptions(snes));
579: PetscCall(VecGetArrayRead(ctx->coords, &coords));
580: PetscCall(VecGetArray(v, &a));
581: PetscCall(DMPlexVecGetClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices));
582: PetscCheck(4 * 2 == coordSize, ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid closure size %" PetscInt_FMT " should be %d", coordSize, 4 * 2);
583: PetscCall(DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x));
584: PetscCall(SNESSetFunction(snes, NULL, NULL, vertices));
585: PetscCall(SNESSetJacobian(snes, NULL, NULL, NULL, vertices));
586: PetscCall(VecGetArray(real, &xi));
587: xi[0] = coords[p * ctx->dim + 0];
588: xi[1] = coords[p * ctx->dim + 1];
589: PetscCall(VecRestoreArray(real, &xi));
590: PetscCall(SNESSolve(snes, real, ref));
591: PetscCall(VecGetArray(ref, &xi));
592: xir[0] = PetscRealPart(xi[0]);
593: xir[1] = PetscRealPart(xi[1]);
594: if (4 * dof == xSize) {
595: 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];
596: } else if (dof == xSize) {
597: for (PetscInt comp = 0; comp < dof; ++comp) a[p * dof + comp] = x[0 * dof + comp];
598: } else {
599: PetscCheck(fem, ctx->comm, PETSC_ERR_ARG_WRONG, "Cannot have a higher order interpolant if the discretization is not PetscFE");
600: xir[0] = 2.0 * xir[0] - 1.0;
601: xir[1] = 2.0 * xir[1] - 1.0;
602: PetscCall(PetscFEComputeTabulation(fem, 1, xir, 0, T));
603: for (PetscInt comp = 0; comp < dof; ++comp) {
604: a[p * dof + comp] = 0.0;
605: for (PetscInt d = 0; d < xSize / dof; ++d) a[p * dof + comp] += x[d * dof + comp] * T->T[0][d * dof + comp];
606: }
607: }
608: PetscCall(VecRestoreArray(ref, &xi));
609: PetscCall(DMPlexVecRestoreClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices));
610: PetscCall(DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x));
611: PetscCall(PetscTabulationDestroy(&T));
612: PetscCall(VecRestoreArray(v, &a));
613: PetscCall(VecRestoreArrayRead(ctx->coords, &coords));
615: PetscCall(SNESDestroy(&snes));
616: PetscCall(VecDestroy(&r));
617: PetscCall(VecDestroy(&ref));
618: PetscCall(VecDestroy(&real));
619: PetscCall(MatDestroy(&J));
620: PetscFunctionReturn(PETSC_SUCCESS);
621: }
623: static inline PetscErrorCode HexMap_Private(SNES snes, Vec Xref, Vec Xreal, void *ctx)
624: {
625: const PetscScalar *vertices = (const PetscScalar *)ctx;
626: const PetscScalar x0 = vertices[0];
627: const PetscScalar y0 = vertices[1];
628: const PetscScalar z0 = vertices[2];
629: const PetscScalar x1 = vertices[9];
630: const PetscScalar y1 = vertices[10];
631: const PetscScalar z1 = vertices[11];
632: const PetscScalar x2 = vertices[6];
633: const PetscScalar y2 = vertices[7];
634: const PetscScalar z2 = vertices[8];
635: const PetscScalar x3 = vertices[3];
636: const PetscScalar y3 = vertices[4];
637: const PetscScalar z3 = vertices[5];
638: const PetscScalar x4 = vertices[12];
639: const PetscScalar y4 = vertices[13];
640: const PetscScalar z4 = vertices[14];
641: const PetscScalar x5 = vertices[15];
642: const PetscScalar y5 = vertices[16];
643: const PetscScalar z5 = vertices[17];
644: const PetscScalar x6 = vertices[18];
645: const PetscScalar y6 = vertices[19];
646: const PetscScalar z6 = vertices[20];
647: const PetscScalar x7 = vertices[21];
648: const PetscScalar y7 = vertices[22];
649: const PetscScalar z7 = vertices[23];
650: const PetscScalar f_1 = x1 - x0;
651: const PetscScalar g_1 = y1 - y0;
652: const PetscScalar h_1 = z1 - z0;
653: const PetscScalar f_3 = x3 - x0;
654: const PetscScalar g_3 = y3 - y0;
655: const PetscScalar h_3 = z3 - z0;
656: const PetscScalar f_4 = x4 - x0;
657: const PetscScalar g_4 = y4 - y0;
658: const PetscScalar h_4 = z4 - z0;
659: const PetscScalar f_01 = x2 - x1 - x3 + x0;
660: const PetscScalar g_01 = y2 - y1 - y3 + y0;
661: const PetscScalar h_01 = z2 - z1 - z3 + z0;
662: const PetscScalar f_12 = x7 - x3 - x4 + x0;
663: const PetscScalar g_12 = y7 - y3 - y4 + y0;
664: const PetscScalar h_12 = z7 - z3 - z4 + z0;
665: const PetscScalar f_02 = x5 - x1 - x4 + x0;
666: const PetscScalar g_02 = y5 - y1 - y4 + y0;
667: const PetscScalar h_02 = z5 - z1 - z4 + z0;
668: const PetscScalar f_012 = x6 - x0 + x1 - x2 + x3 + x4 - x5 - x7;
669: const PetscScalar g_012 = y6 - y0 + y1 - y2 + y3 + y4 - y5 - y7;
670: const PetscScalar h_012 = z6 - z0 + z1 - z2 + z3 + z4 - z5 - z7;
671: const PetscScalar *ref;
672: PetscScalar *real;
674: PetscFunctionBegin;
675: PetscCall(VecGetArrayRead(Xref, &ref));
676: PetscCall(VecGetArray(Xreal, &real));
677: {
678: const PetscScalar p0 = ref[0];
679: const PetscScalar p1 = ref[1];
680: const PetscScalar p2 = ref[2];
682: 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;
683: 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;
684: 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;
685: }
686: PetscCall(PetscLogFlops(114));
687: PetscCall(VecRestoreArrayRead(Xref, &ref));
688: PetscCall(VecRestoreArray(Xreal, &real));
689: PetscFunctionReturn(PETSC_SUCCESS);
690: }
692: static inline PetscErrorCode HexJacobian_Private(SNES snes, Vec Xref, Mat J, Mat M, void *ctx)
693: {
694: const PetscScalar *vertices = (const PetscScalar *)ctx;
695: const PetscScalar x0 = vertices[0];
696: const PetscScalar y0 = vertices[1];
697: const PetscScalar z0 = vertices[2];
698: const PetscScalar x1 = vertices[9];
699: const PetscScalar y1 = vertices[10];
700: const PetscScalar z1 = vertices[11];
701: const PetscScalar x2 = vertices[6];
702: const PetscScalar y2 = vertices[7];
703: const PetscScalar z2 = vertices[8];
704: const PetscScalar x3 = vertices[3];
705: const PetscScalar y3 = vertices[4];
706: const PetscScalar z3 = vertices[5];
707: const PetscScalar x4 = vertices[12];
708: const PetscScalar y4 = vertices[13];
709: const PetscScalar z4 = vertices[14];
710: const PetscScalar x5 = vertices[15];
711: const PetscScalar y5 = vertices[16];
712: const PetscScalar z5 = vertices[17];
713: const PetscScalar x6 = vertices[18];
714: const PetscScalar y6 = vertices[19];
715: const PetscScalar z6 = vertices[20];
716: const PetscScalar x7 = vertices[21];
717: const PetscScalar y7 = vertices[22];
718: const PetscScalar z7 = vertices[23];
719: const PetscScalar f_xy = x2 - x1 - x3 + x0;
720: const PetscScalar g_xy = y2 - y1 - y3 + y0;
721: const PetscScalar h_xy = z2 - z1 - z3 + z0;
722: const PetscScalar f_yz = x7 - x3 - x4 + x0;
723: const PetscScalar g_yz = y7 - y3 - y4 + y0;
724: const PetscScalar h_yz = z7 - z3 - z4 + z0;
725: const PetscScalar f_xz = x5 - x1 - x4 + x0;
726: const PetscScalar g_xz = y5 - y1 - y4 + y0;
727: const PetscScalar h_xz = z5 - z1 - z4 + z0;
728: const PetscScalar f_xyz = x6 - x0 + x1 - x2 + x3 + x4 - x5 - x7;
729: const PetscScalar g_xyz = y6 - y0 + y1 - y2 + y3 + y4 - y5 - y7;
730: const PetscScalar h_xyz = z6 - z0 + z1 - z2 + z3 + z4 - z5 - z7;
731: const PetscScalar *ref;
733: PetscFunctionBegin;
734: PetscCall(VecGetArrayRead(Xref, &ref));
735: {
736: const PetscScalar x = ref[0];
737: const PetscScalar y = ref[1];
738: const PetscScalar z = ref[2];
739: const PetscInt rows[3] = {0, 1, 2};
740: PetscScalar values[9];
742: values[0] = (x1 - x0 + f_xy * y + f_xz * z + f_xyz * y * z) / 2.0;
743: values[1] = (x3 - x0 + f_xy * x + f_yz * z + f_xyz * x * z) / 2.0;
744: values[2] = (x4 - x0 + f_yz * y + f_xz * x + f_xyz * x * y) / 2.0;
745: values[3] = (y1 - y0 + g_xy * y + g_xz * z + g_xyz * y * z) / 2.0;
746: values[4] = (y3 - y0 + g_xy * x + g_yz * z + g_xyz * x * z) / 2.0;
747: values[5] = (y4 - y0 + g_yz * y + g_xz * x + g_xyz * x * y) / 2.0;
748: values[6] = (z1 - z0 + h_xy * y + h_xz * z + h_xyz * y * z) / 2.0;
749: values[7] = (z3 - z0 + h_xy * x + h_yz * z + h_xyz * x * z) / 2.0;
750: values[8] = (z4 - z0 + h_yz * y + h_xz * x + h_xyz * x * y) / 2.0;
752: PetscCall(MatSetValues(J, 3, rows, 3, rows, values, INSERT_VALUES));
753: }
754: PetscCall(PetscLogFlops(152));
755: PetscCall(VecRestoreArrayRead(Xref, &ref));
756: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
757: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
758: PetscFunctionReturn(PETSC_SUCCESS);
759: }
761: static inline PetscErrorCode DMInterpolate_Hex_Private(DMInterpolationInfo ctx, DM dm, PetscInt p, Vec xLocal, Vec v)
762: {
763: const PetscInt c = ctx->cells[p];
764: DM dmCoord;
765: SNES snes;
766: KSP ksp;
767: PC pc;
768: Vec coordsLocal, r, ref, real;
769: Mat J;
770: const PetscScalar *coords;
771: PetscScalar *a;
772: PetscScalar *x = NULL, *vertices = NULL;
773: PetscScalar *xi;
774: PetscReal xir[3];
775: PetscInt coordSize, xSize;
777: PetscFunctionBegin;
778: PetscCall(DMGetCoordinatesLocal(dm, &coordsLocal));
779: PetscCall(DMGetCoordinateDM(dm, &dmCoord));
780: PetscCall(SNESCreate(PETSC_COMM_SELF, &snes));
781: PetscCall(SNESSetOptionsPrefix(snes, "hex_interp_"));
782: PetscCall(VecCreate(PETSC_COMM_SELF, &r));
783: PetscCall(VecSetSizes(r, 3, 3));
784: PetscCall(VecSetType(r, dm->vectype));
785: PetscCall(VecDuplicate(r, &ref));
786: PetscCall(VecDuplicate(r, &real));
787: PetscCall(MatCreate(PETSC_COMM_SELF, &J));
788: PetscCall(MatSetSizes(J, 3, 3, 3, 3));
789: PetscCall(MatSetType(J, MATSEQDENSE));
790: PetscCall(MatSetUp(J));
791: PetscCall(SNESSetFunction(snes, r, HexMap_Private, NULL));
792: PetscCall(SNESSetJacobian(snes, J, J, HexJacobian_Private, NULL));
793: PetscCall(SNESGetKSP(snes, &ksp));
794: PetscCall(KSPGetPC(ksp, &pc));
795: PetscCall(PCSetType(pc, PCLU));
796: PetscCall(SNESSetFromOptions(snes));
798: PetscCall(VecGetArrayRead(ctx->coords, &coords));
799: PetscCall(VecGetArray(v, &a));
800: PetscCall(DMPlexVecGetClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices));
801: PetscCheck(8 * 3 == coordSize, ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid coordinate closure size %" PetscInt_FMT " should be %d", coordSize, 8 * 3);
802: PetscCall(DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x));
803: 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);
804: PetscCall(SNESSetFunction(snes, NULL, NULL, vertices));
805: PetscCall(SNESSetJacobian(snes, NULL, NULL, NULL, vertices));
806: PetscCall(VecGetArray(real, &xi));
807: xi[0] = coords[p * ctx->dim + 0];
808: xi[1] = coords[p * ctx->dim + 1];
809: xi[2] = coords[p * ctx->dim + 2];
810: PetscCall(VecRestoreArray(real, &xi));
811: PetscCall(SNESSolve(snes, real, ref));
812: PetscCall(VecGetArray(ref, &xi));
813: xir[0] = PetscRealPart(xi[0]);
814: xir[1] = PetscRealPart(xi[1]);
815: xir[2] = PetscRealPart(xi[2]);
816: if (8 * ctx->dof == xSize) {
817: for (PetscInt comp = 0; comp < ctx->dof; ++comp) {
818: 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]) +
819: 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];
820: }
821: } else {
822: for (PetscInt comp = 0; comp < ctx->dof; ++comp) a[p * ctx->dof + comp] = x[0 * ctx->dof + comp];
823: }
824: PetscCall(VecRestoreArray(ref, &xi));
825: PetscCall(DMPlexVecRestoreClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices));
826: PetscCall(DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x));
827: PetscCall(VecRestoreArray(v, &a));
828: PetscCall(VecRestoreArrayRead(ctx->coords, &coords));
830: PetscCall(SNESDestroy(&snes));
831: PetscCall(VecDestroy(&r));
832: PetscCall(VecDestroy(&ref));
833: PetscCall(VecDestroy(&real));
834: PetscCall(MatDestroy(&J));
835: PetscFunctionReturn(PETSC_SUCCESS);
836: }
838: /*@C
839: DMInterpolationEvaluate - Using the input from `dm` and `x`, calculates interpolated field values at the interpolation points.
841: Input Parameters:
842: + ctx - The `DMInterpolationInfo` context obtained with `DMInterpolationCreate()`
843: . dm - The `DM`
844: - x - The local vector containing the field to be interpolated, can be created with `DMCreateGlobalVector()`
846: Output Parameter:
847: . v - The vector containing the interpolated values, obtained with `DMInterpolationGetVector()`
849: Level: beginner
851: .seealso: [](ch_dmbase), `DM`, `DMInterpolationInfo`, `DMInterpolationGetVector()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`, `DMInterpolationGetCoordinates()`
852: @*/
853: PetscErrorCode DMInterpolationEvaluate(DMInterpolationInfo ctx, DM dm, Vec x, Vec v)
854: {
855: PetscDS ds;
856: PetscInt n, p, Nf, field;
857: PetscBool useDS = PETSC_FALSE;
859: PetscFunctionBegin;
863: PetscCall(VecGetLocalSize(v, &n));
864: 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);
865: if (!n) PetscFunctionReturn(PETSC_SUCCESS);
866: PetscCall(DMGetDS(dm, &ds));
867: if (ds) {
868: useDS = PETSC_TRUE;
869: PetscCall(PetscDSGetNumFields(ds, &Nf));
870: for (field = 0; field < Nf; ++field) {
871: PetscObject obj;
872: PetscClassId id;
874: PetscCall(PetscDSGetDiscretization(ds, field, &obj));
875: PetscCall(PetscObjectGetClassId(obj, &id));
876: if (id != PETSCFE_CLASSID && id != PETSCFV_CLASSID) {
877: useDS = PETSC_FALSE;
878: break;
879: }
880: }
881: }
882: if (useDS) {
883: const PetscScalar *coords;
884: PetscScalar *interpolant;
885: PetscInt cdim, d;
887: PetscCall(DMGetCoordinateDim(dm, &cdim));
888: PetscCall(VecGetArrayRead(ctx->coords, &coords));
889: PetscCall(VecGetArrayWrite(v, &interpolant));
890: for (p = 0; p < ctx->n; ++p) {
891: PetscReal pcoords[3], xi[3];
892: PetscScalar *xa = NULL;
893: PetscInt coff = 0, foff = 0, clSize;
895: if (ctx->cells[p] < 0) continue;
896: for (d = 0; d < cdim; ++d) pcoords[d] = PetscRealPart(coords[p * cdim + d]);
897: PetscCall(DMPlexCoordinatesToReference(dm, ctx->cells[p], 1, pcoords, xi));
898: PetscCall(DMPlexVecGetClosure(dm, NULL, x, ctx->cells[p], &clSize, &xa));
899: for (field = 0; field < Nf; ++field) {
900: PetscTabulation T;
901: PetscObject obj;
902: PetscClassId id;
904: PetscCall(PetscDSGetDiscretization(ds, field, &obj));
905: PetscCall(PetscObjectGetClassId(obj, &id));
906: if (id == PETSCFE_CLASSID) {
907: PetscFE fe = (PetscFE)obj;
909: PetscCall(PetscFECreateTabulation(fe, 1, 1, xi, 0, &T));
910: {
911: const PetscReal *basis = T->T[0];
912: const PetscInt Nb = T->Nb;
913: const PetscInt Nc = T->Nc;
915: for (PetscInt fc = 0; fc < Nc; ++fc) {
916: interpolant[p * ctx->dof + coff + fc] = 0.0;
917: for (PetscInt f = 0; f < Nb; ++f) interpolant[p * ctx->dof + coff + fc] += xa[foff + f] * basis[(0 * Nb + f) * Nc + fc];
918: }
919: coff += Nc;
920: foff += Nb;
921: }
922: PetscCall(PetscTabulationDestroy(&T));
923: } else if (id == PETSCFV_CLASSID) {
924: PetscFV fv = (PetscFV)obj;
925: PetscInt Nc;
927: // TODO Could use reconstruction if available
928: PetscCall(PetscFVGetNumComponents(fv, &Nc));
929: for (PetscInt fc = 0; fc < Nc; ++fc) interpolant[p * ctx->dof + coff + fc] = xa[foff + fc];
930: coff += Nc;
931: foff += Nc;
932: }
933: }
934: PetscCall(DMPlexVecRestoreClosure(dm, NULL, x, ctx->cells[p], &clSize, &xa));
935: PetscCheck(coff == ctx->dof, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Total components %" PetscInt_FMT " != %" PetscInt_FMT " dof specified for interpolation", coff, ctx->dof);
936: PetscCheck(foff == clSize, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Total FE/FV space size %" PetscInt_FMT " != %" PetscInt_FMT " closure size", foff, clSize);
937: }
938: PetscCall(VecRestoreArrayRead(ctx->coords, &coords));
939: PetscCall(VecRestoreArrayWrite(v, &interpolant));
940: } else {
941: for (PetscInt p = 0; p < ctx->n; ++p) {
942: const PetscInt cell = ctx->cells[p];
943: DMPolytopeType ct;
945: PetscCall(DMPlexGetCellType(dm, cell, &ct));
946: switch (ct) {
947: case DM_POLYTOPE_SEGMENT:
948: PetscCall(DMInterpolate_Segment_Private(ctx, dm, p, x, v));
949: break;
950: case DM_POLYTOPE_TRIANGLE:
951: PetscCall(DMInterpolate_Triangle_Private(ctx, dm, p, x, v));
952: break;
953: case DM_POLYTOPE_QUADRILATERAL:
954: PetscCall(DMInterpolate_Quad_Private(ctx, dm, p, x, v));
955: break;
956: case DM_POLYTOPE_TETRAHEDRON:
957: PetscCall(DMInterpolate_Tetrahedron_Private(ctx, dm, p, x, v));
958: break;
959: case DM_POLYTOPE_HEXAHEDRON:
960: PetscCall(DMInterpolate_Hex_Private(ctx, dm, cell, x, v));
961: break;
962: default:
963: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for cell type %s", DMPolytopeTypes[PetscMax(0, PetscMin(ct, DM_NUM_POLYTOPES))]);
964: }
965: }
966: }
967: PetscFunctionReturn(PETSC_SUCCESS);
968: }
970: /*@C
971: DMInterpolationDestroy - Destroys a `DMInterpolationInfo` context
973: Collective
975: Input Parameter:
976: . ctx - the context
978: Level: beginner
980: .seealso: [](ch_dmbase), `DM`, `DMInterpolationInfo`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`
981: @*/
982: PetscErrorCode DMInterpolationDestroy(DMInterpolationInfo *ctx)
983: {
984: PetscFunctionBegin;
985: PetscAssertPointer(ctx, 1);
986: PetscCall(VecDestroy(&(*ctx)->coords));
987: PetscCall(PetscFree((*ctx)->points));
988: PetscCall(PetscFree((*ctx)->cells));
989: PetscCall(PetscFree(*ctx));
990: *ctx = NULL;
991: PetscFunctionReturn(PETSC_SUCCESS);
992: }