1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/
2 #include <petsc/private/snesimpl.h> /*I "petscsnes.h" I*/
3 #include <petscds.h>
4 #include <petsc/private/petscimpl.h>
5 #include <petsc/private/petscfeimpl.h>
6
7 /*@
8 DMInterpolationCreate - Creates a `DMInterpolationInfo` context
9
10 Collective
11
12 Input Parameter:
13 . comm - the communicator
14
15 Output Parameter:
16 . ctx - the context
17
18 Level: beginner
19
20 Developer Note:
21 The naming is incorrect, either the object should be named `DMInterpolation` or all the routines should begin with `DMInterpolationInfo`
22
23 .seealso: [](ch_dmbase), `DM`, `DMInterpolationInfo`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationDestroy()`
24 @*/
DMInterpolationCreate(MPI_Comm comm,DMInterpolationInfo * ctx)25 PetscErrorCode DMInterpolationCreate(MPI_Comm comm, DMInterpolationInfo *ctx)
26 {
27 PetscFunctionBegin;
28 PetscAssertPointer(ctx, 2);
29 PetscCall(PetscNew(ctx));
30
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 }
40
41 /*@
42 DMInterpolationSetDim - Sets the spatial dimension for the interpolation context
43
44 Not Collective
45
46 Input Parameters:
47 + ctx - the context
48 - dim - the spatial dimension
49
50 Level: intermediate
51
52 .seealso: [](ch_dmbase), `DM`, `DMInterpolationInfo`, `DMInterpolationGetDim()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`
53 @*/
DMInterpolationSetDim(DMInterpolationInfo ctx,PetscInt dim)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 }
61
62 /*@
63 DMInterpolationGetDim - Gets the spatial dimension for the interpolation context
64
65 Not Collective
66
67 Input Parameter:
68 . ctx - the context
69
70 Output Parameter:
71 . dim - the spatial dimension
72
73 Level: intermediate
74
75 .seealso: [](ch_dmbase), `DM`, `DMInterpolationInfo`, `DMInterpolationSetDim()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`
76 @*/
DMInterpolationGetDim(DMInterpolationInfo ctx,PetscInt * dim)77 PetscErrorCode DMInterpolationGetDim(DMInterpolationInfo ctx, PetscInt *dim)
78 {
79 PetscFunctionBegin;
80 PetscAssertPointer(dim, 2);
81 *dim = ctx->dim;
82 PetscFunctionReturn(PETSC_SUCCESS);
83 }
84
85 /*@
86 DMInterpolationSetDof - Sets the number of fields interpolated at a point for the interpolation context
87
88 Not Collective
89
90 Input Parameters:
91 + ctx - the context
92 - dof - the number of fields
93
94 Level: intermediate
95
96 .seealso: [](ch_dmbase), `DM`, `DMInterpolationInfo`, `DMInterpolationGetDof()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`
97 @*/
DMInterpolationSetDof(DMInterpolationInfo ctx,PetscInt dof)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 }
105
106 /*@
107 DMInterpolationGetDof - Gets the number of fields interpolated at a point for the interpolation context
108
109 Not Collective
110
111 Input Parameter:
112 . ctx - the context
113
114 Output Parameter:
115 . dof - the number of fields
116
117 Level: intermediate
118
119 .seealso: [](ch_dmbase), `DM`, `DMInterpolationInfo`, `DMInterpolationSetDof()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`
120 @*/
DMInterpolationGetDof(DMInterpolationInfo ctx,PetscInt * dof)121 PetscErrorCode DMInterpolationGetDof(DMInterpolationInfo ctx, PetscInt *dof)
122 {
123 PetscFunctionBegin;
124 PetscAssertPointer(dof, 2);
125 *dof = ctx->dof;
126 PetscFunctionReturn(PETSC_SUCCESS);
127 }
128
129 /*@
130 DMInterpolationAddPoints - Add points at which we will interpolate the fields
131
132 Not Collective
133
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
138
139 Level: intermediate
140
141 Note:
142 The input coordinate information is copied into the object.
143
144 .seealso: [](ch_dmbase), `DM`, `DMInterpolationInfo`, `DMInterpolationSetDim()`, `DMInterpolationEvaluate()`, `DMInterpolationCreate()`
145 @*/
DMInterpolationAddPoints(DMInterpolationInfo ctx,PetscInt n,PetscReal points[])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;
152
153 PetscCall(PetscMalloc1(n * ctx->dim, &ctx->points));
154 PetscCall(PetscArraycpy(ctx->points, points, n * ctx->dim));
155 PetscFunctionReturn(PETSC_SUCCESS);
156 }
157
158 /*@
159 DMInterpolationSetUp - Compute spatial indices for point location during interpolation
160
161 Collective
162
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
168
169 Level: intermediate
170
171 .seealso: [](ch_dmbase), `DM`, `DMInterpolationInfo`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`
172 @*/
DMInterpolationSetUp(DMInterpolationInfo ctx,DM dm,PetscBool redundantPoints,PetscBool ignoreOutsideDomain)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;
190
191 PetscFunctionBegin;
192 PetscValidHeaderSpecific(dm, DM_CLASSID, 2);
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;
259
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;
266
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 }
285
286 /*@
287 DMInterpolationGetCoordinates - Gets a `Vec` with the coordinates of each interpolation point
288
289 Collective
290
291 Input Parameter:
292 . ctx - the context
293
294 Output Parameter:
295 . coordinates - the coordinates of interpolation points
296
297 Level: intermediate
298
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.
302
303 .seealso: [](ch_dmbase), `DM`, `DMInterpolationInfo`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`
304 @*/
DMInterpolationGetCoordinates(DMInterpolationInfo ctx,Vec * coordinates)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 }
313
314 /*@
315 DMInterpolationGetVector - Gets a `Vec` which can hold all the interpolated field values
316
317 Collective
318
319 Input Parameter:
320 . ctx - the context
321
322 Output Parameter:
323 . v - a vector capable of holding the interpolated field values
324
325 Level: intermediate
326
327 Note:
328 This vector should be returned using `DMInterpolationRestoreVector()`.
329
330 .seealso: [](ch_dmbase), `DM`, `DMInterpolationInfo`, `DMInterpolationRestoreVector()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`
331 @*/
DMInterpolationGetVector(DMInterpolationInfo ctx,Vec * v)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 }
343
344 /*@
345 DMInterpolationRestoreVector - Returns a `Vec` which can hold all the interpolated field values
346
347 Collective
348
349 Input Parameters:
350 + ctx - the context
351 - v - a vector capable of holding the interpolated field values
352
353 Level: intermediate
354
355 .seealso: [](ch_dmbase), `DM`, `DMInterpolationInfo`, `DMInterpolationGetVector()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`
356 @*/
DMInterpolationRestoreVector(DMInterpolationInfo ctx,Vec * v)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 }
365
DMInterpolate_Segment_Private(DMInterpolationInfo ctx,DM dm,PetscInt p,Vec xLocal,Vec v)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;
375
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 }
393
DMInterpolate_Triangle_Private(DMInterpolationInfo ctx,DM dm,PetscInt p,Vec xLocal,Vec v)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];
402
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 }
422
DMInterpolate_Tetrahedron_Private(DMInterpolationInfo ctx,DM dm,PetscInt p,Vec xLocal,Vec v)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];
432
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 }
452
QuadMap_Private(SNES snes,Vec Xref,Vec Xreal,PetscCtx ctx)453 static inline PetscErrorCode QuadMap_Private(SNES snes, Vec Xref, Vec Xreal, PetscCtx 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;
472
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];
479
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 }
488
489 #include <petsc/private/dmimpl.h>
QuadJacobian_Private(SNES snes,Vec Xref,Mat J,Mat M,PetscCtx ctx)490 static inline PetscErrorCode QuadJacobian_Private(SNES snes, Vec Xref, Mat J, Mat M, PetscCtx 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;
504
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];
512
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 }
525
DMInterpolate_Quad_Private(DMInterpolationInfo ctx,DM dm,PetscInt p,Vec xLocal,Vec v)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;
545
546 PetscFunctionBegin;
547 PetscCall(DMGetNumFields(dm, &Nf));
548 if (Nf) {
549 PetscObject obj;
550 PetscClassId id;
551
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));
578
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));
614
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 }
622
HexMap_Private(SNES snes,Vec Xref,Vec Xreal,PetscCtx ctx)623 static inline PetscErrorCode HexMap_Private(SNES snes, Vec Xref, Vec Xreal, PetscCtx 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;
673
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];
681
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 }
691
HexJacobian_Private(SNES snes,Vec Xref,Mat J,Mat M,PetscCtx ctx)692 static inline PetscErrorCode HexJacobian_Private(SNES snes, Vec Xref, Mat J, Mat M, PetscCtx 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;
732
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];
741
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;
751
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 }
760
DMInterpolate_Hex_Private(DMInterpolationInfo ctx,DM dm,PetscInt p,Vec xLocal,Vec v)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;
776
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));
797
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));
829
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 }
837
838 /*@
839 DMInterpolationEvaluate - Using the input from `dm` and `x`, calculates interpolated field values at the interpolation points.
840
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()`
845
846 Output Parameter:
847 . v - The vector containing the interpolated values, obtained with `DMInterpolationGetVector()`
848
849 Level: beginner
850
851 .seealso: [](ch_dmbase), `DM`, `DMInterpolationInfo`, `DMInterpolationGetVector()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`, `DMInterpolationGetCoordinates()`
852 @*/
DMInterpolationEvaluate(DMInterpolationInfo ctx,DM dm,Vec x,Vec v)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;
858
859 PetscFunctionBegin;
860 PetscValidHeaderSpecific(dm, DM_CLASSID, 2);
861 PetscValidHeaderSpecific(x, VEC_CLASSID, 3);
862 PetscValidHeaderSpecific(v, VEC_CLASSID, 4);
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;
873
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;
886
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;
894
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;
903
904 PetscCall(PetscDSGetDiscretization(ds, field, &obj));
905 PetscCall(PetscObjectGetClassId(obj, &id));
906 if (id == PETSCFE_CLASSID) {
907 PetscFE fe = (PetscFE)obj;
908
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;
914
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;
926
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;
944
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 }
969
970 /*@
971 DMInterpolationDestroy - Destroys a `DMInterpolationInfo` context
972
973 Collective
974
975 Input Parameter:
976 . ctx - the context
977
978 Level: beginner
979
980 .seealso: [](ch_dmbase), `DM`, `DMInterpolationInfo`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`
981 @*/
DMInterpolationDestroy(DMInterpolationInfo * ctx)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 }
993