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