Actual source code: dminterpolatesnes.c

      1: #include <petsc/private/dmpleximpl.h>
      2: #include <petsc/private/snesimpl.h>
      3: #include <petscds.h>
      4: #include <petsc/private/petscimpl.h>
      5: #include <petsc/private/petscfeimpl.h>

      7: /*@
      8:   DMInterpolationCreate - Creates a `DMInterpolationInfo` context

     10:   Collective

     12:   Input Parameter:
     13: . comm - the communicator

     15:   Output Parameter:
     16: . ctx - the context

     18:   Level: beginner

     20:   Developer Note:
     21:   The naming is incorrect, either the object should be named `DMInterpolation` or all the routines should begin with `DMInterpolationInfo`

     23: .seealso: [](ch_dmbase), `DM`, `DMInterpolationInfo`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationDestroy()`
     24: @*/
     25: PetscErrorCode DMInterpolationCreate(MPI_Comm comm, DMInterpolationInfo *ctx)
     26: {
     27:   PetscFunctionBegin;
     28:   PetscAssertPointer(ctx, 2);
     29:   PetscCall(PetscNew(ctx));

     31:   (*ctx)->comm   = comm;
     32:   (*ctx)->dim    = -1;
     33:   (*ctx)->nInput = 0;
     34:   (*ctx)->points = NULL;
     35:   (*ctx)->cells  = NULL;
     36:   (*ctx)->n      = -1;
     37:   (*ctx)->coords = NULL;
     38:   PetscFunctionReturn(PETSC_SUCCESS);
     39: }

     41: /*@
     42:   DMInterpolationSetDim - Sets the spatial dimension for the interpolation context

     44:   Not Collective

     46:   Input Parameters:
     47: + ctx - the context
     48: - dim - the spatial dimension

     50:   Level: intermediate

     52: .seealso: [](ch_dmbase), `DM`, `DMInterpolationInfo`, `DMInterpolationGetDim()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`
     53: @*/
     54: PetscErrorCode DMInterpolationSetDim(DMInterpolationInfo ctx, PetscInt dim)
     55: {
     56:   PetscFunctionBegin;
     57:   PetscCheck(!(dim < 1) && !(dim > 3), ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension for points: %" PetscInt_FMT, dim);
     58:   ctx->dim = dim;
     59:   PetscFunctionReturn(PETSC_SUCCESS);
     60: }

     62: /*@
     63:   DMInterpolationGetDim - Gets the spatial dimension for the interpolation context

     65:   Not Collective

     67:   Input Parameter:
     68: . ctx - the context

     70:   Output Parameter:
     71: . dim - the spatial dimension

     73:   Level: intermediate

     75: .seealso: [](ch_dmbase), `DM`, `DMInterpolationInfo`, `DMInterpolationSetDim()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`
     76: @*/
     77: PetscErrorCode DMInterpolationGetDim(DMInterpolationInfo ctx, PetscInt *dim)
     78: {
     79:   PetscFunctionBegin;
     80:   PetscAssertPointer(dim, 2);
     81:   *dim = ctx->dim;
     82:   PetscFunctionReturn(PETSC_SUCCESS);
     83: }

     85: /*@
     86:   DMInterpolationSetDof - Sets the number of fields interpolated at a point for the interpolation context

     88:   Not Collective

     90:   Input Parameters:
     91: + ctx - the context
     92: - dof - the number of fields

     94:   Level: intermediate

     96: .seealso: [](ch_dmbase), `DM`, `DMInterpolationInfo`, `DMInterpolationGetDof()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`
     97: @*/
     98: PetscErrorCode DMInterpolationSetDof(DMInterpolationInfo ctx, PetscInt dof)
     99: {
    100:   PetscFunctionBegin;
    101:   PetscCheck(dof >= 1, ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid number of components: %" PetscInt_FMT, dof);
    102:   ctx->dof = dof;
    103:   PetscFunctionReturn(PETSC_SUCCESS);
    104: }

    106: /*@
    107:   DMInterpolationGetDof - Gets the number of fields interpolated at a point for the interpolation context

    109:   Not Collective

    111:   Input Parameter:
    112: . ctx - the context

    114:   Output Parameter:
    115: . dof - the number of fields

    117:   Level: intermediate

    119: .seealso: [](ch_dmbase), `DM`, `DMInterpolationInfo`, `DMInterpolationSetDof()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`
    120: @*/
    121: PetscErrorCode DMInterpolationGetDof(DMInterpolationInfo ctx, PetscInt *dof)
    122: {
    123:   PetscFunctionBegin;
    124:   PetscAssertPointer(dof, 2);
    125:   *dof = ctx->dof;
    126:   PetscFunctionReturn(PETSC_SUCCESS);
    127: }

    129: /*@
    130:   DMInterpolationAddPoints - Add points at which we will interpolate the fields

    132:   Not Collective

    134:   Input Parameters:
    135: + ctx    - the context
    136: . n      - the number of points
    137: - points - the coordinates for each point, an array of size `n` * dim

    139:   Level: intermediate

    141:   Note:
    142:   The input coordinate information is copied into the object.

    144: .seealso: [](ch_dmbase), `DM`, `DMInterpolationInfo`, `DMInterpolationSetDim()`, `DMInterpolationEvaluate()`, `DMInterpolationCreate()`
    145: @*/
    146: PetscErrorCode DMInterpolationAddPoints(DMInterpolationInfo ctx, PetscInt n, PetscReal points[])
    147: {
    148:   PetscFunctionBegin;
    149:   PetscCheck(ctx->dim >= 0, ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The spatial dimension has not been set");
    150:   PetscCheck(!ctx->points, ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "Cannot add points multiple times yet");
    151:   ctx->nInput = n;

    153:   PetscCall(PetscMalloc1(n * ctx->dim, &ctx->points));
    154:   PetscCall(PetscArraycpy(ctx->points, points, n * ctx->dim));
    155:   PetscFunctionReturn(PETSC_SUCCESS);
    156: }

    158: /*@
    159:   DMInterpolationSetUp - Compute spatial indices for point location during interpolation

    161:   Collective

    163:   Input Parameters:
    164: + ctx                 - the context
    165: . dm                  - the `DM` for the function space used for interpolation
    166: . redundantPoints     - If `PETSC_TRUE`, all processes are passing in the same array of points. Otherwise, points need to be communicated among processes.
    167: - ignoreOutsideDomain - If `PETSC_TRUE`, ignore points outside the domain, otherwise return an error

    169:   Level: intermediate

    171: .seealso: [](ch_dmbase), `DM`, `DMInterpolationInfo`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`
    172: @*/
    173: PetscErrorCode DMInterpolationSetUp(DMInterpolationInfo ctx, DM dm, PetscBool redundantPoints, PetscBool ignoreOutsideDomain)
    174: {
    175:   MPI_Comm           comm = ctx->comm;
    176:   PetscScalar       *a;
    177:   PetscInt           p, q, i;
    178:   PetscMPIInt        rank, size;
    179:   Vec                pointVec;
    180:   PetscSF            cellSF;
    181:   PetscLayout        layout;
    182:   PetscReal         *globalPoints;
    183:   PetscScalar       *globalPointsScalar;
    184:   const PetscInt    *ranges;
    185:   PetscMPIInt       *counts, *displs;
    186:   const PetscSFNode *foundCells;
    187:   const PetscInt    *foundPoints;
    188:   PetscMPIInt       *foundProcs, *globalProcs, in;
    189:   PetscInt           n, N, numFound;

    191:   PetscFunctionBegin;
    193:   PetscCallMPI(MPI_Comm_size(comm, &size));
    194:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
    195:   PetscCheck(ctx->dim >= 0, comm, PETSC_ERR_ARG_WRONGSTATE, "The spatial dimension has not been set");
    196:   /* Locate points */
    197:   n = ctx->nInput;
    198:   if (!redundantPoints) {
    199:     PetscCall(PetscLayoutCreate(comm, &layout));
    200:     PetscCall(PetscLayoutSetBlockSize(layout, 1));
    201:     PetscCall(PetscLayoutSetLocalSize(layout, n));
    202:     PetscCall(PetscLayoutSetUp(layout));
    203:     PetscCall(PetscLayoutGetSize(layout, &N));
    204:     /* Communicate all points to all processes */
    205:     PetscCall(PetscMalloc3(N * ctx->dim, &globalPoints, size, &counts, size, &displs));
    206:     PetscCall(PetscLayoutGetRanges(layout, &ranges));
    207:     for (p = 0; p < size; ++p) {
    208:       PetscCall(PetscMPIIntCast((ranges[p + 1] - ranges[p]) * ctx->dim, &counts[p]));
    209:       PetscCall(PetscMPIIntCast(ranges[p] * ctx->dim, &displs[p]));
    210:     }
    211:     PetscCall(PetscMPIIntCast(n * ctx->dim, &in));
    212:     PetscCallMPI(MPI_Allgatherv(ctx->points, in, MPIU_REAL, globalPoints, counts, displs, MPIU_REAL, comm));
    213:   } else {
    214:     N            = n;
    215:     globalPoints = ctx->points;
    216:     counts = displs = NULL;
    217:     layout          = NULL;
    218:   }
    219: #if 0
    220:   PetscCall(PetscMalloc3(N,&foundCells,N,&foundProcs,N,&globalProcs));
    221:   /* foundCells[p] = m->locatePoint(&globalPoints[p*ctx->dim]); */
    222: #else
    223:   #if defined(PETSC_USE_COMPLEX)
    224:   PetscCall(PetscMalloc1(N * ctx->dim, &globalPointsScalar));
    225:   for (i = 0; i < N * ctx->dim; i++) globalPointsScalar[i] = globalPoints[i];
    226:   #else
    227:   globalPointsScalar = globalPoints;
    228:   #endif
    229:   PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, ctx->dim, N * ctx->dim, globalPointsScalar, &pointVec));
    230:   PetscCall(PetscMalloc2(N, &foundProcs, N, &globalProcs));
    231:   for (p = 0; p < N; ++p) foundProcs[p] = size;
    232:   cellSF = NULL;
    233:   PetscCall(DMLocatePoints(dm, pointVec, DM_POINTLOCATION_REMOVE, &cellSF));
    234:   PetscCall(PetscSFGetGraph(cellSF, NULL, &numFound, &foundPoints, &foundCells));
    235: #endif
    236:   for (p = 0; p < numFound; ++p) {
    237:     if (foundCells[p].index >= 0) foundProcs[foundPoints ? foundPoints[p] : p] = rank;
    238:   }
    239:   /* Let the lowest rank process own each point */
    240:   PetscCallMPI(MPIU_Allreduce(foundProcs, globalProcs, N, MPI_INT, MPI_MIN, comm));
    241:   ctx->n = 0;
    242:   for (p = 0; p < N; ++p) {
    243:     if (globalProcs[p] == size) {
    244:       PetscCheck(ignoreOutsideDomain, comm, PETSC_ERR_PLIB, "Point %" PetscInt_FMT ": %g %g %g not located in mesh", p, (double)globalPoints[p * ctx->dim + 0], (double)(ctx->dim > 1 ? globalPoints[p * ctx->dim + 1] : 0),
    245:                  (double)(ctx->dim > 2 ? globalPoints[p * ctx->dim + 2] : 0));
    246:       if (rank == 0) ++ctx->n;
    247:     } else if (globalProcs[p] == rank) ++ctx->n;
    248:   }
    249:   /* Create coordinates vector and array of owned cells */
    250:   PetscCall(PetscMalloc1(ctx->n, &ctx->cells));
    251:   PetscCall(VecCreate(comm, &ctx->coords));
    252:   PetscCall(VecSetSizes(ctx->coords, ctx->n * ctx->dim, PETSC_DECIDE));
    253:   PetscCall(VecSetBlockSize(ctx->coords, ctx->dim));
    254:   PetscCall(VecSetType(ctx->coords, VECSTANDARD));
    255:   PetscCall(VecGetArray(ctx->coords, &a));
    256:   for (p = 0, q = 0, i = 0; p < N; ++p) {
    257:     if (globalProcs[p] == rank) {
    258:       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: /*@
    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: /*@
    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: /*@
    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: /*@
    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: /*@
    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: }