utils/dm/dmadapt.c

#include <petscdmadaptor.h> /*I "petscdmadaptor.h" I*/
#include <petscdmplex.h>
#include <petscdmforest.h>
#include <petscds.h>
#include <petscblaslapack.h>
#include <petscsnes.h>

#include <petsc/private/dmadaptorimpl.h>
#include <petsc/private/dmpleximpl.h>
#include <petsc/private/petscfeimpl.h>

PetscClassId DMADAPTOR_CLASSID;

PetscFunctionList DMAdaptorList              = NULL;
PetscBool         DMAdaptorRegisterAllCalled = PETSC_FALSE;

/*@C
  DMAdaptorRegister - Adds a new adaptor component implementation

  Not Collective

  Input Parameters:
+ name        - The name of a new user-defined creation routine
- create_func - The creation routine

  Example Usage:
.vb
  DMAdaptorRegister("my_adaptor", MyAdaptorCreate);
.ve

  Then, your adaptor type can be chosen with the procedural interface via
.vb
  DMAdaptorCreate(MPI_Comm, DMAdaptor *);
  DMAdaptorSetType(DMAdaptor, "my_adaptor");
.ve
  or at runtime via the option
.vb
  -adaptor_type my_adaptor
.ve

  Level: advanced

  Note:
  `DMAdaptorRegister()` may be called multiple times to add several user-defined adaptors

.seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMAdaptor`, `DMAdaptorRegisterAll()`, `DMAdaptorRegisterDestroy()`
@*/
PetscErrorCode DMAdaptorRegister(const char name[], PetscErrorCode (*create_func)(DMAdaptor))
{
  PetscFunctionBegin;
  PetscCall(DMInitializePackage());
  PetscCall(PetscFunctionListAdd(&DMAdaptorList, name, create_func));
  PetscFunctionReturn(PETSC_SUCCESS);
}

PETSC_EXTERN PetscErrorCode DMAdaptorCreate_Gradient(DMAdaptor);
PETSC_EXTERN PetscErrorCode DMAdaptorCreate_Flux(DMAdaptor);

/*@C
  DMAdaptorRegisterAll - Registers all of the adaptor components in the `DM` package.

  Not Collective

  Level: advanced

.seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMAdaptorType`, `DMRegisterAll()`, `DMAdaptorRegisterDestroy()`
@*/
PetscErrorCode DMAdaptorRegisterAll(void)
{
  PetscFunctionBegin;
  if (DMAdaptorRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
  DMAdaptorRegisterAllCalled = PETSC_TRUE;

  PetscCall(DMAdaptorRegister(DMADAPTORGRADIENT, DMAdaptorCreate_Gradient));
  PetscCall(DMAdaptorRegister(DMADAPTORFLUX, DMAdaptorCreate_Flux));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  DMAdaptorRegisterDestroy - This function destroys the registered `DMAdaptorType`. It is called from `PetscFinalize()`.

  Not collective

  Level: developer

.seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMRegisterAll()`, `DMAdaptorType`, `PetscInitialize()`
@*/
PetscErrorCode DMAdaptorRegisterDestroy(void)
{
  PetscFunctionBegin;
  PetscCall(PetscFunctionListDestroy(&DMAdaptorList));
  DMAdaptorRegisterAllCalled = PETSC_FALSE;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  DMAdaptorCreate - Create a `DMAdaptor` object. Its purpose is to construct a adaptation `DMLabel` or metric `Vec` that can be used to modify the `DM`.

  Collective

  Input Parameter:
. comm - The communicator for the `DMAdaptor` object

  Output Parameter:
. adaptor - The `DMAdaptor` object

  Level: beginner

.seealso: [](ch_dmbase), `DM`, `DMAdaptor`, `DMAdaptorDestroy()`, `DMAdaptorAdapt()`, `PetscConvEst`, `PetscConvEstCreate()`
@*/
PetscErrorCode DMAdaptorCreate(MPI_Comm comm, DMAdaptor *adaptor)
{
  VecTaggerBox refineBox, coarsenBox;

  PetscFunctionBegin;
  PetscAssertPointer(adaptor, 2);
  PetscCall(PetscSysInitializePackage());

  PetscCall(PetscHeaderCreate(*adaptor, DMADAPTOR_CLASSID, "DMAdaptor", "DM Adaptor", "DMAdaptor", comm, DMAdaptorDestroy, DMAdaptorView));
  (*adaptor)->monitor          = PETSC_FALSE;
  (*adaptor)->adaptCriterion   = DM_ADAPTATION_NONE;
  (*adaptor)->numSeq           = 1;
  (*adaptor)->Nadapt           = -1;
  (*adaptor)->refinementFactor = 2.0;
  refineBox.min = refineBox.max = PETSC_MAX_REAL;
  PetscCall(VecTaggerCreate(PetscObjectComm((PetscObject)*adaptor), &(*adaptor)->refineTag));
  PetscCall(PetscObjectSetOptionsPrefix((PetscObject)(*adaptor)->refineTag, "refine_"));
  PetscCall(VecTaggerSetType((*adaptor)->refineTag, VECTAGGERABSOLUTE));
  PetscCall(VecTaggerAbsoluteSetBox((*adaptor)->refineTag, &refineBox));
  coarsenBox.min = coarsenBox.max = PETSC_MAX_REAL;
  PetscCall(VecTaggerCreate(PetscObjectComm((PetscObject)*adaptor), &(*adaptor)->coarsenTag));
  PetscCall(PetscObjectSetOptionsPrefix((PetscObject)(*adaptor)->coarsenTag, "coarsen_"));
  PetscCall(VecTaggerSetType((*adaptor)->coarsenTag, VECTAGGERABSOLUTE));
  PetscCall(VecTaggerAbsoluteSetBox((*adaptor)->coarsenTag, &coarsenBox));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  DMAdaptorDestroy - Destroys a `DMAdaptor` object

  Collective

  Input Parameter:
. adaptor - The `DMAdaptor` object

  Level: beginner

.seealso: [](ch_dmbase), `DM`, `DMAdaptor`, `DMAdaptorCreate()`, `DMAdaptorAdapt()`
@*/
PetscErrorCode DMAdaptorDestroy(DMAdaptor *adaptor)
{
  PetscFunctionBegin;
  if (!*adaptor) PetscFunctionReturn(PETSC_SUCCESS);
  PetscValidHeaderSpecific(*adaptor, DMADAPTOR_CLASSID, 1);
  if (--((PetscObject)*adaptor)->refct > 0) {
    *adaptor = NULL;
    PetscFunctionReturn(PETSC_SUCCESS);
  }
  PetscCall(VecTaggerDestroy(&(*adaptor)->refineTag));
  PetscCall(VecTaggerDestroy(&(*adaptor)->coarsenTag));
  PetscCall(PetscFree2((*adaptor)->exactSol, (*adaptor)->exactCtx));
  PetscCall(PetscHeaderDestroy(adaptor));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  DMAdaptorSetType - Sets the particular implementation for a adaptor.

  Collective

  Input Parameters:
+ adaptor - The `DMAdaptor`
- method  - The name of the adaptor type

  Options Database Key:
. -adaptor_type <type> - Sets the adaptor type; see `DMAdaptorType`

  Level: intermediate

.seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMAdaptor`, `DMAdaptorType`, `DMAdaptorGetType()`, `DMAdaptorCreate()`
@*/
PetscErrorCode DMAdaptorSetType(DMAdaptor adaptor, DMAdaptorType method)
{
  PetscErrorCode (*r)(DMAdaptor);
  PetscBool match;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(adaptor, DMADAPTOR_CLASSID, 1);
  PetscCall(PetscObjectTypeCompare((PetscObject)adaptor, method, &match));
  if (match) PetscFunctionReturn(PETSC_SUCCESS);

  PetscCall(DMAdaptorRegisterAll());
  PetscCall(PetscFunctionListFind(DMAdaptorList, method, &r));
  PetscCheck(r, PetscObjectComm((PetscObject)adaptor), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown DMAdaptor type: %s", method);

  PetscTryTypeMethod(adaptor, destroy);
  PetscCall(PetscMemzero(adaptor->ops, sizeof(*adaptor->ops)));
  PetscCall(PetscObjectChangeTypeName((PetscObject)adaptor, method));
  PetscCall((*r)(adaptor));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  DMAdaptorGetType - Gets the type name (as a string) from the adaptor.

  Not Collective

  Input Parameter:
. adaptor - The `DMAdaptor`

  Output Parameter:
. type - The `DMAdaptorType` name

  Level: intermediate

.seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMAdaptor`, `DMAdaptorType`, `DMAdaptorSetType()`, `DMAdaptorCreate()`
@*/
PetscErrorCode DMAdaptorGetType(DMAdaptor adaptor, DMAdaptorType *type)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(adaptor, DMADAPTOR_CLASSID, 1);
  PetscAssertPointer(type, 2);
  PetscCall(DMAdaptorRegisterAll());
  *type = ((PetscObject)adaptor)->type_name;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  DMAdaptorSetFromOptions - Sets properties of a `DMAdaptor` object from values in the options database

  Collective

  Input Parameter:
. adaptor - The `DMAdaptor` object

  Options Database Keys:
+ -adaptor_monitor <bool>              - Monitor the adaptation process
. -adaptor_sequence_num <num>          - Number of adaptations to generate an optimal grid
. -adaptor_target_num <num>            - Set the target number of vertices N_adapt, -1 for automatic determination
. -adaptor_refinement_factor <r>       - Set r such that N_adapt = r^dim N_orig
- -adaptor_mixed_setup_function <func> - Set the fnction func that sets up the mixed problem

  Level: beginner

.seealso: [](ch_dmbase), `DM`, `DMAdaptor`, `DMAdaptorCreate()`, `DMAdaptorAdapt()`
@*/
PetscErrorCode DMAdaptorSetFromOptions(DMAdaptor adaptor)
{
  char        typeName[PETSC_MAX_PATH_LEN];
  const char *defName = DMADAPTORGRADIENT;
  char        funcname[PETSC_MAX_PATH_LEN];
  PetscBool   flg;

  PetscFunctionBegin;
  PetscObjectOptionsBegin((PetscObject)adaptor);
  PetscCall(PetscOptionsFList("-adaptor_type", "DMAdaptor", "DMAdaptorSetType", DMAdaptorList, defName, typeName, 1024, &flg));
  if (flg) PetscCall(DMAdaptorSetType(adaptor, typeName));
  else if (!((PetscObject)adaptor)->type_name) PetscCall(DMAdaptorSetType(adaptor, defName));
  PetscCall(PetscOptionsBool("-adaptor_monitor", "Monitor the adaptation process", "DMAdaptorMonitor", adaptor->monitor, &adaptor->monitor, NULL));
  PetscCall(PetscOptionsInt("-adaptor_sequence_num", "Number of adaptations to generate an optimal grid", "DMAdaptorSetSequenceLength", adaptor->numSeq, &adaptor->numSeq, NULL));
  PetscCall(PetscOptionsInt("-adaptor_target_num", "Set the target number of vertices N_adapt, -1 for automatic determination", "DMAdaptor", adaptor->Nadapt, &adaptor->Nadapt, NULL));
  PetscCall(PetscOptionsReal("-adaptor_refinement_factor", "Set r such that N_adapt = r^dim N_orig", "DMAdaptor", adaptor->refinementFactor, &adaptor->refinementFactor, NULL));
  PetscCall(PetscOptionsString("-adaptor_mixed_setup_function", "Function to setup the mixed problem", "DMAdaptorSetMixedSetupFunction", funcname, funcname, sizeof(funcname), &flg));
  if (flg) {
    PetscErrorCode (*setupFunc)(DMAdaptor, DM);

    PetscCall(PetscDLSym(NULL, funcname, (void **)&setupFunc));
    PetscCheck(setupFunc, PetscObjectComm((PetscObject)adaptor), PETSC_ERR_ARG_WRONG, "Could not locate function %s", funcname);
    PetscCall(DMAdaptorSetMixedSetupFunction(adaptor, setupFunc));
  }
  PetscOptionsEnd();
  PetscCall(VecTaggerSetFromOptions(adaptor->refineTag));
  PetscCall(VecTaggerSetFromOptions(adaptor->coarsenTag));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  DMAdaptorView - Views a `DMAdaptor` object

  Collective

  Input Parameters:
+ adaptor - The `DMAdaptor` object
- viewer  - The `PetscViewer` object

  Level: beginner

.seealso: [](ch_dmbase), `DM`, `DMAdaptor`, `DMAdaptorCreate()`, `DMAdaptorAdapt()`
@*/
PetscErrorCode DMAdaptorView(DMAdaptor adaptor, PetscViewer viewer)
{
  PetscFunctionBegin;
  PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)adaptor, viewer));
  PetscCall(PetscViewerASCIIPrintf(viewer, "DM Adaptor\n"));
  PetscCall(PetscViewerASCIIPrintf(viewer, "  sequence length: %" PetscInt_FMT "\n", adaptor->numSeq));
  PetscCall(VecTaggerView(adaptor->refineTag, viewer));
  PetscCall(VecTaggerView(adaptor->coarsenTag, viewer));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  DMAdaptorGetSolver - Gets the solver used to produce discrete solutions

  Not Collective

  Input Parameter:
. adaptor - The `DMAdaptor` object

  Output Parameter:
. snes - The solver

  Level: intermediate

.seealso: [](ch_dmbase), `DM`, `DMAdaptor`, `DMAdaptorSetSolver()`, `DMAdaptorCreate()`, `DMAdaptorAdapt()`
@*/
PetscErrorCode DMAdaptorGetSolver(DMAdaptor adaptor, SNES *snes)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(adaptor, DMADAPTOR_CLASSID, 1);
  PetscAssertPointer(snes, 2);
  *snes = adaptor->snes;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  DMAdaptorSetSolver - Sets the solver used to produce discrete solutions

  Not Collective

  Input Parameters:
+ adaptor - The `DMAdaptor` object
- snes    - The solver, this MUST have an attached `DM`/`PetscDS`, so that the exact solution can be computed

  Level: intermediate

.seealso: [](ch_dmbase), `DMAdaptor`, `DMAdaptorGetSolver()`, `DMAdaptorCreate()`, `DMAdaptorAdapt()`
@*/
PetscErrorCode DMAdaptorSetSolver(DMAdaptor adaptor, SNES snes)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(adaptor, DMADAPTOR_CLASSID, 1);
  PetscValidHeaderSpecific(snes, SNES_CLASSID, 2);
  adaptor->snes = snes;
  PetscCall(SNESGetDM(adaptor->snes, &adaptor->idm));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  DMAdaptorGetSequenceLength - Gets the number of sequential adaptations used by an adapter

  Not Collective

  Input Parameter:
. adaptor - The `DMAdaptor` object

  Output Parameter:
. num - The number of adaptations

  Level: intermediate

.seealso: [](ch_dmbase), `DMAdaptor`, `DMAdaptorSetSequenceLength()`, `DMAdaptorCreate()`, `DMAdaptorAdapt()`
@*/
PetscErrorCode DMAdaptorGetSequenceLength(DMAdaptor adaptor, PetscInt *num)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(adaptor, DMADAPTOR_CLASSID, 1);
  PetscAssertPointer(num, 2);
  *num = adaptor->numSeq;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  DMAdaptorSetSequenceLength - Sets the number of sequential adaptations

  Not Collective

  Input Parameters:
+ adaptor - The `DMAdaptor` object
- num     - The number of adaptations

  Level: intermediate

.seealso: [](ch_dmbase), `DMAdaptorGetSequenceLength()`, `DMAdaptorCreate()`, `DMAdaptorAdapt()`
@*/
PetscErrorCode DMAdaptorSetSequenceLength(DMAdaptor adaptor, PetscInt num)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(adaptor, DMADAPTOR_CLASSID, 1);
  adaptor->numSeq = num;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode DMAdaptorTransferSolution_Exact_Private(DMAdaptor adaptor, DM dm, Vec u, DM adm, Vec au, void *ctx)
{
  PetscFunctionBeginUser;
  PetscCall(DMProjectFunction(adm, 0.0, adaptor->exactSol, adaptor->exactCtx, INSERT_ALL_VALUES, au));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  DMAdaptorSetUp - After the solver is specified, creates data structures for controlling adaptivity

  Collective

  Input Parameter:
. adaptor - The `DMAdaptor` object

  Level: beginner

.seealso: [](ch_dmbase), `DMAdaptor`, `DMAdaptorCreate()`, `DMAdaptorAdapt()`
@*/
PetscErrorCode DMAdaptorSetUp(DMAdaptor adaptor)
{
  PetscDS  prob;
  PetscInt Nf, f;

  PetscFunctionBegin;
  PetscCall(DMGetDS(adaptor->idm, &prob));
  PetscCall(VecTaggerSetUp(adaptor->refineTag));
  PetscCall(VecTaggerSetUp(adaptor->coarsenTag));
  PetscCall(PetscDSGetNumFields(prob, &Nf));
  PetscCall(PetscMalloc2(Nf, &adaptor->exactSol, Nf, &adaptor->exactCtx));
  for (f = 0; f < Nf; ++f) {
    PetscCall(PetscDSGetExactSolution(prob, f, &adaptor->exactSol[f], &adaptor->exactCtx[f]));
    /* TODO Have a flag that forces projection rather than using the exact solution */
    if (adaptor->exactSol[0]) PetscCall(DMAdaptorSetTransferFunction(adaptor, DMAdaptorTransferSolution_Exact_Private));
  }
  PetscTryTypeMethod(adaptor, setup);
  PetscFunctionReturn(PETSC_SUCCESS);
}

PetscErrorCode DMAdaptorGetTransferFunction(DMAdaptor adaptor, PetscErrorCode (**tfunc)(DMAdaptor, DM, Vec, DM, Vec, void *))
{
  PetscFunctionBegin;
  *tfunc = adaptor->ops->transfersolution;
  PetscFunctionReturn(PETSC_SUCCESS);
}

PetscErrorCode DMAdaptorSetTransferFunction(DMAdaptor adaptor, PetscErrorCode (*tfunc)(DMAdaptor, DM, Vec, DM, Vec, void *))
{
  PetscFunctionBegin;
  adaptor->ops->transfersolution = tfunc;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode DMAdaptorPreAdapt(DMAdaptor adaptor, Vec locX)
{
  DM           plex;
  PetscDS      prob;
  PetscObject  obj;
  PetscClassId id;
  PetscBool    isForest;

  PetscFunctionBegin;
  PetscCall(DMConvert(adaptor->idm, DMPLEX, &plex));
  PetscCall(DMGetDS(adaptor->idm, &prob));
  PetscCall(PetscDSGetDiscretization(prob, 0, &obj));
  PetscCall(PetscObjectGetClassId(obj, &id));
  PetscCall(DMIsForest(adaptor->idm, &isForest));
  if (adaptor->adaptCriterion == DM_ADAPTATION_NONE) {
    if (isForest) {
      adaptor->adaptCriterion = DM_ADAPTATION_LABEL;
    }
#if defined(PETSC_HAVE_PRAGMATIC)
    else {
      adaptor->adaptCriterion = DM_ADAPTATION_METRIC;
    }
#elif defined(PETSC_HAVE_MMG)
    else {
      adaptor->adaptCriterion = DM_ADAPTATION_METRIC;
    }
#elif defined(PETSC_HAVE_PARMMG)
    else {
      adaptor->adaptCriterion = DM_ADAPTATION_METRIC;
    }
#else
    else {
      adaptor->adaptCriterion = DM_ADAPTATION_LABEL;
    }
#endif
  }
  if (id == PETSCFV_CLASSID) {
    adaptor->femType = PETSC_FALSE;
  } else {
    adaptor->femType = PETSC_TRUE;
  }
  if (adaptor->femType) {
    /* Compute local solution bc */
    PetscCall(DMPlexInsertBoundaryValues(plex, PETSC_TRUE, locX, 0.0, adaptor->faceGeom, adaptor->cellGeom, NULL));
  } else {
    PetscFV      fvm = (PetscFV)obj;
    PetscLimiter noneLimiter;
    Vec          grad;

    PetscCall(PetscFVGetComputeGradients(fvm, &adaptor->computeGradient));
    PetscCall(PetscFVSetComputeGradients(fvm, PETSC_TRUE));
    /* Use no limiting when reconstructing gradients for adaptivity */
    PetscCall(PetscFVGetLimiter(fvm, &adaptor->limiter));
    PetscCall(PetscObjectReference((PetscObject)adaptor->limiter));
    PetscCall(PetscLimiterCreate(PetscObjectComm((PetscObject)fvm), &noneLimiter));
    PetscCall(PetscLimiterSetType(noneLimiter, PETSCLIMITERNONE));
    PetscCall(PetscFVSetLimiter(fvm, noneLimiter));
    /* Get FVM data */
    PetscCall(DMPlexGetDataFVM(plex, fvm, &adaptor->cellGeom, &adaptor->faceGeom, &adaptor->gradDM));
    PetscCall(VecGetDM(adaptor->cellGeom, &adaptor->cellDM));
    PetscCall(VecGetArrayRead(adaptor->cellGeom, &adaptor->cellGeomArray));
    /* Compute local solution bc */
    PetscCall(DMPlexInsertBoundaryValues(plex, PETSC_TRUE, locX, 0.0, adaptor->faceGeom, adaptor->cellGeom, NULL));
    /* Compute gradients */
    PetscCall(DMCreateGlobalVector(adaptor->gradDM, &grad));
    PetscCall(DMPlexReconstructGradientsFVM(plex, locX, grad));
    PetscCall(DMGetLocalVector(adaptor->gradDM, &adaptor->cellGrad));
    PetscCall(DMGlobalToLocalBegin(adaptor->gradDM, grad, INSERT_VALUES, adaptor->cellGrad));
    PetscCall(DMGlobalToLocalEnd(adaptor->gradDM, grad, INSERT_VALUES, adaptor->cellGrad));
    PetscCall(VecDestroy(&grad));
    PetscCall(VecGetArrayRead(adaptor->cellGrad, &adaptor->cellGradArray));
  }
  PetscCall(DMDestroy(&plex));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode DMAdaptorTransferSolution(DMAdaptor adaptor, DM dm, Vec x, DM adm, Vec ax)
{
  PetscReal time = 0.0;
  Mat       interp;
  void     *ctx;

  PetscFunctionBegin;
  PetscCall(DMGetApplicationContext(dm, &ctx));
  if (adaptor->ops->transfersolution) PetscUseTypeMethod(adaptor, transfersolution, dm, x, adm, ax, ctx);
  else {
    switch (adaptor->adaptCriterion) {
    case DM_ADAPTATION_LABEL:
      PetscCall(DMForestTransferVec(dm, x, adm, ax, PETSC_TRUE, time));
      break;
    case DM_ADAPTATION_REFINE:
    case DM_ADAPTATION_METRIC:
      PetscCall(DMCreateInterpolation(dm, adm, &interp, NULL));
      PetscCall(MatInterpolate(interp, x, ax));
      PetscCall(DMInterpolate(dm, interp, adm));
      PetscCall(MatDestroy(&interp));
      break;
    default:
      SETERRQ(PetscObjectComm((PetscObject)adaptor), PETSC_ERR_SUP, "No built-in projection for this adaptation criterion: %d", adaptor->adaptCriterion);
    }
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode DMAdaptorPostAdapt(DMAdaptor adaptor)
{
  PetscDS      prob;
  PetscObject  obj;
  PetscClassId id;

  PetscFunctionBegin;
  PetscCall(DMGetDS(adaptor->idm, &prob));
  PetscCall(PetscDSGetDiscretization(prob, 0, &obj));
  PetscCall(PetscObjectGetClassId(obj, &id));
  if (id == PETSCFV_CLASSID) {
    PetscFV fvm = (PetscFV)obj;

    PetscCall(PetscFVSetComputeGradients(fvm, adaptor->computeGradient));
    /* Restore original limiter */
    PetscCall(PetscFVSetLimiter(fvm, adaptor->limiter));

    PetscCall(VecRestoreArrayRead(adaptor->cellGeom, &adaptor->cellGeomArray));
    PetscCall(VecRestoreArrayRead(adaptor->cellGrad, &adaptor->cellGradArray));
    PetscCall(DMRestoreLocalVector(adaptor->gradDM, &adaptor->cellGrad));
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*
  DMAdaptorComputeCellErrorIndicator_Gradient - Use the integrated gradient as an error indicator in the `DMAdaptor`

  Input Parameters:
+ adaptor  - The `DMAdaptor` object
. dim      - The topological dimension
. cell     - The cell
. field    - The field integrated over the cell
. gradient - The gradient integrated over the cell
. cg       - A `PetscFVCellGeom` struct
- ctx      - A user context

  Output Parameter:
. errInd   - The error indicator

  Developer Note:
  Some of the input arguments are absurdly specialized to special situations, it is not clear this is a good general API

.seealso: [](ch_dmbase), `DMAdaptor`
*/
static PetscErrorCode DMAdaptorComputeCellErrorIndicator_Gradient(DMAdaptor adaptor, PetscInt dim, PetscInt Nc, const PetscScalar *field, const PetscScalar *gradient, const PetscFVCellGeom *cg, PetscReal *errInd, void *ctx)
{
  PetscReal err = 0.;
  PetscInt  c, d;

  PetscFunctionBeginHot;
  for (c = 0; c < Nc; c++) {
    for (d = 0; d < dim; ++d) err += PetscSqr(PetscRealPart(gradient[c * dim + d]));
  }
  *errInd = cg->volume * err;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode DMAdaptorComputeErrorIndicator_Gradient(DMAdaptor adaptor, Vec locX, Vec errVec)
{
  DM              dm, plex, edm, eplex;
  PetscDS         ds;
  PetscObject     obj;
  PetscClassId    id;
  void           *ctx;
  PetscQuadrature quad;
  PetscScalar    *earray;
  PetscReal       minMaxInd[2] = {PETSC_MAX_REAL, PETSC_MIN_REAL}, minMaxIndGlobal[2];
  PetscInt        dim, cdim, cStart, cEnd, Nf, Nc;

  PetscFunctionBegin;
  PetscCall(VecGetDM(locX, &dm));
  PetscCall(DMConvert(dm, DMPLEX, &plex));
  PetscCall(VecGetDM(errVec, &edm));
  PetscCall(DMConvert(edm, DMPLEX, &eplex));
  PetscCall(DMGetDimension(plex, &dim));
  PetscCall(DMGetCoordinateDim(plex, &cdim));
  PetscCall(DMGetApplicationContext(plex, &ctx));
  PetscCall(DMGetDS(plex, &ds));
  PetscCall(PetscDSGetNumFields(ds, &Nf));
  PetscCall(PetscDSGetDiscretization(ds, 0, &obj));
  PetscCall(PetscObjectGetClassId(obj, &id));

  PetscCall(VecGetArray(errVec, &earray));
  PetscCall(DMPlexGetSimplexOrBoxCells(plex, 0, &cStart, &cEnd));
  for (PetscInt cell = cStart; cell < cEnd; ++cell) {
    PetscScalar *eval;
    PetscReal    errInd = 0.;

    if (id == PETSCFV_CLASSID) {
      PetscFV            fv = (PetscFV)obj;
      const PetscScalar *pointSols;
      const PetscScalar *pointSol;
      const PetscScalar *pointGrad;
      PetscFVCellGeom   *cg;

      PetscCall(PetscFVGetNumComponents(fv, &Nc));
      PetscCall(VecGetArrayRead(locX, &pointSols));
      PetscCall(DMPlexPointLocalRead(plex, cell, pointSols, (void *)&pointSol));
      PetscCall(DMPlexPointLocalRead(adaptor->gradDM, cell, adaptor->cellGradArray, (void *)&pointGrad));
      PetscCall(DMPlexPointLocalRead(adaptor->cellDM, cell, adaptor->cellGeomArray, &cg));
      PetscUseTypeMethod(adaptor, computecellerrorindicator, dim, Nc, pointSol, pointGrad, cg, &errInd, ctx);
      PetscCall(VecRestoreArrayRead(locX, &pointSols));
    } else {
      PetscFE          fe = (PetscFE)obj;
      PetscScalar     *x  = NULL, *field, *gradient, *interpolant, *interpolantGrad;
      PetscFVCellGeom  cg;
      PetscFEGeom      fegeom;
      const PetscReal *quadWeights;
      PetscReal       *coords;
      PetscInt         Nb, Nq, qNc;

      fegeom.dim      = dim;
      fegeom.dimEmbed = cdim;
      PetscCall(PetscFEGetNumComponents(fe, &Nc));
      PetscCall(PetscFEGetQuadrature(fe, &quad));
      PetscCall(PetscFEGetDimension(fe, &Nb));
      PetscCall(PetscQuadratureGetData(quad, NULL, &qNc, &Nq, NULL, &quadWeights));
      PetscCall(PetscCalloc6(Nc, &field, cdim * Nc, &gradient, cdim * Nq, &coords, Nq, &fegeom.detJ, cdim * cdim * Nq, &fegeom.J, cdim * cdim * Nq, &fegeom.invJ));
      PetscCall(PetscMalloc2(Nc, &interpolant, cdim * Nc, &interpolantGrad));
      PetscCall(DMPlexComputeCellGeometryFEM(plex, cell, quad, coords, fegeom.J, fegeom.invJ, fegeom.detJ));
      PetscCall(DMPlexComputeCellGeometryFVM(plex, cell, &cg.volume, NULL, NULL));
      PetscCall(PetscArrayzero(gradient, cdim * Nc));
      PetscCall(DMPlexVecGetClosure(plex, NULL, locX, cell, NULL, &x));
      for (PetscInt f = 0; f < Nf; ++f) {
        PetscInt qc = 0;

        PetscCall(PetscDSGetDiscretization(ds, f, &obj));
        PetscCall(PetscArrayzero(interpolant, Nc));
        PetscCall(PetscArrayzero(interpolantGrad, cdim * Nc));
        for (PetscInt q = 0; q < Nq; ++q) {
          PetscCall(PetscFEInterpolateFieldAndGradient_Static((PetscFE)obj, 1, x, &fegeom, q, interpolant, interpolantGrad));
          for (PetscInt fc = 0; fc < Nc; ++fc) {
            const PetscReal wt = quadWeights[q * qNc + qc + fc];

            field[fc] += interpolant[fc] * wt * fegeom.detJ[q];
            for (PetscInt d = 0; d < cdim; ++d) gradient[fc * cdim + d] += interpolantGrad[fc * dim + d] * wt * fegeom.detJ[q];
          }
        }
        qc += Nc;
      }
      PetscCall(PetscFree2(interpolant, interpolantGrad));
      PetscCall(DMPlexVecRestoreClosure(plex, NULL, locX, cell, NULL, &x));
      for (PetscInt fc = 0; fc < Nc; ++fc) {
        field[fc] /= cg.volume;
        for (PetscInt d = 0; d < cdim; ++d) gradient[fc * cdim + d] /= cg.volume;
      }
      PetscUseTypeMethod(adaptor, computecellerrorindicator, dim, Nc, field, gradient, &cg, &errInd, ctx);
      PetscCall(PetscFree6(field, gradient, coords, fegeom.detJ, fegeom.J, fegeom.invJ));
    }
    PetscCall(DMPlexPointGlobalRef(eplex, cell, earray, (void *)&eval));
    eval[0]      = errInd;
    minMaxInd[0] = PetscMin(minMaxInd[0], errInd);
    minMaxInd[1] = PetscMax(minMaxInd[1], errInd);
  }
  PetscCall(VecRestoreArray(errVec, &earray));
  PetscCall(DMDestroy(&plex));
  PetscCall(DMDestroy(&eplex));
  PetscCall(PetscGlobalMinMaxReal(PetscObjectComm((PetscObject)adaptor), minMaxInd, minMaxIndGlobal));
  PetscCall(PetscInfo(adaptor, "DMAdaptor: error indicator range (%g, %g)\n", (double)minMaxIndGlobal[0], (double)minMaxIndGlobal[1]));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode DMAdaptorComputeErrorIndicator_Flux(DMAdaptor adaptor, Vec lu, Vec errVec)
{
  DM    dm, mdm;
  SNES  msnes;
  Vec   mu, lmu;
  void *ctx;

  PetscFunctionBegin;
  PetscCall(VecGetDM(lu, &dm));

  // Set up and solve mixed problem
  PetscCall(DMClone(dm, &mdm));
  PetscCall(SNESCreate(PetscObjectComm((PetscObject)mdm), &msnes));
  PetscCall(SNESSetDM(msnes, mdm));
  PetscCall(SNESAppendOptionsPrefix(msnes, "mixed_"));

  PetscTryTypeMethod(adaptor, mixedsetup, mdm);
  PetscCall(DMGetApplicationContext(dm, &ctx));
  PetscCall(DMPlexSetSNESLocalFEM(mdm, PETSC_FALSE, ctx));
  PetscCall(SNESSetFromOptions(msnes));

  PetscCall(DMCreateGlobalVector(mdm, &mu));
  PetscCall(PetscObjectSetName((PetscObject)mu, "Mixed Solution"));
  PetscCall(VecSet(mu, 0.0));
  PetscCall(SNESSolve(msnes, NULL, mu));
  PetscCall(VecViewFromOptions(mu, (PetscObject)adaptor, "-adapt_mixed_sol_vec_view"));

  PetscCall(DMGetLocalVector(mdm, &lmu));
  PetscCall(DMGlobalToLocal(mdm, mu, INSERT_VALUES, lmu));
  PetscCall(DMPlexInsertBoundaryValues(mdm, PETSC_TRUE, lmu, 0.0, NULL, NULL, NULL));
  PetscCall(DMPlexComputeL2FluxDiffVecLocal(lu, 0, lmu, 0, errVec));
  PetscCall(DMRestoreLocalVector(mdm, &lmu));
  PetscCall(VecDestroy(&mu));
  PetscCall(SNESDestroy(&msnes));
  PetscCall(DMDestroy(&mdm));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static void identityFunc(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f[])
{
  PetscInt i, j;

  for (i = 0; i < dim; ++i) {
    for (j = 0; j < dim; ++j) f[i + dim * j] = u[i + dim * j];
  }
}

static PetscErrorCode DMAdaptorAdapt_Sequence_Private(DMAdaptor adaptor, Vec inx, PetscBool doSolve, DM *adm, Vec *ax)
{
  PetscDS  ds;
  void    *ctx;
  MPI_Comm comm;
  PetscInt numAdapt = adaptor->numSeq, adaptIter;
  PetscInt dim, coordDim, Nf;

  PetscFunctionBegin;
  PetscCall(DMViewFromOptions(adaptor->idm, NULL, "-dm_adapt_pre_view"));
  PetscCall(VecViewFromOptions(inx, NULL, "-sol_adapt_pre_view"));
  PetscCall(PetscObjectGetComm((PetscObject)adaptor, &comm));
  PetscCall(DMGetDimension(adaptor->idm, &dim));
  PetscCall(DMGetCoordinateDim(adaptor->idm, &coordDim));
  PetscCall(DMGetApplicationContext(adaptor->idm, &ctx));
  PetscCall(DMGetDS(adaptor->idm, &ds));
  PetscCall(PetscDSGetNumFields(ds, &Nf));
  PetscCheck(Nf != 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot refine with no fields present!");

  /* Adapt until nothing changes */
  /* Adapt for a specified number of iterates */
  for (adaptIter = 0; adaptIter < numAdapt - 1; ++adaptIter) PetscCall(PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_(comm)));
  for (adaptIter = 0; adaptIter < numAdapt; ++adaptIter) {
    PetscBool adapted = PETSC_FALSE;
    DM        dm      = adaptIter ? *adm : adaptor->idm, odm;
    Vec       x       = adaptIter ? *ax : inx, locX, ox;

    PetscCall(DMGetLocalVector(dm, &locX));
    PetscCall(DMAdaptorPreAdapt(adaptor, locX));
    if (doSolve) {
      SNES snes;

      PetscCall(DMAdaptorGetSolver(adaptor, &snes));
      PetscCall(SNESSolve(snes, NULL, adaptIter ? *ax : x));
    }
    PetscCall(DMGlobalToLocalBegin(dm, adaptIter ? *ax : x, INSERT_VALUES, locX));
    PetscCall(DMGlobalToLocalEnd(dm, adaptIter ? *ax : x, INSERT_VALUES, locX));
    PetscCall(VecViewFromOptions(adaptIter ? *ax : x, (PetscObject)adaptor, "-adapt_primal_sol_vec_view"));
    /* PetscCall(DMAdaptorMonitor(adaptor));
       Print iterate, memory used, DM, solution */
    switch (adaptor->adaptCriterion) {
    case DM_ADAPTATION_REFINE:
      PetscCall(DMRefine(dm, comm, &odm));
      PetscCheck(odm, comm, PETSC_ERR_ARG_INCOMP, "DMRefine() did not perform any refinement, cannot continue grid sequencing");
      adapted = PETSC_TRUE;
      break;
    case DM_ADAPTATION_LABEL: {
      /* Adapt DM
           Create local solution
           Reconstruct gradients (FVM) or solve adjoint equation (FEM)
           Produce cellwise error indicator */
      DM             edm, plex;
      PetscFE        efe;
      DMLabel        adaptLabel;
      IS             refineIS, coarsenIS;
      Vec            errVec;
      DMPolytopeType ct;
      PetscInt       nRefine, nCoarsen, cStart;

      PetscCall(DMLabelCreate(PETSC_COMM_SELF, "adapt", &adaptLabel));

      // TODO Move this creation to PreAdapt
      PetscCall(DMClone(dm, &edm));
      PetscCall(DMConvert(edm, DMPLEX, &plex));
      PetscCall(DMPlexGetHeightStratum(plex, 0, &cStart, NULL));
      PetscCall(DMPlexGetCellType(plex, cStart, &ct));
      PetscCall(DMDestroy(&plex));
      PetscCall(PetscFECreateLagrangeByCell(PETSC_COMM_SELF, dim, 1, ct, 0, PETSC_DEFAULT, &efe));
      PetscCall(PetscObjectSetName((PetscObject)efe, "Error"));
      PetscCall(DMSetField(edm, 0, NULL, (PetscObject)efe));
      PetscCall(PetscFEDestroy(&efe));
      PetscCall(DMCreateDS(edm));
      PetscCall(DMGetGlobalVector(edm, &errVec));
      PetscCall(PetscObjectSetName((PetscObject)errVec, "Error Estimator"));

      PetscUseTypeMethod(adaptor, computeerrorindicator, locX, errVec);
      PetscCall(VecViewFromOptions(errVec, (PetscObject)adaptor, "-adapt_error_vec_view"));

      // Compute IS from VecTagger
      PetscCall(VecTaggerComputeIS(adaptor->refineTag, errVec, &refineIS, NULL));
      PetscCall(VecTaggerComputeIS(adaptor->coarsenTag, errVec, &coarsenIS, NULL));
      PetscCall(ISViewFromOptions(refineIS, (PetscObject)adaptor->refineTag, "-is_view"));
      PetscCall(ISViewFromOptions(coarsenIS, (PetscObject)adaptor->coarsenTag, "-is_view"));
      PetscCall(ISGetSize(refineIS, &nRefine));
      PetscCall(ISGetSize(coarsenIS, &nCoarsen));
      PetscCall(PetscInfo(adaptor, "DMAdaptor: numRefine %" PetscInt_FMT ", numCoarsen %" PetscInt_FMT "\n", nRefine, nCoarsen));
      if (nRefine) PetscCall(DMLabelSetStratumIS(adaptLabel, DM_ADAPT_REFINE, refineIS));
      if (nCoarsen) PetscCall(DMLabelSetStratumIS(adaptLabel, DM_ADAPT_COARSEN, coarsenIS));
      PetscCall(ISDestroy(&coarsenIS));
      PetscCall(ISDestroy(&refineIS));
      PetscCall(DMRestoreGlobalVector(edm, &errVec));
      PetscCall(DMDestroy(&edm));
      // Adapt DM from label
      if (nRefine || nCoarsen) {
        char        oprefix[PETSC_MAX_PATH_LEN];
        const char *p;
        PetscBool   flg;

        PetscCall(PetscOptionsHasName(NULL, adaptor->hdr.prefix, "-adapt_vec_view", &flg));
        if (flg) {
          Vec ref;

          PetscCall(DMPlexCreateLabelField(dm, adaptLabel, &ref));
          PetscCall(VecViewFromOptions(ref, (PetscObject)adaptor, "-adapt_vec_view"));
          PetscCall(VecDestroy(&ref));
        }

        PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &p));
        PetscCall(PetscStrncpy(oprefix, p, PETSC_MAX_PATH_LEN));
        PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)dm, "adapt_"));
        PetscCall(DMAdaptLabel(dm, adaptLabel, &odm));
        PetscCall(PetscObjectSetOptionsPrefix((PetscObject)dm, oprefix));
        PetscCall(PetscObjectSetOptionsPrefix((PetscObject)odm, oprefix));
        adapted = PETSC_TRUE;
      }
      PetscCall(DMLabelDestroy(&adaptLabel));
    } break;
    case DM_ADAPTATION_METRIC: {
      DM        dmGrad, dmHess, dmMetric, dmDet;
      Vec       xGrad, xHess, metric, determinant;
      PetscReal N;
      DMLabel   bdLabel = NULL, rgLabel = NULL;
      PetscBool higherOrder = PETSC_FALSE;
      PetscInt  Nd          = coordDim * coordDim, f, vStart, vEnd;
      void (**funcs)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]);

      PetscCall(PetscMalloc(1, &funcs));
      funcs[0] = identityFunc;

      /*     Setup finite element spaces */
      PetscCall(DMClone(dm, &dmGrad));
      PetscCall(DMClone(dm, &dmHess));
      PetscCheck(Nf <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Adaptation with multiple fields not yet considered"); // TODO
      for (f = 0; f < Nf; ++f) {
        PetscFE         fe, feGrad, feHess;
        PetscDualSpace  Q;
        PetscSpace      space;
        DM              K;
        PetscQuadrature q;
        PetscInt        Nc, qorder, p;
        const char     *prefix;

        PetscCall(PetscDSGetDiscretization(ds, f, (PetscObject *)&fe));
        PetscCall(PetscFEGetNumComponents(fe, &Nc));
        PetscCheck(Nc <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Adaptation with multiple components not yet considered"); // TODO
        PetscCall(PetscFEGetBasisSpace(fe, &space));
        PetscCall(PetscSpaceGetDegree(space, NULL, &p));
        if (p > 1) higherOrder = PETSC_TRUE;
        PetscCall(PetscFEGetDualSpace(fe, &Q));
        PetscCall(PetscDualSpaceGetDM(Q, &K));
        PetscCall(DMPlexGetDepthStratum(K, 0, &vStart, &vEnd));
        PetscCall(PetscFEGetQuadrature(fe, &q));
        PetscCall(PetscQuadratureGetOrder(q, &qorder));
        PetscCall(PetscObjectGetOptionsPrefix((PetscObject)fe, &prefix));
        PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)dmGrad), dim, Nc * coordDim, PETSC_TRUE, prefix, qorder, &feGrad));
        PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)dmHess), dim, Nc * Nd, PETSC_TRUE, prefix, qorder, &feHess));
        PetscCall(DMSetField(dmGrad, f, NULL, (PetscObject)feGrad));
        PetscCall(DMSetField(dmHess, f, NULL, (PetscObject)feHess));
        PetscCall(DMCreateDS(dmGrad));
        PetscCall(DMCreateDS(dmHess));
        PetscCall(PetscFEDestroy(&feGrad));
        PetscCall(PetscFEDestroy(&feHess));
      }
      /*     Compute vertexwise gradients from cellwise gradients */
      PetscCall(DMCreateLocalVector(dmGrad, &xGrad));
      PetscCall(VecViewFromOptions(locX, NULL, "-sol_adapt_loc_pre_view"));
      PetscCall(DMPlexComputeGradientClementInterpolant(dm, locX, xGrad));
      PetscCall(VecViewFromOptions(xGrad, NULL, "-adapt_gradient_view"));
      /*     Compute vertexwise Hessians from cellwise Hessians */
      PetscCall(DMCreateLocalVector(dmHess, &xHess));
      PetscCall(DMPlexComputeGradientClementInterpolant(dmGrad, xGrad, xHess));
      PetscCall(VecViewFromOptions(xHess, NULL, "-adapt_hessian_view"));
      PetscCall(VecDestroy(&xGrad));
      PetscCall(DMDestroy(&dmGrad));
      /*     Compute L-p normalized metric */
      PetscCall(DMClone(dm, &dmMetric));
      N = adaptor->Nadapt >= 0 ? adaptor->Nadapt : PetscPowRealInt(adaptor->refinementFactor, dim) * ((PetscReal)(vEnd - vStart));
      if (adaptor->monitor) {
        PetscMPIInt rank, size;
        PetscCallMPI(MPI_Comm_rank(comm, &size));
        PetscCallMPI(MPI_Comm_rank(comm, &rank));
        PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] N_orig: %" PetscInt_FMT " N_adapt: %g\n", rank, vEnd - vStart, (double)N));
      }
      PetscCall(DMPlexMetricSetTargetComplexity(dmMetric, (PetscReal)N));
      if (higherOrder) {
        /*   Project Hessian into P1 space, if required */
        PetscCall(DMPlexMetricCreate(dmMetric, 0, &metric));
        PetscCall(DMProjectFieldLocal(dmMetric, 0.0, xHess, funcs, INSERT_ALL_VALUES, metric));
        PetscCall(VecDestroy(&xHess));
        xHess = metric;
      }
      PetscCall(PetscFree(funcs));
      PetscCall(DMPlexMetricCreate(dmMetric, 0, &metric));
      PetscCall(DMPlexMetricDeterminantCreate(dmMetric, 0, &determinant, &dmDet));
      PetscCall(DMPlexMetricNormalize(dmMetric, xHess, PETSC_TRUE, PETSC_TRUE, metric, determinant));
      PetscCall(VecDestroy(&determinant));
      PetscCall(DMDestroy(&dmDet));
      PetscCall(VecDestroy(&xHess));
      PetscCall(DMDestroy(&dmHess));
      /*     Adapt DM from metric */
      PetscCall(DMGetLabel(dm, "marker", &bdLabel));
      PetscCall(DMAdaptMetric(dm, metric, bdLabel, rgLabel, &odm));
      adapted = PETSC_TRUE;
      /*     Cleanup */
      PetscCall(VecDestroy(&metric));
      PetscCall(DMDestroy(&dmMetric));
    } break;
    default:
      SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid adaptation type: %d", adaptor->adaptCriterion);
    }
    PetscCall(DMAdaptorPostAdapt(adaptor));
    PetscCall(DMRestoreLocalVector(dm, &locX));
    /* If DM was adapted, replace objects and recreate solution */
    if (adapted) {
      const char *name;

      PetscCall(PetscObjectGetName((PetscObject)dm, &name));
      PetscCall(PetscObjectSetName((PetscObject)odm, name));
      /* Reconfigure solver */
      PetscCall(SNESReset(adaptor->snes));
      PetscCall(SNESSetDM(adaptor->snes, odm));
      PetscCall(DMAdaptorSetSolver(adaptor, adaptor->snes));
      PetscCall(DMPlexSetSNESLocalFEM(odm, PETSC_FALSE, ctx));
      PetscCall(SNESSetFromOptions(adaptor->snes));
      /* Transfer system */
      PetscCall(DMCopyDisc(dm, odm));
      /* Transfer solution */
      PetscCall(DMCreateGlobalVector(odm, &ox));
      PetscCall(PetscObjectGetName((PetscObject)x, &name));
      PetscCall(PetscObjectSetName((PetscObject)ox, name));
      PetscCall(DMAdaptorTransferSolution(adaptor, dm, x, odm, ox));
      /* Cleanup adaptivity info */
      if (adaptIter > 0) PetscCall(PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_(comm)));
      PetscCall(DMForestSetAdaptivityForest(dm, NULL)); /* clear internal references to the previous dm */
      PetscCall(DMDestroy(&dm));
      PetscCall(VecDestroy(&x));
      *adm = odm;
      *ax  = ox;
    } else {
      *adm      = dm;
      *ax       = x;
      adaptIter = numAdapt;
    }
    if (adaptIter < numAdapt - 1) {
      PetscCall(DMViewFromOptions(odm, NULL, "-dm_adapt_iter_view"));
      PetscCall(VecViewFromOptions(ox, NULL, "-sol_adapt_iter_view"));
    }
  }
  PetscCall(DMViewFromOptions(*adm, NULL, "-dm_adapt_view"));
  PetscCall(VecViewFromOptions(*ax, NULL, "-sol_adapt_view"));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  DMAdaptorAdapt - Creates a new `DM` that is adapted to the problem

  Not Collective

  Input Parameters:
+ adaptor  - The `DMAdaptor` object
. x        - The global approximate solution
- strategy - The adaptation strategy, see `DMAdaptationStrategy`

  Output Parameters:
+ adm - The adapted `DM`
- ax  - The adapted solution

  Options Database Keys:
+ -snes_adapt <strategy> - initial, sequential, multigrid
. -adapt_gradient_view   - View the Clement interpolant of the solution gradient
. -adapt_hessian_view    - View the Clement interpolant of the solution Hessian
- -adapt_metric_view     - View the metric tensor for adaptive mesh refinement

  Level: intermediate

.seealso: [](ch_dmbase), `DMAdaptor`, `DMAdaptationStrategy`, `DMAdaptorSetSolver()`, `DMAdaptorCreate()`
@*/
PetscErrorCode DMAdaptorAdapt(DMAdaptor adaptor, Vec x, DMAdaptationStrategy strategy, DM *adm, Vec *ax)
{
  PetscFunctionBegin;
  switch (strategy) {
  case DM_ADAPTATION_INITIAL:
    PetscCall(DMAdaptorAdapt_Sequence_Private(adaptor, x, PETSC_FALSE, adm, ax));
    break;
  case DM_ADAPTATION_SEQUENTIAL:
    PetscCall(DMAdaptorAdapt_Sequence_Private(adaptor, x, PETSC_TRUE, adm, ax));
    break;
  default:
    SETERRQ(PetscObjectComm((PetscObject)adaptor), PETSC_ERR_ARG_WRONG, "Unrecognized adaptation strategy %d", strategy);
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  DMAdaptorGetMixedSetupFunction - Get the function setting up the mixed problem, if it exists

  Not Collective

  Input Parameter:
. adaptor - the `DMAdaptor`

  Output Parameter:
. setupFunc - the function setting up the mixed problem, or `NULL`

  Level: advanced

.seealso: `DMAdaptor`, `DMAdaptorSetMixedSetupFunction()`, `DMAdaptorAdapt()`
@*/
PetscErrorCode DMAdaptorGetMixedSetupFunction(DMAdaptor adaptor, PetscErrorCode (**setupFunc)(DMAdaptor, DM))
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(adaptor, DMADAPTOR_CLASSID, 1);
  PetscAssertPointer(setupFunc, 2);
  *setupFunc = adaptor->ops->mixedsetup;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  DMAdaptorSetMixedSetupFunction - Set the function setting up the mixed problem

  Not Collective

  Input Parameters:
+ adaptor   - the `DMAdaptor`
- setupFunc - the function setting up the mixed problem

  Level: advanced

.seealso: `DMAdaptor`, `DMAdaptorGetMixedSetupFunction()`, `DMAdaptorAdapt()`
@*/
PetscErrorCode DMAdaptorSetMixedSetupFunction(DMAdaptor adaptor, PetscErrorCode (*setupFunc)(DMAdaptor, DM))
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(adaptor, DMADAPTOR_CLASSID, 1);
  PetscValidFunction(setupFunc, 2);
  adaptor->ops->mixedsetup = setupFunc;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode DMAdaptorInitialize_Gradient(DMAdaptor adaptor)
{
  PetscFunctionBegin;
  adaptor->ops->computeerrorindicator     = DMAdaptorComputeErrorIndicator_Gradient;
  adaptor->ops->computecellerrorindicator = DMAdaptorComputeCellErrorIndicator_Gradient;
  PetscFunctionReturn(PETSC_SUCCESS);
}

PETSC_EXTERN PetscErrorCode DMAdaptorCreate_Gradient(DMAdaptor adaptor)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(adaptor, DMADAPTOR_CLASSID, 1);
  adaptor->data = NULL;

  PetscCall(DMAdaptorInitialize_Gradient(adaptor));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode DMAdaptorInitialize_Flux(DMAdaptor adaptor)
{
  PetscFunctionBegin;
  adaptor->ops->computeerrorindicator = DMAdaptorComputeErrorIndicator_Flux;
  PetscFunctionReturn(PETSC_SUCCESS);
}

PETSC_EXTERN PetscErrorCode DMAdaptorCreate_Flux(DMAdaptor adaptor)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(adaptor, DMADAPTOR_CLASSID, 1);
  adaptor->data = NULL;

  PetscCall(DMAdaptorInitialize_Flux(adaptor));
  PetscFunctionReturn(PETSC_SUCCESS);
}