#include <petsc/private/petscdsimpl.h> /*I "petscds.h" I*/

PetscClassId PETSCDS_CLASSID = 0;

PetscFunctionList PetscDSList              = NULL;
PetscBool         PetscDSRegisterAllCalled = PETSC_FALSE;

/*@C
  PetscDSRegister - Adds a new `PetscDS` implementation

  Not Collective; No Fortran Support

  Input Parameters:
+ sname    - The name of a new user-defined creation routine
- function - The creation routine itself

  Example Usage:
.vb
    PetscDSRegister("my_ds", MyPetscDSCreate);
.ve

  Then, your PetscDS type can be chosen with the procedural interface via
.vb
    PetscDSCreate(MPI_Comm, PetscDS *);
    PetscDSSetType(PetscDS, "my_ds");
.ve
  or at runtime via the option
.vb
    -petscds_type my_ds
.ve

  Level: advanced

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

.seealso: `PetscDSType`, `PetscDS`, `PetscDSRegisterAll()`, `PetscDSRegisterDestroy()`
@*/
PetscErrorCode PetscDSRegister(const char sname[], PetscErrorCode (*function)(PetscDS))
{
  PetscFunctionBegin;
  PetscCall(PetscFunctionListAdd(&PetscDSList, sname, function));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSSetType - Builds a particular `PetscDS`

  Collective; No Fortran Support

  Input Parameters:
+ prob - The `PetscDS` object
- name - The `PetscDSType`

  Options Database Key:
. -petscds_type <type> - Sets the PetscDS type; use -help for a list of available types

  Level: intermediate

.seealso: `PetscDSType`, `PetscDS`, `PetscDSGetType()`, `PetscDSCreate()`
@*/
PetscErrorCode PetscDSSetType(PetscDS prob, PetscDSType name)
{
  PetscErrorCode (*r)(PetscDS);
  PetscBool match;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
  PetscCall(PetscObjectTypeCompare((PetscObject)prob, name, &match));
  if (match) PetscFunctionReturn(PETSC_SUCCESS);

  PetscCall(PetscDSRegisterAll());
  PetscCall(PetscFunctionListFind(PetscDSList, name, &r));
  PetscCheck(r, PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscDS type: %s", name);

  PetscTryTypeMethod(prob, destroy);
  prob->ops->destroy = NULL;

  PetscCall((*r)(prob));
  PetscCall(PetscObjectChangeTypeName((PetscObject)prob, name));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSGetType - Gets the `PetscDSType` name (as a string) from the `PetscDS`

  Not Collective; No Fortran Support

  Input Parameter:
. prob - The `PetscDS`

  Output Parameter:
. name - The `PetscDSType` name

  Level: intermediate

.seealso: `PetscDSType`, `PetscDS`, `PetscDSSetType()`, `PetscDSCreate()`
@*/
PetscErrorCode PetscDSGetType(PetscDS prob, PetscDSType *name)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
  PetscAssertPointer(name, 2);
  PetscCall(PetscDSRegisterAll());
  *name = ((PetscObject)prob)->type_name;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PetscDSView_Ascii(PetscDS ds, PetscViewer viewer)
{
  PetscViewerFormat  format;
  const PetscScalar *constants;
  PetscInt           Nf, numConstants, f;

  PetscFunctionBegin;
  PetscCall(PetscDSGetNumFields(ds, &Nf));
  PetscCall(PetscViewerGetFormat(viewer, &format));
  PetscCall(PetscViewerASCIIPrintf(viewer, "Discrete System with %" PetscInt_FMT " fields\n", Nf));
  PetscCall(PetscViewerASCIIPushTab(viewer));
  PetscCall(PetscViewerASCIIPrintf(viewer, "  cell total dim %" PetscInt_FMT " total comp %" PetscInt_FMT "\n", ds->totDim, ds->totComp));
  if (ds->isCohesive) PetscCall(PetscViewerASCIIPrintf(viewer, "  cohesive cell\n"));
  for (f = 0; f < Nf; ++f) {
    DSBoundary      b;
    PetscObject     obj;
    PetscClassId    id;
    PetscQuadrature q;
    const char     *name;
    PetscInt        Nc, Nq, Nqc;

    PetscCall(PetscDSGetDiscretization(ds, f, &obj));
    PetscCall(PetscObjectGetClassId(obj, &id));
    PetscCall(PetscObjectGetName(obj, &name));
    PetscCall(PetscViewerASCIIPrintf(viewer, "Field %s", name ? name : "<unknown>"));
    PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
    if (id == PETSCFE_CLASSID) {
      PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
      PetscCall(PetscFEGetQuadrature((PetscFE)obj, &q));
      PetscCall(PetscViewerASCIIPrintf(viewer, " FEM"));
    } else if (id == PETSCFV_CLASSID) {
      PetscCall(PetscFVGetNumComponents((PetscFV)obj, &Nc));
      PetscCall(PetscFVGetQuadrature((PetscFV)obj, &q));
      PetscCall(PetscViewerASCIIPrintf(viewer, " FVM"));
    } else SETERRQ(PetscObjectComm((PetscObject)ds), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %" PetscInt_FMT, f);
    if (Nc > 1) PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT " components", Nc));
    else PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT " component ", Nc));
    if (ds->implicit[f]) PetscCall(PetscViewerASCIIPrintf(viewer, " (implicit)"));
    else PetscCall(PetscViewerASCIIPrintf(viewer, " (explicit)"));
    if (q) {
      PetscCall(PetscQuadratureGetData(q, NULL, &Nqc, &Nq, NULL, NULL));
      PetscCall(PetscViewerASCIIPrintf(viewer, " (Nq %" PetscInt_FMT " Nqc %" PetscInt_FMT ")", Nq, Nqc));
    }
    PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT "-jet", ds->jetDegree[f]));
    PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
    PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
    PetscCall(PetscViewerASCIIPushTab(viewer));
    if (id == PETSCFE_CLASSID) PetscCall(PetscFEView((PetscFE)obj, viewer));
    else if (id == PETSCFV_CLASSID) PetscCall(PetscFVView((PetscFV)obj, viewer));
    PetscCall(PetscViewerASCIIPopTab(viewer));

    for (b = ds->boundary; b; b = b->next) {
      char    *name;
      PetscInt c, i;

      if (b->field != f) continue;
      PetscCall(PetscViewerASCIIPushTab(viewer));
      PetscCall(PetscViewerASCIIPrintf(viewer, "Boundary %s (%s) %s\n", b->name, b->lname, DMBoundaryConditionTypes[b->type]));
      if (!b->Nc) {
        PetscCall(PetscViewerASCIIPrintf(viewer, "  all components\n"));
      } else {
        PetscCall(PetscViewerASCIIPrintf(viewer, "  components: "));
        PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
        for (c = 0; c < b->Nc; ++c) {
          if (c > 0) PetscCall(PetscViewerASCIIPrintf(viewer, ", "));
          PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT, b->comps[c]));
        }
        PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
        PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
      }
      PetscCall(PetscViewerASCIIPrintf(viewer, "  values: "));
      PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
      for (i = 0; i < b->Nv; ++i) {
        if (i > 0) PetscCall(PetscViewerASCIIPrintf(viewer, ", "));
        PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT, b->values[i]));
      }
      PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
      PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
#if defined(__clang__)
      PETSC_PRAGMA_DIAGNOSTIC_IGNORED_BEGIN("-Wformat-pedantic")
#elif defined(__GNUC__) || defined(__GNUG__)
      PETSC_PRAGMA_DIAGNOSTIC_IGNORED_BEGIN("-Wformat")
#endif
      if (b->func) {
        PetscCall(PetscDLAddr(b->func, &name));
        if (name) PetscCall(PetscViewerASCIIPrintf(viewer, "  func: %s\n", name));
        else PetscCall(PetscViewerASCIIPrintf(viewer, "  func: %p\n", b->func));
        PetscCall(PetscFree(name));
      }
      if (b->func_t) {
        PetscCall(PetscDLAddr(b->func_t, &name));
        if (name) PetscCall(PetscViewerASCIIPrintf(viewer, "  func_t: %s\n", name));
        else PetscCall(PetscViewerASCIIPrintf(viewer, "  func_t: %p\n", b->func_t));
        PetscCall(PetscFree(name));
      }
      PETSC_PRAGMA_DIAGNOSTIC_IGNORED_END()
      PetscCall(PetscWeakFormView(b->wf, viewer));
      PetscCall(PetscViewerASCIIPopTab(viewer));
    }
  }
  PetscCall(PetscDSGetConstants(ds, &numConstants, &constants));
  if (numConstants) {
    PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " constants\n", numConstants));
    PetscCall(PetscViewerASCIIPushTab(viewer));
    for (f = 0; f < numConstants; ++f) PetscCall(PetscViewerASCIIPrintf(viewer, "%g\n", (double)PetscRealPart(constants[f])));
    PetscCall(PetscViewerASCIIPopTab(viewer));
  }
  PetscCall(PetscWeakFormView(ds->wf, viewer));
  PetscCall(PetscViewerASCIIPopTab(viewer));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSViewFromOptions - View a `PetscDS` based on values in the options database

  Collective

  Input Parameters:
+ A    - the `PetscDS` object
. obj  - Optional object that provides the options prefix used in the search of the options database
- name - command line option

  Level: intermediate

.seealso: `PetscDSType`, `PetscDS`, `PetscDSView()`, `PetscObjectViewFromOptions()`, `PetscDSCreate()`
@*/
PetscErrorCode PetscDSViewFromOptions(PetscDS A, PetscObject obj, const char name[])
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(A, PETSCDS_CLASSID, 1);
  PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSView - Views a `PetscDS`

  Collective

  Input Parameters:
+ prob - the `PetscDS` object to view
- v    - the viewer

  Level: developer

.seealso: `PetscDSType`, `PetscDS`, `PetscViewer`, `PetscDSDestroy()`, `PetscDSViewFromOptions()`
@*/
PetscErrorCode PetscDSView(PetscDS prob, PetscViewer v)
{
  PetscBool isascii;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
  if (!v) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)prob), &v));
  else PetscValidHeaderSpecific(v, PETSC_VIEWER_CLASSID, 2);
  PetscCall(PetscObjectTypeCompare((PetscObject)v, PETSCVIEWERASCII, &isascii));
  if (isascii) PetscCall(PetscDSView_Ascii(prob, v));
  PetscTryTypeMethod(prob, view, v);
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSSetFromOptions - sets parameters in a `PetscDS` from the options database

  Collective

  Input Parameter:
. prob - the `PetscDS` object to set options for

  Options Database Keys:
+ -petscds_type <type>     - Set the `PetscDS` type
. -petscds_view <view opt> - View the `PetscDS`
. -petscds_jac_pre         - Turn formation of a separate Jacobian preconditioner on or off
. -bc_<name> <ids>         - Specify a list of label ids for a boundary condition
- -bc_<name>_comp <comps>  - Specify a list of field components to constrain for a boundary condition

  Level: intermediate

.seealso: `PetscDS`, `PetscDSView()`
@*/
PetscErrorCode PetscDSSetFromOptions(PetscDS prob)
{
  DSBoundary  b;
  const char *defaultType;
  char        name[256];
  PetscBool   flg;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
  if (!((PetscObject)prob)->type_name) {
    defaultType = PETSCDSBASIC;
  } else {
    defaultType = ((PetscObject)prob)->type_name;
  }
  PetscCall(PetscDSRegisterAll());

  PetscObjectOptionsBegin((PetscObject)prob);
  for (b = prob->boundary; b; b = b->next) {
    char      optname[1024];
    PetscInt  ids[1024], len = 1024;
    PetscBool flg;

    PetscCall(PetscSNPrintf(optname, sizeof(optname), "-bc_%s", b->name));
    PetscCall(PetscMemzero(ids, sizeof(ids)));
    PetscCall(PetscOptionsIntArray(optname, "List of boundary IDs", "", ids, &len, &flg));
    if (flg) {
      b->Nv = len;
      PetscCall(PetscFree(b->values));
      PetscCall(PetscMalloc1(len, &b->values));
      PetscCall(PetscArraycpy(b->values, ids, len));
      PetscCall(PetscWeakFormRewriteKeys(b->wf, b->label, len, b->values));
    }
    len = 1024;
    PetscCall(PetscSNPrintf(optname, sizeof(optname), "-bc_%s_comp", b->name));
    PetscCall(PetscMemzero(ids, sizeof(ids)));
    PetscCall(PetscOptionsIntArray(optname, "List of boundary field components", "", ids, &len, &flg));
    if (flg) {
      b->Nc = len;
      PetscCall(PetscFree(b->comps));
      PetscCall(PetscMalloc1(len, &b->comps));
      PetscCall(PetscArraycpy(b->comps, ids, len));
    }
  }
  PetscCall(PetscOptionsFList("-petscds_type", "Discrete System", "PetscDSSetType", PetscDSList, defaultType, name, 256, &flg));
  if (flg) {
    PetscCall(PetscDSSetType(prob, name));
  } else if (!((PetscObject)prob)->type_name) {
    PetscCall(PetscDSSetType(prob, defaultType));
  }
  PetscCall(PetscOptionsBool("-petscds_jac_pre", "Discrete System", "PetscDSUseJacobianPreconditioner", prob->useJacPre, &prob->useJacPre, &flg));
  PetscCall(PetscOptionsBool("-petscds_force_quad", "Discrete System", "PetscDSSetForceQuad", prob->forceQuad, &prob->forceQuad, &flg));
  PetscCall(PetscOptionsInt("-petscds_print_integrate", "Discrete System", "", prob->printIntegrate, &prob->printIntegrate, NULL));
  PetscTryTypeMethod(prob, setfromoptions);
  /* process any options handlers added with PetscObjectAddOptionsHandler() */
  PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)prob, PetscOptionsObject));
  PetscOptionsEnd();
  if (prob->Nf) PetscCall(PetscDSViewFromOptions(prob, NULL, "-petscds_view"));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSSetUp - Construct data structures for the `PetscDS`

  Collective

  Input Parameter:
. prob - the `PetscDS` object to setup

  Level: developer

.seealso: `PetscDS`, `PetscDSView()`, `PetscDSDestroy()`
@*/
PetscErrorCode PetscDSSetUp(PetscDS prob)
{
  const PetscInt Nf          = prob->Nf;
  PetscBool      hasH        = PETSC_FALSE;
  PetscInt       maxOrder[4] = {-2, -2, -2, -2};
  PetscInt       dim, dimEmbed, NbMax = 0, NcMax = 0, NqMax = 0, NsMax = 1, f;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
  if (prob->setup) PetscFunctionReturn(PETSC_SUCCESS);
  /* Calculate sizes */
  PetscCall(PetscDSGetSpatialDimension(prob, &dim));
  PetscCall(PetscDSGetCoordinateDimension(prob, &dimEmbed));
  prob->totDim = prob->totComp = 0;
  PetscCall(PetscMalloc2(Nf, &prob->Nc, Nf, &prob->Nb));
  PetscCall(PetscCalloc2(Nf + 1, &prob->off, Nf + 1, &prob->offDer));
  PetscCall(PetscCalloc6(Nf + 1, &prob->offCohesive[0], Nf + 1, &prob->offCohesive[1], Nf + 1, &prob->offCohesive[2], Nf + 1, &prob->offDerCohesive[0], Nf + 1, &prob->offDerCohesive[1], Nf + 1, &prob->offDerCohesive[2]));
  PetscCall(PetscMalloc2(Nf, &prob->T, Nf, &prob->Tf));
  if (prob->forceQuad) {
    // Note: This assumes we have one kind of cell at each dimension.
    //       We can fix this by having quadrature hold the celltype
    PetscQuadrature maxQuad[4] = {NULL, NULL, NULL, NULL};

    for (f = 0; f < Nf; ++f) {
      PetscObject     obj;
      PetscClassId    id;
      PetscQuadrature q = NULL, fq = NULL;
      PetscInt        dim = -1, order = -1, forder = -1;

      PetscCall(PetscDSGetDiscretization(prob, f, &obj));
      if (!obj) continue;
      PetscCall(PetscObjectGetClassId(obj, &id));
      if (id == PETSCFE_CLASSID) {
        PetscFE fe = (PetscFE)obj;

        PetscCall(PetscFEGetQuadrature(fe, &q));
        PetscCall(PetscFEGetFaceQuadrature(fe, &fq));
      } else if (id == PETSCFV_CLASSID) {
        PetscFV fv = (PetscFV)obj;

        PetscCall(PetscFVGetQuadrature(fv, &q));
      }
      if (q) {
        PetscCall(PetscQuadratureGetData(q, &dim, NULL, NULL, NULL, NULL));
        PetscCall(PetscQuadratureGetOrder(q, &order));
        if (order > maxOrder[dim]) {
          maxOrder[dim] = order;
          maxQuad[dim]  = q;
        }
      }
      if (fq) {
        PetscCall(PetscQuadratureGetData(fq, &dim, NULL, NULL, NULL, NULL));
        PetscCall(PetscQuadratureGetOrder(fq, &forder));
        if (forder > maxOrder[dim]) {
          maxOrder[dim] = forder;
          maxQuad[dim]  = fq;
        }
      }
    }
    for (f = 0; f < Nf; ++f) {
      PetscObject     obj;
      PetscClassId    id;
      PetscQuadrature q;
      PetscInt        dim;

      PetscCall(PetscDSGetDiscretization(prob, f, &obj));
      if (!obj) continue;
      PetscCall(PetscObjectGetClassId(obj, &id));
      if (id == PETSCFE_CLASSID) {
        PetscFE fe = (PetscFE)obj;

        PetscCall(PetscFEGetQuadrature(fe, &q));
        PetscCall(PetscQuadratureGetData(q, &dim, NULL, NULL, NULL, NULL));
        PetscCall(PetscFESetQuadrature(fe, maxQuad[dim]));
        PetscCall(PetscFESetFaceQuadrature(fe, dim ? maxQuad[dim - 1] : NULL));
      } else if (id == PETSCFV_CLASSID) {
        PetscFV fv = (PetscFV)obj;

        PetscCall(PetscFVGetQuadrature(fv, &q));
        PetscCall(PetscQuadratureGetData(q, &dim, NULL, NULL, NULL, NULL));
        PetscCall(PetscFVSetQuadrature(fv, maxQuad[dim]));
      }
    }
  }
  for (f = 0; f < Nf; ++f) {
    PetscObject     obj;
    PetscClassId    id;
    PetscQuadrature q  = NULL;
    PetscInt        Nq = 0, Nb, Nc;

    PetscCall(PetscDSGetDiscretization(prob, f, &obj));
    if (prob->jetDegree[f] > 1) hasH = PETSC_TRUE;
    if (!obj) {
      /* Empty mesh */
      Nb = Nc    = 0;
      prob->T[f] = prob->Tf[f] = NULL;
    } else {
      PetscCall(PetscObjectGetClassId(obj, &id));
      if (id == PETSCFE_CLASSID) {
        PetscFE fe = (PetscFE)obj;

        PetscCall(PetscFEGetQuadrature(fe, &q));
        {
          PetscQuadrature fq;
          PetscInt        dim, order;

          PetscCall(PetscQuadratureGetData(q, &dim, NULL, NULL, NULL, NULL));
          PetscCall(PetscQuadratureGetOrder(q, &order));
          if (maxOrder[dim] < 0) maxOrder[dim] = order;
          PetscCheck(order == maxOrder[dim], PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Field %" PetscInt_FMT " cell quadrature order %" PetscInt_FMT " != %" PetscInt_FMT " DS cell quadrature order", f, order, maxOrder[dim]);
          PetscCall(PetscFEGetFaceQuadrature(fe, &fq));
          if (fq) {
            PetscCall(PetscQuadratureGetData(fq, &dim, NULL, NULL, NULL, NULL));
            PetscCall(PetscQuadratureGetOrder(fq, &order));
            if (maxOrder[dim] < 0) maxOrder[dim] = order;
            PetscCheck(order == maxOrder[dim], PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Field %" PetscInt_FMT " face quadrature order %" PetscInt_FMT " != %" PetscInt_FMT " DS face quadrature order", f, order, maxOrder[dim]);
          }
        }
        PetscCall(PetscFEGetDimension(fe, &Nb));
        PetscCall(PetscFEGetNumComponents(fe, &Nc));
        PetscCall(PetscFEGetCellTabulation(fe, prob->jetDegree[f], &prob->T[f]));
        PetscCall(PetscFEGetFaceTabulation(fe, prob->jetDegree[f], &prob->Tf[f]));
      } else if (id == PETSCFV_CLASSID) {
        PetscFV fv = (PetscFV)obj;

        PetscCall(PetscFVGetQuadrature(fv, &q));
        PetscCall(PetscFVGetNumComponents(fv, &Nc));
        Nb = Nc;
        PetscCall(PetscFVGetCellTabulation(fv, &prob->T[f]));
        /* TODO: should PetscFV also have face tabulation? Otherwise there will be a null pointer in prob->basisFace */
      } else SETERRQ(PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %" PetscInt_FMT, f);
    }
    prob->Nc[f]                    = Nc;
    prob->Nb[f]                    = Nb;
    prob->off[f + 1]               = Nc + prob->off[f];
    prob->offDer[f + 1]            = Nc * dim + prob->offDer[f];
    prob->offCohesive[0][f + 1]    = (prob->cohesive[f] ? Nc : Nc * 2) + prob->offCohesive[0][f];
    prob->offDerCohesive[0][f + 1] = (prob->cohesive[f] ? Nc : Nc * 2) * dimEmbed + prob->offDerCohesive[0][f];
    prob->offCohesive[1][f]        = (prob->cohesive[f] ? 0 : Nc) + prob->offCohesive[0][f];
    prob->offDerCohesive[1][f]     = (prob->cohesive[f] ? 0 : Nc) * dimEmbed + prob->offDerCohesive[0][f];
    prob->offCohesive[2][f + 1]    = (prob->cohesive[f] ? Nc : Nc * 2) + prob->offCohesive[2][f];
    prob->offDerCohesive[2][f + 1] = (prob->cohesive[f] ? Nc : Nc * 2) * dimEmbed + prob->offDerCohesive[2][f];
    if (q) PetscCall(PetscQuadratureGetData(q, NULL, NULL, &Nq, NULL, NULL));
    NqMax = PetscMax(NqMax, Nq);
    NbMax = PetscMax(NbMax, Nb);
    NcMax = PetscMax(NcMax, Nc);
    prob->totDim += Nb;
    prob->totComp += Nc;
    /* There are two faces for all fields on a cohesive cell, except for cohesive fields */
    if (prob->isCohesive && !prob->cohesive[f]) prob->totDim += Nb;
  }
  prob->offCohesive[1][Nf]    = prob->offCohesive[0][Nf];
  prob->offDerCohesive[1][Nf] = prob->offDerCohesive[0][Nf];
  /* Allocate works space */
  NsMax = 2; /* A non-cohesive discretizations can be used on a cohesive cell, so we need this extra workspace for all DS */
  PetscCall(PetscMalloc3(NsMax * prob->totComp, &prob->u, NsMax * prob->totComp, &prob->u_t, NsMax * prob->totComp * dimEmbed + (hasH ? NsMax * prob->totComp * dimEmbed * dimEmbed : 0), &prob->u_x));
  PetscCall(PetscMalloc5(dimEmbed, &prob->x, NbMax * NcMax, &prob->basisReal, NbMax * NcMax * dimEmbed, &prob->basisDerReal, NbMax * NcMax, &prob->testReal, NbMax * NcMax * dimEmbed, &prob->testDerReal));
  PetscCall(PetscMalloc6(NsMax * NqMax * NcMax, &prob->f0, NsMax * NqMax * NcMax * dimEmbed, &prob->f1, NsMax * NsMax * NqMax * NcMax * NcMax, &prob->g0, NsMax * NsMax * NqMax * NcMax * NcMax * dimEmbed, &prob->g1, NsMax * NsMax * NqMax * NcMax * NcMax * dimEmbed,
                         &prob->g2, NsMax * NsMax * NqMax * NcMax * NcMax * dimEmbed * dimEmbed, &prob->g3));
  PetscTryTypeMethod(prob, setup);
  prob->setup = PETSC_TRUE;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PetscDSDestroyStructs_Static(PetscDS prob)
{
  PetscFunctionBegin;
  PetscCall(PetscFree2(prob->Nc, prob->Nb));
  PetscCall(PetscFree2(prob->off, prob->offDer));
  PetscCall(PetscFree6(prob->offCohesive[0], prob->offCohesive[1], prob->offCohesive[2], prob->offDerCohesive[0], prob->offDerCohesive[1], prob->offDerCohesive[2]));
  PetscCall(PetscFree2(prob->T, prob->Tf));
  PetscCall(PetscFree3(prob->u, prob->u_t, prob->u_x));
  PetscCall(PetscFree5(prob->x, prob->basisReal, prob->basisDerReal, prob->testReal, prob->testDerReal));
  PetscCall(PetscFree6(prob->f0, prob->f1, prob->g0, prob->g1, prob->g2, prob->g3));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PetscDSEnlarge_Static(PetscDS prob, PetscInt NfNew)
{
  PetscObject         *tmpd;
  PetscBool           *tmpi;
  PetscInt            *tmpk;
  PetscBool           *tmpc;
  PetscPointFn       **tmpup;
  PetscSimplePointFn **tmpexactSol, **tmpexactSol_t, **tmplowerBound, **tmpupperBound;
  void               **tmpexactCtx, **tmpexactCtx_t, **tmplowerCtx, **tmpupperCtx;
  void               **tmpctx;
  PetscInt             Nf = prob->Nf, f;

  PetscFunctionBegin;
  if (Nf >= NfNew) PetscFunctionReturn(PETSC_SUCCESS);
  prob->setup = PETSC_FALSE;
  PetscCall(PetscDSDestroyStructs_Static(prob));
  PetscCall(PetscMalloc4(NfNew, &tmpd, NfNew, &tmpi, NfNew, &tmpc, NfNew, &tmpk));
  for (f = 0; f < Nf; ++f) {
    tmpd[f] = prob->disc[f];
    tmpi[f] = prob->implicit[f];
    tmpc[f] = prob->cohesive[f];
    tmpk[f] = prob->jetDegree[f];
  }
  for (f = Nf; f < NfNew; ++f) {
    tmpd[f] = NULL;
    tmpi[f] = PETSC_TRUE, tmpc[f] = PETSC_FALSE;
    tmpk[f] = 1;
  }
  PetscCall(PetscFree4(prob->disc, prob->implicit, prob->cohesive, prob->jetDegree));
  PetscCall(PetscWeakFormSetNumFields(prob->wf, NfNew));
  prob->Nf        = NfNew;
  prob->disc      = tmpd;
  prob->implicit  = tmpi;
  prob->cohesive  = tmpc;
  prob->jetDegree = tmpk;
  PetscCall(PetscCalloc2(NfNew, &tmpup, NfNew, &tmpctx));
  for (f = 0; f < Nf; ++f) tmpup[f] = prob->update[f];
  for (f = 0; f < Nf; ++f) tmpctx[f] = prob->ctx[f];
  for (f = Nf; f < NfNew; ++f) tmpup[f] = NULL;
  for (f = Nf; f < NfNew; ++f) tmpctx[f] = NULL;
  PetscCall(PetscFree2(prob->update, prob->ctx));
  prob->update = tmpup;
  prob->ctx    = tmpctx;
  PetscCall(PetscCalloc4(NfNew, &tmpexactSol, NfNew, &tmpexactCtx, NfNew, &tmpexactSol_t, NfNew, &tmpexactCtx_t));
  PetscCall(PetscCalloc4(NfNew, &tmplowerBound, NfNew, &tmplowerCtx, NfNew, &tmpupperBound, NfNew, &tmpupperCtx));
  for (f = 0; f < Nf; ++f) tmpexactSol[f] = prob->exactSol[f];
  for (f = 0; f < Nf; ++f) tmpexactCtx[f] = prob->exactCtx[f];
  for (f = 0; f < Nf; ++f) tmpexactSol_t[f] = prob->exactSol_t[f];
  for (f = 0; f < Nf; ++f) tmpexactCtx_t[f] = prob->exactCtx_t[f];
  for (f = 0; f < Nf; ++f) tmplowerBound[f] = prob->lowerBound[f];
  for (f = 0; f < Nf; ++f) tmplowerCtx[f] = prob->lowerCtx[f];
  for (f = 0; f < Nf; ++f) tmpupperBound[f] = prob->upperBound[f];
  for (f = 0; f < Nf; ++f) tmpupperCtx[f] = prob->upperCtx[f];
  for (f = Nf; f < NfNew; ++f) tmpexactSol[f] = NULL;
  for (f = Nf; f < NfNew; ++f) tmpexactCtx[f] = NULL;
  for (f = Nf; f < NfNew; ++f) tmpexactSol_t[f] = NULL;
  for (f = Nf; f < NfNew; ++f) tmpexactCtx_t[f] = NULL;
  for (f = Nf; f < NfNew; ++f) tmplowerBound[f] = NULL;
  for (f = Nf; f < NfNew; ++f) tmplowerCtx[f] = NULL;
  for (f = Nf; f < NfNew; ++f) tmpupperBound[f] = NULL;
  for (f = Nf; f < NfNew; ++f) tmpupperCtx[f] = NULL;
  PetscCall(PetscFree4(prob->exactSol, prob->exactCtx, prob->exactSol_t, prob->exactCtx_t));
  PetscCall(PetscFree4(prob->lowerBound, prob->lowerCtx, prob->upperBound, prob->upperCtx));
  prob->exactSol   = tmpexactSol;
  prob->exactCtx   = tmpexactCtx;
  prob->exactSol_t = tmpexactSol_t;
  prob->exactCtx_t = tmpexactCtx_t;
  prob->lowerBound = tmplowerBound;
  prob->lowerCtx   = tmplowerCtx;
  prob->upperBound = tmpupperBound;
  prob->upperCtx   = tmpupperCtx;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSDestroy - Destroys a `PetscDS` object

  Collective

  Input Parameter:
. ds - the `PetscDS` object to destroy

  Level: developer

.seealso: `PetscDSView()`
@*/
PetscErrorCode PetscDSDestroy(PetscDS *ds)
{
  PetscInt f;

  PetscFunctionBegin;
  if (!*ds) PetscFunctionReturn(PETSC_SUCCESS);
  PetscValidHeaderSpecific(*ds, PETSCDS_CLASSID, 1);

  if (--((PetscObject)*ds)->refct > 0) {
    *ds = NULL;
    PetscFunctionReturn(PETSC_SUCCESS);
  }
  ((PetscObject)*ds)->refct = 0;
  if ((*ds)->subprobs) {
    PetscInt dim, d;

    PetscCall(PetscDSGetSpatialDimension(*ds, &dim));
    for (d = 0; d < dim; ++d) PetscCall(PetscDSDestroy(&(*ds)->subprobs[d]));
  }
  PetscCall(PetscFree((*ds)->subprobs));
  PetscCall(PetscDSDestroyStructs_Static(*ds));
  for (f = 0; f < (*ds)->Nf; ++f) PetscCall(PetscObjectDereference((*ds)->disc[f]));
  PetscCall(PetscFree4((*ds)->disc, (*ds)->implicit, (*ds)->cohesive, (*ds)->jetDegree));
  PetscCall(PetscWeakFormDestroy(&(*ds)->wf));
  PetscCall(PetscFree2((*ds)->update, (*ds)->ctx));
  PetscCall(PetscFree4((*ds)->exactSol, (*ds)->exactCtx, (*ds)->exactSol_t, (*ds)->exactCtx_t));
  PetscCall(PetscFree4((*ds)->lowerBound, (*ds)->lowerCtx, (*ds)->upperBound, (*ds)->upperCtx));
  PetscTryTypeMethod(*ds, destroy);
  PetscCall(PetscDSDestroyBoundary(*ds));
  PetscCall(PetscFree((*ds)->constants));
  for (PetscInt c = 0; c < DM_NUM_POLYTOPES; ++c) {
    const PetscInt Na = DMPolytopeTypeGetNumArrangements((DMPolytopeType)c);
    if ((*ds)->quadPerm[c])
      for (PetscInt o = 0; o < Na; ++o) PetscCall(ISDestroy(&(*ds)->quadPerm[c][o]));
    PetscCall(PetscFree((*ds)->quadPerm[c]));
    (*ds)->quadPerm[c] = NULL;
  }
  PetscCall(PetscHeaderDestroy(ds));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSCreate - Creates an empty `PetscDS` object. The type can then be set with `PetscDSSetType()`.

  Collective

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

  Output Parameter:
. ds - The `PetscDS` object

  Level: beginner

.seealso: `PetscDS`, `PetscDSSetType()`, `PETSCDSBASIC`, `PetscDSType`, `PetscDSDestroy()`
@*/
PetscErrorCode PetscDSCreate(MPI_Comm comm, PetscDS *ds)
{
  PetscDS p;

  PetscFunctionBegin;
  PetscAssertPointer(ds, 2);
  PetscCall(PetscDSInitializePackage());

  PetscCall(PetscHeaderCreate(p, PETSCDS_CLASSID, "PetscDS", "Discrete System", "PetscDS", comm, PetscDSDestroy, PetscDSView));
  p->Nf               = 0;
  p->setup            = PETSC_FALSE;
  p->numConstants     = 0;
  p->numFuncConstants = 3; // Row and col fields, cell size
  p->dimEmbed         = -1;
  p->useJacPre        = PETSC_TRUE;
  p->forceQuad        = PETSC_TRUE;
  PetscCall(PetscMalloc1(p->numConstants + p->numFuncConstants, &p->constants));
  PetscCall(PetscWeakFormCreate(comm, &p->wf));
  PetscCall(PetscArrayzero(p->quadPerm, DM_NUM_POLYTOPES));
  *ds = p;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSGetNumFields - Returns the number of fields in the `PetscDS`

  Not Collective

  Input Parameter:
. prob - The `PetscDS` object

  Output Parameter:
. Nf - The number of fields

  Level: beginner

.seealso: `PetscDS`, `PetscDSGetSpatialDimension()`, `PetscDSCreate()`
@*/
PetscErrorCode PetscDSGetNumFields(PetscDS prob, PetscInt *Nf)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
  PetscAssertPointer(Nf, 2);
  *Nf = prob->Nf;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSGetSpatialDimension - Returns the spatial dimension of the `PetscDS`, meaning the topological dimension of the discretizations

  Not Collective

  Input Parameter:
. prob - The `PetscDS` object

  Output Parameter:
. dim - The spatial dimension

  Level: beginner

.seealso: `PetscDS`, `PetscDSGetCoordinateDimension()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
@*/
PetscErrorCode PetscDSGetSpatialDimension(PetscDS prob, PetscInt *dim)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
  PetscAssertPointer(dim, 2);
  *dim = 0;
  if (prob->Nf) {
    PetscObject  obj;
    PetscClassId id;

    PetscCall(PetscDSGetDiscretization(prob, 0, &obj));
    if (obj) {
      PetscCall(PetscObjectGetClassId(obj, &id));
      if (id == PETSCFE_CLASSID) PetscCall(PetscFEGetSpatialDimension((PetscFE)obj, dim));
      else if (id == PETSCFV_CLASSID) PetscCall(PetscFVGetSpatialDimension((PetscFV)obj, dim));
      else SETERRQ(PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", 0);
    }
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSGetCoordinateDimension - Returns the coordinate dimension of the `PetscDS`, meaning the dimension of the space into which the discretiaztions are embedded

  Not Collective

  Input Parameter:
. prob - The `PetscDS` object

  Output Parameter:
. dimEmbed - The coordinate dimension

  Level: beginner

.seealso: `PetscDS`, `PetscDSSetCoordinateDimension()`, `PetscDSGetSpatialDimension()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
@*/
PetscErrorCode PetscDSGetCoordinateDimension(PetscDS prob, PetscInt *dimEmbed)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
  PetscAssertPointer(dimEmbed, 2);
  PetscCheck(prob->dimEmbed >= 0, PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_WRONGSTATE, "No coordinate dimension set for this DS");
  *dimEmbed = prob->dimEmbed;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSSetCoordinateDimension - Set the coordinate dimension of the `PetscDS`, meaning the dimension of the space into which the discretiaztions are embedded

  Logically Collective

  Input Parameters:
+ prob     - The `PetscDS` object
- dimEmbed - The coordinate dimension

  Level: beginner

.seealso: `PetscDS`, `PetscDSGetCoordinateDimension()`, `PetscDSGetSpatialDimension()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
@*/
PetscErrorCode PetscDSSetCoordinateDimension(PetscDS prob, PetscInt dimEmbed)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
  PetscCheck(dimEmbed >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Coordinate dimension must be non-negative, not %" PetscInt_FMT, dimEmbed);
  prob->dimEmbed = dimEmbed;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSGetForceQuad - Returns the flag to force matching quadratures among the field discretizations

  Not collective

  Input Parameter:
. ds - The `PetscDS` object

  Output Parameter:
. forceQuad - The flag

  Level: intermediate

.seealso: `PetscDS`, `PetscDSSetForceQuad()`, `PetscDSGetDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
@*/
PetscErrorCode PetscDSGetForceQuad(PetscDS ds, PetscBool *forceQuad)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  PetscAssertPointer(forceQuad, 2);
  *forceQuad = ds->forceQuad;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSSetForceQuad - Set the flag to force matching quadratures among the field discretizations

  Logically collective on ds

  Input Parameters:
+ ds        - The `PetscDS` object
- forceQuad - The flag

  Level: intermediate

.seealso: `PetscDS`, `PetscDSGetForceQuad()`, `PetscDSGetDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
@*/
PetscErrorCode PetscDSSetForceQuad(PetscDS ds, PetscBool forceQuad)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  ds->forceQuad = forceQuad;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSIsCohesive - Returns the flag indicating that this `PetscDS` is for a cohesive cell

  Not Collective

  Input Parameter:
. ds - The `PetscDS` object

  Output Parameter:
. isCohesive - The flag

  Level: developer

.seealso: `PetscDS`, `PetscDSGetNumCohesive()`, `PetscDSGetCohesive()`, `PetscDSSetCohesive()`, `PetscDSCreate()`
@*/
PetscErrorCode PetscDSIsCohesive(PetscDS ds, PetscBool *isCohesive)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  PetscAssertPointer(isCohesive, 2);
  *isCohesive = ds->isCohesive;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSGetNumCohesive - Returns the number of cohesive fields, meaning those defined on the interior of a cohesive cell

  Not Collective

  Input Parameter:
. ds - The `PetscDS` object

  Output Parameter:
. numCohesive - The number of cohesive fields

  Level: developer

.seealso: `PetscDS`, `PetscDSSetCohesive()`, `PetscDSCreate()`
@*/
PetscErrorCode PetscDSGetNumCohesive(PetscDS ds, PetscInt *numCohesive)
{
  PetscInt f;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  PetscAssertPointer(numCohesive, 2);
  *numCohesive = 0;
  for (f = 0; f < ds->Nf; ++f) *numCohesive += ds->cohesive[f] ? 1 : 0;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSGetCohesive - Returns the flag indicating that a field is cohesive, meaning it is defined on the interior of a cohesive cell

  Not Collective

  Input Parameters:
+ ds - The `PetscDS` object
- f  - The field index

  Output Parameter:
. isCohesive - The flag

  Level: developer

.seealso: `PetscDS`, `PetscDSSetCohesive()`, `PetscDSIsCohesive()`, `PetscDSCreate()`
@*/
PetscErrorCode PetscDSGetCohesive(PetscDS ds, PetscInt f, PetscBool *isCohesive)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  PetscAssertPointer(isCohesive, 3);
  PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
  *isCohesive = ds->cohesive[f];
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSSetCohesive - Set the flag indicating that a field is cohesive, meaning it is defined on the interior of a cohesive cell

  Not Collective

  Input Parameters:
+ ds         - The `PetscDS` object
. f          - The field index
- isCohesive - The flag for a cohesive field

  Level: developer

.seealso: `PetscDS`, `PetscDSGetCohesive()`, `PetscDSIsCohesive()`, `PetscDSCreate()`
@*/
PetscErrorCode PetscDSSetCohesive(PetscDS ds, PetscInt f, PetscBool isCohesive)
{
  PetscInt i;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
  ds->cohesive[f] = isCohesive;
  ds->isCohesive  = PETSC_FALSE;
  for (i = 0; i < ds->Nf; ++i) ds->isCohesive = ds->isCohesive || ds->cohesive[f] ? PETSC_TRUE : PETSC_FALSE;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSGetTotalDimension - Returns the total size of the approximation space for this system

  Not Collective

  Input Parameter:
. prob - The `PetscDS` object

  Output Parameter:
. dim - The total problem dimension

  Level: beginner

.seealso: `PetscDS`, `PetscDSGetNumFields()`, `PetscDSCreate()`
@*/
PetscErrorCode PetscDSGetTotalDimension(PetscDS prob, PetscInt *dim)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
  PetscCall(PetscDSSetUp(prob));
  PetscAssertPointer(dim, 2);
  *dim = prob->totDim;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSGetTotalComponents - Returns the total number of components in this system

  Not Collective

  Input Parameter:
. prob - The `PetscDS` object

  Output Parameter:
. Nc - The total number of components

  Level: beginner

.seealso: `PetscDS`, `PetscDSGetNumFields()`, `PetscDSCreate()`
@*/
PetscErrorCode PetscDSGetTotalComponents(PetscDS prob, PetscInt *Nc)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
  PetscCall(PetscDSSetUp(prob));
  PetscAssertPointer(Nc, 2);
  *Nc = prob->totComp;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSGetDiscretization - Returns the discretization object for the given field

  Not Collective

  Input Parameters:
+ prob - The `PetscDS` object
- f    - The field number

  Output Parameter:
. disc - The discretization object, this can be a `PetscFE` or a `PetscFV`

  Level: beginner

.seealso: `PetscDS`, `PetscFE`, `PetscFV`, `PetscDSSetDiscretization()`, `PetscDSAddDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
@*/
PetscErrorCode PetscDSGetDiscretization(PetscDS prob, PetscInt f, PetscObject *disc)
{
  PetscFunctionBeginHot;
  PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
  PetscAssertPointer(disc, 3);
  PetscCheck(!(f < 0) && !(f >= prob->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, prob->Nf);
  *disc = prob->disc[f];
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSSetDiscretization - Sets the discretization object for the given field

  Not Collective

  Input Parameters:
+ prob - The `PetscDS` object
. f    - The field number
- disc - The discretization object, this can be a `PetscFE` or a `PetscFV`

  Level: beginner

.seealso: `PetscDS`, `PetscFE`, `PetscFV`, `PetscDSGetDiscretization()`, `PetscDSAddDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
@*/
PetscErrorCode PetscDSSetDiscretization(PetscDS prob, PetscInt f, PetscObject disc)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
  if (disc) PetscAssertPointer(disc, 3);
  PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
  PetscCall(PetscDSEnlarge_Static(prob, f + 1));
  PetscCall(PetscObjectDereference(prob->disc[f]));
  prob->disc[f] = disc;
  PetscCall(PetscObjectReference(disc));
  if (disc) {
    PetscClassId id;

    PetscCall(PetscObjectGetClassId(disc, &id));
    if (id == PETSCFE_CLASSID) {
      PetscCall(PetscDSSetImplicit(prob, f, PETSC_TRUE));
    } else if (id == PETSCFV_CLASSID) {
      PetscCall(PetscDSSetImplicit(prob, f, PETSC_FALSE));
    }
    PetscCall(PetscDSSetJetDegree(prob, f, 1));
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSGetWeakForm - Returns the weak form object from within the `PetscDS`

  Not Collective

  Input Parameter:
. ds - The `PetscDS` object

  Output Parameter:
. wf - The weak form object

  Level: beginner

.seealso: `PetscWeakForm`, `PetscDSSetWeakForm()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
@*/
PetscErrorCode PetscDSGetWeakForm(PetscDS ds, PetscWeakForm *wf)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  PetscAssertPointer(wf, 2);
  *wf = ds->wf;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSSetWeakForm - Sets the weak form object to be used by the `PetscDS`

  Not Collective

  Input Parameters:
+ ds - The `PetscDS` object
- wf - The weak form object

  Level: beginner

.seealso: `PetscWeakForm`, `PetscDSGetWeakForm()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
@*/
PetscErrorCode PetscDSSetWeakForm(PetscDS ds, PetscWeakForm wf)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  PetscValidHeaderSpecific(wf, PETSCWEAKFORM_CLASSID, 2);
  PetscCall(PetscObjectDereference((PetscObject)ds->wf));
  ds->wf = wf;
  PetscCall(PetscObjectReference((PetscObject)wf));
  PetscCall(PetscWeakFormSetNumFields(wf, ds->Nf));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSAddDiscretization - Adds a discretization object

  Not Collective

  Input Parameters:
+ prob - The `PetscDS` object
- disc - The discretization object, this can be a `PetscFE` or `PetscFV`

  Level: beginner

.seealso: `PetscWeakForm`, `PetscFE`, `PetscFV`, `PetscDSGetDiscretization()`, `PetscDSSetDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
@*/
PetscErrorCode PetscDSAddDiscretization(PetscDS prob, PetscObject disc)
{
  PetscFunctionBegin;
  PetscCall(PetscDSSetDiscretization(prob, prob->Nf, disc));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSGetQuadrature - Returns the quadrature, which must agree for all fields in the `PetscDS`

  Not Collective

  Input Parameter:
. prob - The `PetscDS` object

  Output Parameter:
. q - The quadrature object

  Level: intermediate

.seealso: `PetscDS`, `PetscQuadrature`, `PetscDSSetImplicit()`, `PetscDSSetDiscretization()`, `PetscDSAddDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
@*/
PetscErrorCode PetscDSGetQuadrature(PetscDS prob, PetscQuadrature *q)
{
  PetscObject  obj;
  PetscClassId id;

  PetscFunctionBegin;
  *q = NULL;
  if (!prob->Nf) PetscFunctionReturn(PETSC_SUCCESS);
  PetscCall(PetscDSGetDiscretization(prob, 0, &obj));
  PetscCall(PetscObjectGetClassId(obj, &id));
  if (id == PETSCFE_CLASSID) PetscCall(PetscFEGetQuadrature((PetscFE)obj, q));
  else if (id == PETSCFV_CLASSID) PetscCall(PetscFVGetQuadrature((PetscFV)obj, q));
  else SETERRQ(PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", 0);
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSGetImplicit - Returns the flag for implicit solve for this field. This is just a guide for `TSARKIMEX`

  Not Collective

  Input Parameters:
+ prob - The `PetscDS` object
- f    - The field number

  Output Parameter:
. implicit - The flag indicating what kind of solve to use for this field

  Level: developer

.seealso: `TSARKIMEX`, `PetscDS`, `PetscDSSetImplicit()`, `PetscDSSetDiscretization()`, `PetscDSAddDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
@*/
PetscErrorCode PetscDSGetImplicit(PetscDS prob, PetscInt f, PetscBool *implicit)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
  PetscAssertPointer(implicit, 3);
  PetscCheck(!(f < 0) && !(f >= prob->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, prob->Nf);
  *implicit = prob->implicit[f];
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSSetImplicit - Set the flag for implicit solve for this field. This is just a guide for `TSARKIMEX`

  Not Collective

  Input Parameters:
+ prob     - The `PetscDS` object
. f        - The field number
- implicit - The flag indicating what kind of solve to use for this field

  Level: developer

.seealso: `TSARKIMEX`, `PetscDSGetImplicit()`, `PetscDSSetDiscretization()`, `PetscDSAddDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
@*/
PetscErrorCode PetscDSSetImplicit(PetscDS prob, PetscInt f, PetscBool implicit)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
  PetscCheck(!(f < 0) && !(f >= prob->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, prob->Nf);
  prob->implicit[f] = implicit;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSGetJetDegree - Returns the highest derivative for this field equation, or the k-jet that the discretization needs to tabulate.

  Not Collective

  Input Parameters:
+ ds - The `PetscDS` object
- f  - The field number

  Output Parameter:
. k - The highest derivative we need to tabulate

  Level: developer

.seealso: `PetscDS`, `PetscDSSetJetDegree()`, `PetscDSSetDiscretization()`, `PetscDSAddDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
@*/
PetscErrorCode PetscDSGetJetDegree(PetscDS ds, PetscInt f, PetscInt *k)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  PetscAssertPointer(k, 3);
  PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
  *k = ds->jetDegree[f];
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSSetJetDegree - Set the highest derivative for this field equation, or the k-jet that the discretization needs to tabulate.

  Not Collective

  Input Parameters:
+ ds - The `PetscDS` object
. f  - The field number
- k  - The highest derivative we need to tabulate

  Level: developer

.seealso: `PetscDS`, `PetscDSGetJetDegree()`, `PetscDSSetDiscretization()`, `PetscDSAddDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
@*/
PetscErrorCode PetscDSSetJetDegree(PetscDS ds, PetscInt f, PetscInt k)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
  ds->jetDegree[f] = k;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  PetscDSGetObjective - Get the pointwise objective function for a given test field that was provided with `PetscDSSetObjective()`

  Not Collective

  Input Parameters:
+ ds - The `PetscDS`
- f  - The test field number

  Output Parameter:
. obj - integrand for the test function term, see `PetscPointFn`

  Level: intermediate

  Note:
  We are using a first order FEM model for the weak form\: $  \int_\Omega \phi\,\mathrm{obj}(u, u_t, \nabla u, x, t)$

.seealso: `PetscPointFn`, `PetscDS`, `PetscDSSetObjective()`, `PetscDSGetResidual()`
@*/
PetscErrorCode PetscDSGetObjective(PetscDS ds, PetscInt f, PetscPointFn **obj)
{
  PetscPointFn **tmp;
  PetscInt       n;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  PetscAssertPointer(obj, 3);
  PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
  PetscCall(PetscWeakFormGetObjective(ds->wf, NULL, 0, f, 0, &n, &tmp));
  *obj = tmp ? tmp[0] : NULL;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  PetscDSSetObjective - Set the pointwise objective function for a given test field

  Not Collective

  Input Parameters:
+ ds  - The `PetscDS`
. f   - The test field number
- obj - integrand for the test function term, see `PetscPointFn`

  Level: intermediate

  Note:
  We are using a first order FEM model for the weak form\: $  \int_\Omega \phi\,\mathrm{obj}(u, u_t, \nabla u, x, t)$

.seealso: `PetscPointFn`, `PetscDS`, `PetscDSGetObjective()`, `PetscDSSetResidual()`
@*/
PetscErrorCode PetscDSSetObjective(PetscDS ds, PetscInt f, PetscPointFn *obj)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  if (obj) PetscValidFunction(obj, 3);
  PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
  PetscCall(PetscWeakFormSetIndexObjective(ds->wf, NULL, 0, f, 0, 0, obj));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  PetscDSGetResidual - Get the pointwise residual function for a given test field

  Not Collective

  Input Parameters:
+ ds - The `PetscDS`
- f  - The test field number

  Output Parameters:
+ f0 - integrand for the test function term, see `PetscPointFn`
- f1 - integrand for the test function gradient term, see `PetscPointFn`

  Level: intermediate

  Note:
  We are using a first order FEM model for the weak form\: $  \int_\Omega \phi f_0(u, u_t, \nabla u, x, t) + \nabla\phi \cdot {\vec f}_1(u, u_t, \nabla u, x, t)$

.seealso: `PetscPointFn`, `PetscDS`, `PetscDSSetResidual()`
@*/
PetscErrorCode PetscDSGetResidual(PetscDS ds, PetscInt f, PetscPointFn **f0, PetscPointFn **f1)
{
  PetscPointFn **tmp0, **tmp1;
  PetscInt       n0, n1;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
  PetscCall(PetscWeakFormGetResidual(ds->wf, NULL, 0, f, 0, &n0, &tmp0, &n1, &tmp1));
  *f0 = tmp0 ? tmp0[0] : NULL;
  *f1 = tmp1 ? tmp1[0] : NULL;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  PetscDSSetResidual - Set the pointwise residual function for a given test field

  Not Collective

  Input Parameters:
+ ds - The `PetscDS`
. f  - The test field number
. f0 - integrand for the test function term, see `PetscPointFn`
- f1 - integrand for the test function gradient term, see `PetscPointFn`

  Level: intermediate

  Note:
  We are using a first order FEM model for the weak form\: $  \int_\Omega \phi f_0(u, u_t, \nabla u, x, t) + \nabla\phi \cdot {\vec f}_1(u, u_t, \nabla u, x, t)$

.seealso: `PetscPointFn`, `PetscDS`, `PetscDSGetResidual()`
@*/
PetscErrorCode PetscDSSetResidual(PetscDS ds, PetscInt f, PetscPointFn *f0, PetscPointFn *f1)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  if (f0) PetscValidFunction(f0, 3);
  if (f1) PetscValidFunction(f1, 4);
  PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
  PetscCall(PetscWeakFormSetIndexResidual(ds->wf, NULL, 0, f, 0, 0, f0, 0, f1));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  PetscDSGetRHSResidual - Get the pointwise RHS residual function for explicit timestepping for a given test field

  Not Collective

  Input Parameters:
+ ds - The `PetscDS`
- f  - The test field number

  Output Parameters:
+ f0 - integrand for the test function term, see `PetscPointFn`
- f1 - integrand for the test function gradient term, see `PetscPointFn`

  Level: intermediate

  Note:
  We are using a first order FEM model for the weak form\: $ \int_\Omega \phi f_0(u, u_t, \nabla u, x, t) + \nabla\phi \cdot {\vec f}_1(u, u_t, \nabla u, x, t)$

.seealso: `PetscPointFn`, `PetscDS`, `PetscDSSetRHSResidual()`
@*/
PetscErrorCode PetscDSGetRHSResidual(PetscDS ds, PetscInt f, PetscPointFn **f0, PetscPointFn **f1)
{
  PetscPointFn **tmp0, **tmp1;
  PetscInt       n0, n1;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
  PetscCall(PetscWeakFormGetResidual(ds->wf, NULL, 0, f, 100, &n0, &tmp0, &n1, &tmp1));
  *f0 = tmp0 ? tmp0[0] : NULL;
  *f1 = tmp1 ? tmp1[0] : NULL;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  PetscDSSetRHSResidual - Set the pointwise residual function for explicit timestepping for a given test field

  Not Collective

  Input Parameters:
+ ds - The `PetscDS`
. f  - The test field number
. f0 - integrand for the test function term, see `PetscPointFn`
- f1 - integrand for the test function gradient term, see `PetscPointFn`

  Level: intermediate

  Note:
  We are using a first order FEM model for the weak form\: $ \int_\Omega \phi f_0(u, u_t, \nabla u, x, t) + \nabla\phi \cdot {\vec f}_1(u, u_t, \nabla u, x, t)$

.seealso: `PetscDS`, `PetscDSGetResidual()`
@*/
PetscErrorCode PetscDSSetRHSResidual(PetscDS ds, PetscInt f, PetscPointFn *f0, PetscPointFn *f1)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  if (f0) PetscValidFunction(f0, 3);
  if (f1) PetscValidFunction(f1, 4);
  PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
  PetscCall(PetscWeakFormSetIndexResidual(ds->wf, NULL, 0, f, 100, 0, f0, 0, f1));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSHasJacobian - Checks that the Jacobian functions have been set

  Not Collective

  Input Parameter:
. ds - The `PetscDS`

  Output Parameter:
. hasJac - flag that indicates the pointwise function for the Jacobian has been set

  Level: intermediate

.seealso: `PetscDS`, `PetscDSGetJacobianPreconditioner()`, `PetscDSSetJacobianPreconditioner()`, `PetscDSGetJacobian()`
@*/
PetscErrorCode PetscDSHasJacobian(PetscDS ds, PetscBool *hasJac)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  PetscCall(PetscWeakFormHasJacobian(ds->wf, hasJac));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  PetscDSGetJacobian - Get the pointwise Jacobian function for given test and basis field

  Not Collective

  Input Parameters:
+ ds - The `PetscDS`
. f  - The test field number
- g  - The field number

  Output Parameters:
+ g0 - integrand for the test and basis function term, see `PetscPointJacFn`
. g1 - integrand for the test function and basis function gradient term, see `PetscPointJacFn`
. g2 - integrand for the test function gradient and basis function term, see `PetscPointJacFn`
- g3 - integrand for the test function gradient and basis function gradient term, see `PetscPointJacFn`

  Level: intermediate

  Note:
  We are using a first order FEM model for the weak form\:

  $$
  \int_\Omega \phi\, g_0(u, u_t, \nabla u, x, t) \psi + \phi\, {\vec g}_1(u, u_t, \nabla u, x, t) \nabla \psi
  + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \nabla \psi
  $$

.seealso: `PetscDS`, `PetscDSSetJacobian()`, `PetscPointJacFn`
@*/
PetscErrorCode PetscDSGetJacobian(PetscDS ds, PetscInt f, PetscInt g, PetscPointJacFn **g0, PetscPointJacFn **g1, PetscPointJacFn **g2, PetscPointJacFn **g3)
{
  PetscPointJacFn **tmp0, **tmp1, **tmp2, **tmp3;
  PetscInt          n0, n1, n2, n3;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
  PetscCheck(!(g < 0) && !(g >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", g, ds->Nf);
  PetscCall(PetscWeakFormGetJacobian(ds->wf, NULL, 0, f, g, 0, &n0, &tmp0, &n1, &tmp1, &n2, &tmp2, &n3, &tmp3));
  *g0 = tmp0 ? tmp0[0] : NULL;
  *g1 = tmp1 ? tmp1[0] : NULL;
  *g2 = tmp2 ? tmp2[0] : NULL;
  *g3 = tmp3 ? tmp3[0] : NULL;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  PetscDSSetJacobian - Set the pointwise Jacobian function for given test and basis fields

  Not Collective

  Input Parameters:
+ ds - The `PetscDS`
. f  - The test field number
. g  - The field number
. g0 - integrand for the test and basis function term, see `PetscPointJacFn`
. g1 - integrand for the test function and basis function gradient term, see `PetscPointJacFn`
. g2 - integrand for the test function gradient and basis function term, see `PetscPointJacFn`
- g3 - integrand for the test function gradient and basis function gradient term, see `PetscPointJacFn`

  Level: intermediate

  Note:
  We are using a first order FEM model for the weak form\:

  $$
  \int_\Omega \phi\, g_0(u, u_t, \nabla u, x, t) \psi + \phi\, {\vec g}_1(u, u_t, \nabla u, x, t) \nabla \psi
  + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \nabla \psi
  $$

.seealso: `PetscDS`, `PetscDSGetJacobian()`, `PetscPointJacFn`
@*/
PetscErrorCode PetscDSSetJacobian(PetscDS ds, PetscInt f, PetscInt g, PetscPointJacFn *g0, PetscPointJacFn *g1, PetscPointJacFn *g2, PetscPointJacFn *g3)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  if (g0) PetscValidFunction(g0, 4);
  if (g1) PetscValidFunction(g1, 5);
  if (g2) PetscValidFunction(g2, 6);
  if (g3) PetscValidFunction(g3, 7);
  PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
  PetscCheck(g >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", g);
  PetscCall(PetscWeakFormSetIndexJacobian(ds->wf, NULL, 0, f, g, 0, 0, g0, 0, g1, 0, g2, 0, g3));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSUseJacobianPreconditioner - Set whether to construct a Jacobian preconditioner

  Not Collective

  Input Parameters:
+ prob      - The `PetscDS`
- useJacPre - flag that enables construction of a Jacobian preconditioner

  Level: intermediate

  Developer Note:
  Should be called `PetscDSSetUseJacobianPreconditioner()`

.seealso: `PetscDS`, `PetscDSGetJacobianPreconditioner()`, `PetscDSSetJacobianPreconditioner()`, `PetscDSGetJacobian()`
@*/
PetscErrorCode PetscDSUseJacobianPreconditioner(PetscDS prob, PetscBool useJacPre)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
  prob->useJacPre = useJacPre;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSHasJacobianPreconditioner - Checks if a Jacobian matrix for constructing a preconditioner has been set

  Not Collective

  Input Parameter:
. ds - The `PetscDS`

  Output Parameter:
. hasJacPre - the flag

  Level: intermediate

.seealso: `PetscDS`, `PetscDSGetJacobianPreconditioner()`, `PetscDSSetJacobianPreconditioner()`, `PetscDSGetJacobian()`
@*/
PetscErrorCode PetscDSHasJacobianPreconditioner(PetscDS ds, PetscBool *hasJacPre)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  *hasJacPre = PETSC_FALSE;
  if (!ds->useJacPre) PetscFunctionReturn(PETSC_SUCCESS);
  PetscCall(PetscWeakFormHasJacobianPreconditioner(ds->wf, hasJacPre));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  PetscDSGetJacobianPreconditioner - Get the pointwise Jacobian function for given test and basis field that constructs the matrix used
  to compute the preconditioner. If this is missing, the system matrix is used to build the preconditioner.

  Not Collective

  Input Parameters:
+ ds - The `PetscDS`
. f  - The test field number
- g  - The field number

  Output Parameters:
+ g0 - integrand for the test and basis function term, see `PetscPointJacFn`
. g1 - integrand for the test function and basis function gradient term, see `PetscPointJacFn`
. g2 - integrand for the test function gradient and basis function term, see `PetscPointJacFn`
- g3 - integrand for the test function gradient and basis function gradient term, see `PetscPointJacFn`

  Level: intermediate

  Note:
  We are using a first order FEM model for the weak form\:

  $$
  \int_\Omega \phi\, g_0(u, u_t, \nabla u, x, t) \psi + \phi\, {\vec g}_1(u, u_t, \nabla u, x, t) \nabla \psi
  + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \nabla \psi
  $$

  Developer Note:
  The name is confusing since the function computes a matrix used to construct the preconditioner, not a preconditioner.

.seealso: `PetscDS`, `PetscDSSetJacobianPreconditioner()`, `PetscDSGetJacobian()`, `PetscPointJacFn`
@*/
PetscErrorCode PetscDSGetJacobianPreconditioner(PetscDS ds, PetscInt f, PetscInt g, PetscPointJacFn **g0, PetscPointJacFn **g1, PetscPointJacFn **g2, PetscPointJacFn **g3)
{
  PetscPointJacFn **tmp0, **tmp1, **tmp2, **tmp3;
  PetscInt          n0, n1, n2, n3;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
  PetscCheck(!(g < 0) && !(g >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", g, ds->Nf);
  PetscCall(PetscWeakFormGetJacobianPreconditioner(ds->wf, NULL, 0, f, g, 0, &n0, &tmp0, &n1, &tmp1, &n2, &tmp2, &n3, &tmp3));
  *g0 = tmp0 ? tmp0[0] : NULL;
  *g1 = tmp1 ? tmp1[0] : NULL;
  *g2 = tmp2 ? tmp2[0] : NULL;
  *g3 = tmp3 ? tmp3[0] : NULL;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  PetscDSSetJacobianPreconditioner - Set the pointwise Jacobian function for given test and basis fields that constructs the matrix used
  to compute the preconditioner. If this is missing, the system matrix is used to build the preconditioner.

  Not Collective

  Input Parameters:
+ ds - The `PetscDS`
. f  - The test field number
. g  - The field number
. g0 - integrand for the test and basis function term, see `PetscPointJacFn`
. g1 - integrand for the test function and basis function gradient term, see `PetscPointJacFn`
. g2 - integrand for the test function gradient and basis function term, see `PetscPointJacFn`
- g3 - integrand for the test function gradient and basis function gradient term, see `PetscPointJacFn`

  Level: intermediate

  Note:
  We are using a first order FEM model for the weak form\:

  $$
  \int_\Omega \phi\, g_0(u, u_t, \nabla u, x, t) \psi + \phi\, {\vec g}_1(u, u_t, \nabla u, x, t) \nabla \psi
  + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \nabla \psi
  $$

  Developer Note:
  The name is confusing since the function computes a matrix used to construct the preconditioner, not a preconditioner.

.seealso: `PetscDS`, `PetscDSGetJacobianPreconditioner()`, `PetscDSSetJacobian()`, `PetscPointJacFn`
@*/
PetscErrorCode PetscDSSetJacobianPreconditioner(PetscDS ds, PetscInt f, PetscInt g, PetscPointJacFn *g0, PetscPointJacFn *g1, PetscPointJacFn *g2, PetscPointJacFn *g3)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  if (g0) PetscValidFunction(g0, 4);
  if (g1) PetscValidFunction(g1, 5);
  if (g2) PetscValidFunction(g2, 6);
  if (g3) PetscValidFunction(g3, 7);
  PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
  PetscCheck(g >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", g);
  PetscCall(PetscWeakFormSetIndexJacobianPreconditioner(ds->wf, NULL, 0, f, g, 0, 0, g0, 0, g1, 0, g2, 0, g3));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSHasDynamicJacobian - Signals that a dynamic Jacobian, $dF/du_t$, has been set

  Not Collective

  Input Parameter:
. ds - The `PetscDS`

  Output Parameter:
. hasDynJac - flag that pointwise function for dynamic Jacobian has been set

  Level: intermediate

.seealso: `PetscDS`, `PetscDSGetDynamicJacobian()`, `PetscDSSetDynamicJacobian()`, `PetscDSGetJacobian()`
@*/
PetscErrorCode PetscDSHasDynamicJacobian(PetscDS ds, PetscBool *hasDynJac)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  PetscCall(PetscWeakFormHasDynamicJacobian(ds->wf, hasDynJac));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  PetscDSGetDynamicJacobian - Get the pointwise dynamic Jacobian, $dF/du_t$, function for given test and basis field

  Not Collective

  Input Parameters:
+ ds - The `PetscDS`
. f  - The test field number
- g  - The field number

  Output Parameters:
+ g0 - integrand for the test and basis function term, see `PetscPointJacFn`
. g1 - integrand for the test function and basis function gradient term, see `PetscPointJacFn`
. g2 - integrand for the test function gradient and basis function term, see `PetscPointJacFn`
- g3 - integrand for the test function gradient and basis function gradient term, see `PetscPointJacFn`

  Level: intermediate

  Note:
  We are using a first order FEM model for the weak form\:

  $$
  \int_\Omega \phi\, g_0(u, u_t, \nabla u, x, t) \psi + \phi\, {\vec g}_1(u, u_t, \nabla u, x, t) \nabla \psi
  + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \nabla \psi
  $$

.seealso: `PetscDS`, `PetscDSSetJacobian()`, `PetscDSSetDynamicJacobian()`, `PetscPointJacFn`
@*/
PetscErrorCode PetscDSGetDynamicJacobian(PetscDS ds, PetscInt f, PetscInt g, PetscPointJacFn **g0, PetscPointJacFn **g1, PetscPointJacFn **g2, PetscPointJacFn **g3)
{
  PetscPointJacFn **tmp0, **tmp1, **tmp2, **tmp3;
  PetscInt          n0, n1, n2, n3;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
  PetscCheck(!(g < 0) && !(g >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", g, ds->Nf);
  PetscCall(PetscWeakFormGetDynamicJacobian(ds->wf, NULL, 0, f, g, 0, &n0, &tmp0, &n1, &tmp1, &n2, &tmp2, &n3, &tmp3));
  *g0 = tmp0 ? tmp0[0] : NULL;
  *g1 = tmp1 ? tmp1[0] : NULL;
  *g2 = tmp2 ? tmp2[0] : NULL;
  *g3 = tmp3 ? tmp3[0] : NULL;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  PetscDSSetDynamicJacobian - Set the pointwise dynamic Jacobian, $dF/du_t$, function for given test and basis fields

  Not Collective

  Input Parameters:
+ ds - The `PetscDS`
. f  - The test field number
. g  - The field number
. g0 - integrand for the test and basis function term, see `PetscPointJacFn`
. g1 - integrand for the test function and basis function gradient term, see `PetscPointJacFn`
. g2 - integrand for the test function gradient and basis function term, see `PetscPointJacFn`
- g3 - integrand for the test function gradient and basis function gradient term, see `PetscPointJacFn`

  Level: intermediate

  Note:
  We are using a first order FEM model for the weak form\:

  $$
  \int_\Omega \phi\, g_0(u, u_t, \nabla u, x, t) \psi + \phi\, {\vec g}_1(u, u_t, \nabla u, x, t) \nabla \psi
  + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \nabla \psi
  $$

.seealso: `PetscDS`, `PetscDSGetDynamicJacobian()`, `PetscDSGetJacobian()`, `PetscPointJacFn`
@*/
PetscErrorCode PetscDSSetDynamicJacobian(PetscDS ds, PetscInt f, PetscInt g, PetscPointJacFn *g0, PetscPointJacFn *g1, PetscPointJacFn *g2, PetscPointJacFn *g3)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  if (g0) PetscValidFunction(g0, 4);
  if (g1) PetscValidFunction(g1, 5);
  if (g2) PetscValidFunction(g2, 6);
  if (g3) PetscValidFunction(g3, 7);
  PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
  PetscCheck(g >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", g);
  PetscCall(PetscWeakFormSetIndexDynamicJacobian(ds->wf, NULL, 0, f, g, 0, 0, g0, 0, g1, 0, g2, 0, g3));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  PetscDSGetRiemannSolver - Returns the Riemann solver for the given field

  Not Collective

  Input Parameters:
+ ds - The `PetscDS` object
- f  - The field number

  Output Parameter:
. r - Riemann solver, see `PetscRiemannFn`

  Level: intermediate

.seealso: `PetscDS`, `PetscRiemannFn`, `PetscDSSetRiemannSolver()`
@*/
PetscErrorCode PetscDSGetRiemannSolver(PetscDS ds, PetscInt f, PetscRiemannFn **r)
{
  PetscRiemannFn **tmp;
  PetscInt         n;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  PetscAssertPointer(r, 3);
  PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
  PetscCall(PetscWeakFormGetRiemannSolver(ds->wf, NULL, 0, f, 0, &n, &tmp));
  *r = tmp ? tmp[0] : NULL;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  PetscDSSetRiemannSolver - Sets the Riemann solver for the given field

  Not Collective

  Input Parameters:
+ ds - The `PetscDS` object
. f  - The field number
- r  - Riemann solver, see `PetscRiemannFn`

  Level: intermediate

.seealso: `PetscDS`, `PetscRiemannFn`, `PetscDSGetRiemannSolver()`
@*/
PetscErrorCode PetscDSSetRiemannSolver(PetscDS ds, PetscInt f, PetscRiemannFn *r)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  if (r) PetscValidFunction(r, 3);
  PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
  PetscCall(PetscWeakFormSetIndexRiemannSolver(ds->wf, NULL, 0, f, 0, 0, r));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  PetscDSGetUpdate - Get the pointwise update function for a given field

  Not Collective

  Input Parameters:
+ ds - The `PetscDS`
- f  - The field number

  Output Parameter:
. update - update function, see `PetscPointFn`

  Level: intermediate

.seealso: `PetscDS`, `PetscPointFn`, `PetscDSSetUpdate()`, `PetscDSSetResidual()`
@*/
PetscErrorCode PetscDSGetUpdate(PetscDS ds, PetscInt f, PetscPointFn **update)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
  if (update) {
    PetscAssertPointer(update, 3);
    *update = ds->update[f];
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  PetscDSSetUpdate - Set the pointwise update function for a given field

  Not Collective

  Input Parameters:
+ ds     - The `PetscDS`
. f      - The field number
- update - update function, see `PetscPointFn`

  Level: intermediate

.seealso: `PetscDS`, `PetscPointFn`, `PetscDSGetResidual()`
@*/
PetscErrorCode PetscDSSetUpdate(PetscDS ds, PetscInt f, PetscPointFn *update)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  if (update) PetscValidFunction(update, 3);
  PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
  PetscCall(PetscDSEnlarge_Static(ds, f + 1));
  ds->update[f] = update;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  PetscDSGetContext - Returns the context that was passed by `PetscDSSetContext()`

  Not Collective

  Input Parameters:
+ ds  - The `PetscDS`
. f   - The field number
- ctx - the context

  Level: intermediate

  Fortran Notes:
  This only works when the context is a Fortran derived type or a `PetscObject`. Define `ctx` with
.vb
  type(tUsertype), pointer :: ctx
.ve

.seealso: `PetscDS`, `PetscPointFn`, `PetscDSSetContext()`
@*/
PetscErrorCode PetscDSGetContext(PetscDS ds, PetscInt f, PetscCtxRt ctx)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
  PetscAssertPointer(ctx, 3);
  *(void **)ctx = ds->ctx[f];
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  PetscDSSetContext - Sets the context that is passed back to some of the pointwise function callbacks used by this `PetscDS`

  Not Collective

  Input Parameters:
+ ds  - The `PetscDS`
. f   - The field number
- ctx - the context

  Level: intermediate

.seealso: `PetscDS`, `PetscPointFn`, `PetscDSGetContext()`
@*/
PetscErrorCode PetscDSSetContext(PetscDS ds, PetscInt f, PetscCtx ctx)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
  PetscCall(PetscDSEnlarge_Static(ds, f + 1));
  ds->ctx[f] = ctx;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  PetscDSGetBdResidual - Get the pointwise boundary residual function for a given test field

  Not Collective

  Input Parameters:
+ ds - The PetscDS
- f  - The test field number

  Output Parameters:
+ f0 - boundary integrand for the test function term, see `PetscBdPointFn`
- f1 - boundary integrand for the test function gradient term, see `PetscBdPointFn`

  Level: intermediate

  Note:
  We are using a first order FEM model for the weak form\:

  $$
  \int_\Gamma \phi {\vec f}_0(u, u_t, \nabla u, x, t) \cdot \hat n + \nabla\phi \cdot {\overleftrightarrow f}_1(u, u_t, \nabla u, x, t) \cdot \hat n
  $$

.seealso: `PetscDS`, `PetscBdPointFn`, `PetscDSSetBdResidual()`
@*/
PetscErrorCode PetscDSGetBdResidual(PetscDS ds, PetscInt f, PetscBdPointFn **f0, PetscBdPointFn **f1)
{
  PetscBdPointFn **tmp0, **tmp1;
  PetscInt         n0, n1;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
  PetscCall(PetscWeakFormGetBdResidual(ds->wf, NULL, 0, f, 0, &n0, &tmp0, &n1, &tmp1));
  *f0 = tmp0 ? tmp0[0] : NULL;
  *f1 = tmp1 ? tmp1[0] : NULL;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  PetscDSSetBdResidual - Get the pointwise boundary residual function for a given test field

  Not Collective

  Input Parameters:
+ ds - The `PetscDS`
. f  - The test field number
. f0 - boundary integrand for the test function term, see `PetscBdPointFn`
- f1 - boundary integrand for the test function gradient term, see `PetscBdPointFn`

  Level: intermediate

  Note:
  We are using a first order FEM model for the weak form\:

  $$
  \int_\Gamma \phi {\vec f}_0(u, u_t, \nabla u, x, t) \cdot \hat n + \nabla\phi \cdot {\overleftrightarrow f}_1(u, u_t, \nabla u, x, t) \cdot \hat n
  $$

.seealso: `PetscDS`, `PetscBdPointFn`, `PetscDSGetBdResidual()`
@*/
PetscErrorCode PetscDSSetBdResidual(PetscDS ds, PetscInt f, PetscBdPointFn *f0, PetscBdPointFn *f1)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
  PetscCall(PetscWeakFormSetIndexBdResidual(ds->wf, NULL, 0, f, 0, 0, f0, 0, f1));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSHasBdJacobian - Indicates that boundary Jacobian functions have been set

  Not Collective

  Input Parameter:
. ds - The `PetscDS`

  Output Parameter:
. hasBdJac - flag that pointwise function for the boundary Jacobian has been set

  Level: intermediate

.seealso: `PetscDS`, `PetscDSHasJacobian()`, `PetscDSSetBdJacobian()`, `PetscDSGetBdJacobian()`
@*/
PetscErrorCode PetscDSHasBdJacobian(PetscDS ds, PetscBool *hasBdJac)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  PetscAssertPointer(hasBdJac, 2);
  PetscCall(PetscWeakFormHasBdJacobian(ds->wf, hasBdJac));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  PetscDSGetBdJacobian - Get the pointwise boundary Jacobian function for given test and basis field

  Not Collective

  Input Parameters:
+ ds - The `PetscDS`
. f  - The test field number
- g  - The field number

  Output Parameters:
+ g0 - integrand for the test and basis function term, see `PetscBdPointJacFn`
. g1 - integrand for the test function and basis function gradient term, see `PetscBdPointJacFn`
. g2 - integrand for the test function gradient and basis function term, see `PetscBdPointJacFn`
- g3 - integrand for the test function gradient and basis function gradient term, see `PetscBdPointJacFn`

  Level: intermediate

  Note:
  We are using a first order FEM model for the weak form\:

  $$
  \int_\Gamma \phi\, {\vec g}_0(u, u_t, \nabla u, x, t) \cdot \hat n \psi + \phi\, {\vec g}_1(u, u_t, \nabla u, x, t) \cdot \hat n \nabla \psi
  + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \cdot \hat n \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \hat n \cdot \nabla \psi
  $$

.seealso: `PetscDS`, `PetscBdPointJacFn`, `PetscDSSetBdJacobian()`
@*/
PetscErrorCode PetscDSGetBdJacobian(PetscDS ds, PetscInt f, PetscInt g, PetscBdPointJacFn **g0, PetscBdPointJacFn **g1, PetscBdPointJacFn **g2, PetscBdPointJacFn **g3)
{
  PetscBdPointJacFn **tmp0, **tmp1, **tmp2, **tmp3;
  PetscInt            n0, n1, n2, n3;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
  PetscCheck(!(g < 0) && !(g >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", g, ds->Nf);
  PetscCall(PetscWeakFormGetBdJacobian(ds->wf, NULL, 0, f, g, 0, &n0, &tmp0, &n1, &tmp1, &n2, &tmp2, &n3, &tmp3));
  *g0 = tmp0 ? tmp0[0] : NULL;
  *g1 = tmp1 ? tmp1[0] : NULL;
  *g2 = tmp2 ? tmp2[0] : NULL;
  *g3 = tmp3 ? tmp3[0] : NULL;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  PetscDSSetBdJacobian - Set the pointwise boundary Jacobian function for given test and basis field

  Not Collective

  Input Parameters:
+ ds - The PetscDS
. f  - The test field number
. g  - The field number
. g0 - integrand for the test and basis function term, see `PetscBdPointJacFn`
. g1 - integrand for the test function and basis function gradient term, see `PetscBdPointJacFn`
. g2 - integrand for the test function gradient and basis function term, see `PetscBdPointJacFn`
- g3 - integrand for the test function gradient and basis function gradient term, see `PetscBdPointJacFn`

  Level: intermediate

  Note:
  We are using a first order FEM model for the weak form\:

  $$
  \int_\Gamma \phi\, {\vec g}_0(u, u_t, \nabla u, x, t) \cdot \hat n \psi + \phi\, {\vec g}_1(u, u_t, \nabla u, x, t) \cdot \hat n \nabla \psi
  + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \cdot \hat n \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \hat n \cdot \nabla \psi
  $$

.seealso: `PetscDS`, `PetscBdPointJacFn`, `PetscDSGetBdJacobian()`
@*/
PetscErrorCode PetscDSSetBdJacobian(PetscDS ds, PetscInt f, PetscInt g, PetscBdPointJacFn *g0, PetscBdPointJacFn *g1, PetscBdPointJacFn *g2, PetscBdPointJacFn *g3)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  if (g0) PetscValidFunction(g0, 4);
  if (g1) PetscValidFunction(g1, 5);
  if (g2) PetscValidFunction(g2, 6);
  if (g3) PetscValidFunction(g3, 7);
  PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
  PetscCheck(g >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", g);
  PetscCall(PetscWeakFormSetIndexBdJacobian(ds->wf, NULL, 0, f, g, 0, 0, g0, 0, g1, 0, g2, 0, g3));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSHasBdJacobianPreconditioner - Signals that boundary Jacobian preconditioner functions have been set with `PetscDSSetBdJacobianPreconditioner()`

  Not Collective

  Input Parameter:
. ds - The `PetscDS`

  Output Parameter:
. hasBdJacPre - flag that pointwise function for the boundary Jacobian matrix to construct the preconditioner has been set

  Level: intermediate

  Developer Note:
  The name is confusing since the function computes a matrix used to construct the preconditioner, not a preconditioner.

.seealso: `PetscDS`, `PetscDSHasJacobian()`, `PetscDSSetBdJacobian()`, `PetscDSGetBdJacobian()`
@*/
PetscErrorCode PetscDSHasBdJacobianPreconditioner(PetscDS ds, PetscBool *hasBdJacPre)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  PetscAssertPointer(hasBdJacPre, 2);
  PetscCall(PetscWeakFormHasBdJacobianPreconditioner(ds->wf, hasBdJacPre));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  PetscDSGetBdJacobianPreconditioner - Get the pointwise boundary Jacobian function for given test and basis field that constructs the
  matrix used to construct the preconditioner

  Not Collective; No Fortran Support

  Input Parameters:
+ ds - The `PetscDS`
. f  - The test field number
- g  - The field number

  Output Parameters:
+ g0 - integrand for the test and basis function term, see `PetscBdPointJacFn`
. g1 - integrand for the test function and basis function gradient term, see `PetscBdPointJacFn`
. g2 - integrand for the test function gradient and basis function term, see `PetscBdPointJacFn`
- g3 - integrand for the test function gradient and basis function gradient term, see `PetscBdPointJacFn`

  Level: intermediate

  Note:
  We are using a first order FEM model for the weak form\:

  $$
  \int_\Gamma \phi\, {\vec g}_0(u, u_t, \nabla u, x, t) \cdot \hat n \psi + \phi\, {\vec g}_1(u, u_t, \nabla u, x, t) \cdot \hat n \nabla \psi
  + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \cdot \hat n \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \hat n \cdot \nabla \psi
  $$

  Developer Note:
  The name is confusing since the function computes a matrix used to construct the preconditioner, not a preconditioner.

.seealso: `PetscDS`, `PetscBdPointJacFn`, `PetscDSSetBdJacobianPreconditioner()`
@*/
PetscErrorCode PetscDSGetBdJacobianPreconditioner(PetscDS ds, PetscInt f, PetscInt g, PetscBdPointJacFn **g0, PetscBdPointJacFn **g1, PetscBdPointJacFn **g2, PetscBdPointJacFn **g3)
{
  PetscBdPointJacFn **tmp0, **tmp1, **tmp2, **tmp3;
  PetscInt            n0, n1, n2, n3;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
  PetscCheck(!(g < 0) && !(g >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", g, ds->Nf);
  PetscCall(PetscWeakFormGetBdJacobianPreconditioner(ds->wf, NULL, 0, f, g, 0, &n0, &tmp0, &n1, &tmp1, &n2, &tmp2, &n3, &tmp3));
  *g0 = tmp0 ? tmp0[0] : NULL;
  *g1 = tmp1 ? tmp1[0] : NULL;
  *g2 = tmp2 ? tmp2[0] : NULL;
  *g3 = tmp3 ? tmp3[0] : NULL;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  PetscDSSetBdJacobianPreconditioner - Set the pointwise boundary Jacobian preconditioner function for given test and basis field that constructs the
  matrix used to construct the preconditioner

  Not Collective; No Fortran Support

  Input Parameters:
+ ds - The `PetscDS`
. f  - The test field number
. g  - The field number
. g0 - integrand for the test and basis function term, see `PetscBdPointJacFn`
. g1 - integrand for the test function and basis function gradient term, see `PetscBdPointJacFn`
. g2 - integrand for the test function gradient and basis function term, see `PetscBdPointJacFn`
- g3 - integrand for the test function gradient and basis function gradient term, see `PetscBdPointJacFn`

  Level: intermediate

  Note:
  We are using a first order FEM model for the weak form\:

  $$
  \int_\Gamma \phi\, {\vec g}_0(u, u_t, \nabla u, x, t) \cdot \hat n \psi + \phi\, {\vec g}_1(u, u_t, \nabla u, x, t) \cdot \hat n \nabla \psi
  + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \cdot \hat n \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \hat n \cdot \nabla \psi
  $$

  Developer Note:
  The name is confusing since the function computes a matrix used to construct the preconditioner, not a preconditioner.

.seealso: `PetscDS`, `PetscBdPointJacFn`, `PetscDSGetBdJacobianPreconditioner()`
@*/
PetscErrorCode PetscDSSetBdJacobianPreconditioner(PetscDS ds, PetscInt f, PetscInt g, PetscBdPointJacFn *g0, PetscBdPointJacFn *g1, PetscBdPointJacFn *g2, PetscBdPointJacFn *g3)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  if (g0) PetscValidFunction(g0, 4);
  if (g1) PetscValidFunction(g1, 5);
  if (g2) PetscValidFunction(g2, 6);
  if (g3) PetscValidFunction(g3, 7);
  PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
  PetscCheck(g >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", g);
  PetscCall(PetscWeakFormSetIndexBdJacobianPreconditioner(ds->wf, NULL, 0, f, g, 0, 0, g0, 0, g1, 0, g2, 0, g3));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  PetscDSGetExactSolution - Get the pointwise exact solution function for a given test field

  Not Collective

  Input Parameters:
+ prob - The `PetscDS`
- f    - The test field number

  Output Parameters:
+ sol - exact solution function for the test field, see `PetscPointExactSolutionFn`
- ctx - exact solution context

  Level: intermediate

.seealso: `PetscDS`, `PetscPointExactSolutionFn`, `PetscDSSetExactSolution()`, `PetscDSGetExactSolutionTimeDerivative()`
@*/
PetscErrorCode PetscDSGetExactSolution(PetscDS prob, PetscInt f, PetscPointExactSolutionFn **sol, void **ctx)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
  PetscCheck(!(f < 0) && !(f >= prob->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, prob->Nf);
  if (sol) {
    PetscAssertPointer(sol, 3);
    *sol = prob->exactSol[f];
  }
  if (ctx) {
    PetscAssertPointer(ctx, 4);
    *ctx = prob->exactCtx[f];
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  PetscDSSetExactSolution - Set the pointwise exact solution function for a given test field

  Not Collective

  Input Parameters:
+ prob - The `PetscDS`
. f    - The test field number
. sol  - solution function for the test fields, see `PetscPointExactSolutionFn`
- ctx  - solution context or `NULL`

  Level: intermediate

.seealso: `PetscDS`, `PetscPointExactSolutionFn`, `PetscDSGetExactSolution()`
@*/
PetscErrorCode PetscDSSetExactSolution(PetscDS prob, PetscInt f, PetscPointExactSolutionFn *sol, PetscCtx ctx)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
  PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
  PetscCall(PetscDSEnlarge_Static(prob, f + 1));
  if (sol) {
    PetscValidFunction(sol, 3);
    prob->exactSol[f] = sol;
  }
  if (ctx) {
    PetscValidFunction(ctx, 4);
    prob->exactCtx[f] = ctx;
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  PetscDSGetExactSolutionTimeDerivative - Get the pointwise time derivative of the exact solution function for a given test field

  Not Collective

  Input Parameters:
+ prob - The `PetscDS`
- f    - The test field number

  Output Parameters:
+ sol - time derivative of the exact solution for the test field, see `PetscPointExactSolutionFn`
- ctx - the exact solution context

  Level: intermediate

.seealso: `PetscDS`, `PetscPointExactSolutionFn`, `PetscDSSetExactSolutionTimeDerivative()`, `PetscDSGetExactSolution()`
@*/
PetscErrorCode PetscDSGetExactSolutionTimeDerivative(PetscDS prob, PetscInt f, PetscPointExactSolutionFn **sol, void **ctx)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
  PetscCheck(!(f < 0) && !(f >= prob->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, prob->Nf);
  if (sol) {
    PetscAssertPointer(sol, 3);
    *sol = prob->exactSol_t[f];
  }
  if (ctx) {
    PetscAssertPointer(ctx, 4);
    *ctx = prob->exactCtx_t[f];
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  PetscDSSetExactSolutionTimeDerivative - Set the pointwise time derivative of the exact solution function for a given test field

  Not Collective

  Input Parameters:
+ prob - The `PetscDS`
. f    - The test field number
. sol  - time derivative of the solution function for the test fields, see `PetscPointExactSolutionFn`
- ctx  - the solution context or `NULL`

  Level: intermediate

.seealso: `PetscDS`, `PetscPointExactSolutionFn`, `PetscDSGetExactSolutionTimeDerivative()`, `PetscDSSetExactSolution()`
@*/
PetscErrorCode PetscDSSetExactSolutionTimeDerivative(PetscDS prob, PetscInt f, PetscPointExactSolutionFn *sol, PetscCtx ctx)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
  PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
  PetscCall(PetscDSEnlarge_Static(prob, f + 1));
  if (sol) {
    PetscValidFunction(sol, 3);
    prob->exactSol_t[f] = sol;
  }
  if (ctx) {
    PetscValidFunction(ctx, 4);
    prob->exactCtx_t[f] = ctx;
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  PetscDSGetLowerBound - Get the pointwise lower bound function for a given field

  Not Collective

  Input Parameters:
+ ds - The PetscDS
- f  - The field number

  Output Parameters:
+ lb  - lower bound function for the field, see `PetscPointBoundFn`
- ctx - lower bound context that was set with `PetscDSSetLowerBound()`

  Level: intermediate

.seealso: `PetscDS`, `PetscPointBoundFn`, `PetscDSSetLowerBound()`, `PetscDSGetUpperBound()`, `PetscDSGetExactSolution()`
@*/
PetscErrorCode PetscDSGetLowerBound(PetscDS ds, PetscInt f, PetscPointBoundFn **lb, void **ctx)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
  if (lb) {
    PetscAssertPointer(lb, 3);
    *lb = ds->lowerBound[f];
  }
  if (ctx) {
    PetscAssertPointer(ctx, 4);
    *ctx = ds->lowerCtx[f];
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  PetscDSSetLowerBound - Set the pointwise lower bound function for a given field

  Not Collective

  Input Parameters:
+ ds  - The `PetscDS`
. f   - The field number
. lb  - lower bound function for the test fields, see `PetscPointBoundFn`
- ctx - lower bound context or `NULL` which will be passed to `lb`

  Level: intermediate

.seealso: `PetscDS`, `PetscPointBoundFn`, `PetscDSGetLowerBound()`, `PetscDSGetUpperBound()`, `PetscDSGetExactSolution()`
@*/
PetscErrorCode PetscDSSetLowerBound(PetscDS ds, PetscInt f, PetscPointBoundFn *lb, PetscCtx ctx)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
  PetscCall(PetscDSEnlarge_Static(ds, f + 1));
  if (lb) {
    PetscValidFunction(lb, 3);
    ds->lowerBound[f] = lb;
  }
  if (ctx) {
    PetscValidFunction(ctx, 4);
    ds->lowerCtx[f] = ctx;
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  PetscDSGetUpperBound - Get the pointwise upper bound function for a given field

  Not Collective

  Input Parameters:
+ ds - The `PetscDS`
- f  - The field number

  Output Parameters:
+ ub  - upper bound function for the field, see `PetscPointBoundFn`
- ctx - upper bound context that was set with `PetscDSSetUpperBound()`

  Level: intermediate

.seealso: `PetscDS`, `PetscPointBoundFn`, `PetscDSSetUpperBound()`, `PetscDSGetLowerBound()`, `PetscDSGetExactSolution()`
@*/
PetscErrorCode PetscDSGetUpperBound(PetscDS ds, PetscInt f, PetscPointBoundFn **ub, void **ctx)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
  if (ub) {
    PetscAssertPointer(ub, 3);
    *ub = ds->upperBound[f];
  }
  if (ctx) {
    PetscAssertPointer(ctx, 4);
    *ctx = ds->upperCtx[f];
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  PetscDSSetUpperBound - Set the pointwise upper bound function for a given field

  Not Collective

  Input Parameters:
+ ds  - The `PetscDS`
. f   - The field number
. ub  - upper bound function for the test fields, see `PetscPointBoundFn`
- ctx - context or `NULL` that will be passed to `ub`

  Level: intermediate

.seealso: `PetscDS`, `PetscPointBoundFn`, `PetscDSGetUpperBound()`, `PetscDSGetLowerBound()`, `PetscDSGetExactSolution()`
@*/
PetscErrorCode PetscDSSetUpperBound(PetscDS ds, PetscInt f, PetscPointBoundFn *ub, PetscCtx ctx)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  PetscCheck(f >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be non-negative", f);
  PetscCall(PetscDSEnlarge_Static(ds, f + 1));
  if (ub) {
    PetscValidFunction(ub, 3);
    ds->upperBound[f] = ub;
  }
  if (ctx) {
    PetscValidFunction(ctx, 4);
    ds->upperCtx[f] = ctx;
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  PetscDSGetConstants - Returns the array of constants passed to point functions from a `PetscDS` object

  Not Collective

  Input Parameter:
. ds - The `PetscDS` object

  Output Parameters:
+ numConstants - The number of constants, or pass in `NULL` if not required
- constants    - The array of constants, `NULL` if there are none

  Level: intermediate

.seealso: `PetscDS`, `PetscDSSetConstants()`, `PetscDSCreate()`
@*/
PetscErrorCode PetscDSGetConstants(PetscDS ds, PeOp PetscInt *numConstants, PeOp const PetscScalar *constants[])
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  if (numConstants) {
    PetscAssertPointer(numConstants, 2);
    *numConstants = ds->numConstants;
  }
  if (constants) {
    PetscAssertPointer(constants, 3);
    *constants = ds->constants;
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  PetscDSSetConstants - Set the array of constants passed to point functions from a `PetscDS`

  Not Collective

  Input Parameters:
+ ds           - The `PetscDS` object
. numConstants - The number of constants
- constants    - The array of constants, `NULL` if there are none

  Level: intermediate

.seealso: `PetscDS`, `PetscDSGetConstants()`, `PetscDSCreate()`
@*/
PetscErrorCode PetscDSSetConstants(PetscDS ds, PetscInt numConstants, PetscScalar constants[])
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  if (numConstants != ds->numConstants) {
    PetscCall(PetscFree(ds->constants));
    ds->numConstants = numConstants;
    PetscCall(PetscMalloc1(ds->numConstants + ds->numFuncConstants, &ds->constants));
  }
  if (ds->numConstants) {
    PetscAssertPointer(constants, 3);
    PetscCall(PetscArraycpy(ds->constants, constants, ds->numConstants));
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  PetscDSSetIntegrationParameters - Set the parameters for a particular integration

  Not Collective

  Input Parameters:
+ ds     - The `PetscDS` object
. fieldI - The test field for a given point function, or `PETSC_DETERMINE`
- fieldJ - The basis field for a given point function, or `PETSC_DETERMINE`

  Level: intermediate

.seealso: `PetscDS`, `PetscDSSetConstants()`, `PetscDSGetConstants()`, `PetscDSCreate()`
@*/
PetscErrorCode PetscDSSetIntegrationParameters(PetscDS ds, PetscInt fieldI, PetscInt fieldJ)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  ds->constants[ds->numConstants]     = fieldI;
  ds->constants[ds->numConstants + 1] = fieldJ;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  PetscDSSetCellParameters - Set the parameters for a particular cell

  Not Collective

  Input Parameters:
+ ds     - The `PetscDS` object
- volume - The cell volume

  Level: intermediate

.seealso: `PetscDS`, `PetscDSSetConstants()`, `PetscDSGetConstants()`, `PetscDSCreate()`
@*/
PetscErrorCode PetscDSSetCellParameters(PetscDS ds, PetscReal volume)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  ds->constants[ds->numConstants + 2] = volume;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSGetFieldIndex - Returns the index of the given field

  Not Collective

  Input Parameters:
+ prob - The `PetscDS` object
- disc - The discretization object

  Output Parameter:
. f - The field number

  Level: beginner

.seealso: `PetscDS`, `PetscGetDiscretization()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
@*/
PetscErrorCode PetscDSGetFieldIndex(PetscDS prob, PetscObject disc, PetscInt *f)
{
  PetscInt g;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
  PetscAssertPointer(f, 3);
  *f = -1;
  for (g = 0; g < prob->Nf; ++g) {
    if (disc == prob->disc[g]) break;
  }
  PetscCheck(g != prob->Nf, PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_WRONG, "Field not found in PetscDS.");
  *f = g;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSGetFieldSize - Returns the size of the given field in the full space basis

  Not Collective

  Input Parameters:
+ prob - The `PetscDS` object
- f    - The field number

  Output Parameter:
. size - The size

  Level: beginner

.seealso: `PetscDS`, `PetscDSGetFieldOffset()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
@*/
PetscErrorCode PetscDSGetFieldSize(PetscDS prob, PetscInt f, PetscInt *size)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
  PetscAssertPointer(size, 3);
  PetscCheck(!(f < 0) && !(f >= prob->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, prob->Nf);
  PetscCall(PetscDSSetUp(prob));
  *size = prob->Nb[f];
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSGetFieldOffset - Returns the offset of the given field in the full space basis

  Not Collective

  Input Parameters:
+ prob - The `PetscDS` object
- f    - The field number

  Output Parameter:
. off - The offset

  Level: beginner

.seealso: `PetscDS`, `PetscDSGetFieldSize()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
@*/
PetscErrorCode PetscDSGetFieldOffset(PetscDS prob, PetscInt f, PetscInt *off)
{
  PetscInt size, g;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
  PetscAssertPointer(off, 3);
  PetscCheck(!(f < 0) && !(f >= prob->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, prob->Nf);
  *off = 0;
  for (g = 0; g < f; ++g) {
    PetscCall(PetscDSGetFieldSize(prob, g, &size));
    *off += size;
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSGetFieldOffsetCohesive - Returns the offset of the given field in the full space basis on a cohesive cell

  Not Collective

  Input Parameters:
+ ds - The `PetscDS` object
- f  - The field number

  Output Parameter:
. off - The offset

  Level: beginner

.seealso: `PetscDS`, `PetscDSGetFieldSize()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
@*/
PetscErrorCode PetscDSGetFieldOffsetCohesive(PetscDS ds, PetscInt f, PetscInt *off)
{
  PetscInt size, g;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  PetscAssertPointer(off, 3);
  PetscCheck(!(f < 0) && !(f >= ds->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, ds->Nf);
  *off = 0;
  for (g = 0; g < f; ++g) {
    PetscBool cohesive;

    PetscCall(PetscDSGetCohesive(ds, g, &cohesive));
    PetscCall(PetscDSGetFieldSize(ds, g, &size));
    *off += cohesive ? size : size * 2;
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSGetDimensions - Returns the size of the approximation space for each field on an evaluation point

  Not Collective

  Input Parameter:
. prob - The `PetscDS` object

  Output Parameter:
. dimensions - The number of dimensions

  Level: beginner

.seealso: `PetscDS`, `PetscDSGetComponentOffsets()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
@*/
PetscErrorCode PetscDSGetDimensions(PetscDS prob, PetscInt *dimensions[])
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
  PetscCall(PetscDSSetUp(prob));
  PetscAssertPointer(dimensions, 2);
  *dimensions = prob->Nb;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSGetComponents - Returns the number of components for each field on an evaluation point

  Not Collective

  Input Parameter:
. prob - The `PetscDS` object

  Output Parameter:
. components - The number of components

  Level: beginner

.seealso: `PetscDS`, `PetscDSGetComponentOffsets()`, `PetscDSGetNumFields()`, `PetscDSCreate()`
@*/
PetscErrorCode PetscDSGetComponents(PetscDS prob, PetscInt *components[])
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
  PetscCall(PetscDSSetUp(prob));
  PetscAssertPointer(components, 2);
  *components = prob->Nc;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSGetComponentOffset - Returns the offset of the given field on an evaluation point

  Not Collective

  Input Parameters:
+ prob - The `PetscDS` object
- f    - The field number

  Output Parameter:
. off - The offset

  Level: beginner

.seealso: `PetscDS`, `PetscDSGetNumFields()`, `PetscDSCreate()`
@*/
PetscErrorCode PetscDSGetComponentOffset(PetscDS prob, PetscInt f, PetscInt *off)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
  PetscAssertPointer(off, 3);
  PetscCheck(!(f < 0) && !(f >= prob->Nf), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, prob->Nf);
  PetscCall(PetscDSSetUp(prob));
  *off = prob->off[f];
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSGetComponentOffsets - Returns the offset of each field on an evaluation point

  Not Collective

  Input Parameter:
. prob - The `PetscDS` object

  Output Parameter:
. offsets - The offsets

  Level: beginner

.seealso: `PetscDS`, `PetscDSGetNumFields()`, `PetscDSCreate()`
@*/
PetscErrorCode PetscDSGetComponentOffsets(PetscDS prob, PetscInt *offsets[])
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
  PetscAssertPointer(offsets, 2);
  PetscCall(PetscDSSetUp(prob));
  *offsets = prob->off;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSGetComponentDerivativeOffsets - Returns the offset of each field derivative on an evaluation point

  Not Collective

  Input Parameter:
. prob - The `PetscDS` object

  Output Parameter:
. offsets - The offsets

  Level: beginner

.seealso: `PetscDS`, `PetscDSGetNumFields()`, `PetscDSCreate()`
@*/
PetscErrorCode PetscDSGetComponentDerivativeOffsets(PetscDS prob, PetscInt *offsets[])
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
  PetscAssertPointer(offsets, 2);
  PetscCall(PetscDSSetUp(prob));
  *offsets = prob->offDer;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSGetComponentOffsetsCohesive - Returns the offset of each field on an evaluation point

  Not Collective

  Input Parameters:
+ ds - The `PetscDS` object
- s  - The cohesive side, 0 for negative, 1 for positive, 2 for cohesive

  Output Parameter:
. offsets - The offsets

  Level: beginner

.seealso: `PetscDS`, `PetscDSGetNumFields()`, `PetscDSCreate()`
@*/
PetscErrorCode PetscDSGetComponentOffsetsCohesive(PetscDS ds, PetscInt s, PetscInt *offsets[])
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  PetscAssertPointer(offsets, 3);
  PetscCheck(ds->isCohesive, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Cohesive offsets are only valid for a cohesive DS");
  PetscCheck(!(s < 0) && !(s > 2), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cohesive side %" PetscInt_FMT " is not in [0, 2]", s);
  PetscCall(PetscDSSetUp(ds));
  *offsets = ds->offCohesive[s];
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSGetComponentDerivativeOffsetsCohesive - Returns the offset of each field derivative on an evaluation point

  Not Collective

  Input Parameters:
+ ds - The `PetscDS` object
- s  - The cohesive side, 0 for negative, 1 for positive, 2 for cohesive

  Output Parameter:
. offsets - The offsets

  Level: beginner

.seealso: `PetscDS`, `PetscDSGetNumFields()`, `PetscDSCreate()`
@*/
PetscErrorCode PetscDSGetComponentDerivativeOffsetsCohesive(PetscDS ds, PetscInt s, PetscInt *offsets[])
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  PetscAssertPointer(offsets, 3);
  PetscCheck(ds->isCohesive, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Cohesive offsets are only valid for a cohesive DS");
  PetscCheck(!(s < 0) && !(s > 2), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cohesive side %" PetscInt_FMT " is not in [0, 2]", s);
  PetscCall(PetscDSSetUp(ds));
  *offsets = ds->offDerCohesive[s];
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  PetscDSGetTabulation - Return the basis tabulation at quadrature points for the volume discretization

  Not Collective

  Input Parameter:
. prob - The `PetscDS` object

  Output Parameter:
. T - The basis function and derivatives tabulation at quadrature points for each field, see `PetscTabulation` for its details

  Level: intermediate

  Note:
  The tabulation is only valid so long as the `PetscDS` has not be destroyed. There is no `PetscDSRestoreTabulation()` in C.

  Fortran Note:
  Use the declaration
.vb
  PetscTabulation, pointer :: tab(:)
.ve
  and access the values using, for example,
.vb
  tab(i)%ptr%K
  tab(i)%ptr%T(j)%ptr
.ve
  where $ i = 1, 2, ..., Nf $ and $ j = 1, 2, ..., tab(i)%ptr%K+1 $.

  Use `PetscDSRestoreTabulation()` to restore the array

  Developer Note:
  The Fortran language syntax does not directly support arrays of pointers, the '%ptr' notation allows mimicking their use in Fortran.

.seealso: `PetscDS`, `PetscTabulation`, `PetscDSCreate()`
@*/
PetscErrorCode PetscDSGetTabulation(PetscDS prob, PetscTabulation *T[]) PeNS
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
  PetscAssertPointer(T, 2);
  PetscCall(PetscDSSetUp(prob));
  *T = prob->T;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  PetscDSGetFaceTabulation - Return the basis tabulation at quadrature points on the faces

  Not Collective

  Input Parameter:
. prob - The `PetscDS` object

  Output Parameter:
. Tf - The basis function and derivative tabulation on each local face at quadrature points for each field

  Level: intermediate

  Note:
  The tabulation is only valid so long as the `PetscDS` has not be destroyed. There is no `PetscDSRestoreFaceTabulation()` in C.

.seealso: `PetscTabulation`, `PetscDS`, `PetscDSGetTabulation()`, `PetscDSCreate()`
@*/
PetscErrorCode PetscDSGetFaceTabulation(PetscDS prob, PetscTabulation *Tf[])
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
  PetscAssertPointer(Tf, 2);
  PetscCall(PetscDSSetUp(prob));
  *Tf = prob->Tf;
  PetscFunctionReturn(PETSC_SUCCESS);
}

PetscErrorCode PetscDSGetEvaluationArrays(PetscDS prob, PetscScalar *u[], PetscScalar *u_t[], PetscScalar *u_x[])
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
  PetscCall(PetscDSSetUp(prob));
  if (u) {
    PetscAssertPointer(u, 2);
    *u = prob->u;
  }
  if (u_t) {
    PetscAssertPointer(u_t, 3);
    *u_t = prob->u_t;
  }
  if (u_x) {
    PetscAssertPointer(u_x, 4);
    *u_x = prob->u_x;
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

PetscErrorCode PetscDSGetWeakFormArrays(PetscDS prob, PetscScalar *f0[], PetscScalar *f1[], PetscScalar *g0[], PetscScalar *g1[], PetscScalar *g2[], PetscScalar *g3[])
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
  PetscCall(PetscDSSetUp(prob));
  if (f0) {
    PetscAssertPointer(f0, 2);
    *f0 = prob->f0;
  }
  if (f1) {
    PetscAssertPointer(f1, 3);
    *f1 = prob->f1;
  }
  if (g0) {
    PetscAssertPointer(g0, 4);
    *g0 = prob->g0;
  }
  if (g1) {
    PetscAssertPointer(g1, 5);
    *g1 = prob->g1;
  }
  if (g2) {
    PetscAssertPointer(g2, 6);
    *g2 = prob->g2;
  }
  if (g3) {
    PetscAssertPointer(g3, 7);
    *g3 = prob->g3;
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

PetscErrorCode PetscDSGetWorkspace(PetscDS prob, PetscReal **x, PetscScalar **basisReal, PetscScalar **basisDerReal, PetscScalar **testReal, PetscScalar **testDerReal)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
  PetscCall(PetscDSSetUp(prob));
  if (x) {
    PetscAssertPointer(x, 2);
    *x = prob->x;
  }
  if (basisReal) {
    PetscAssertPointer(basisReal, 3);
    *basisReal = prob->basisReal;
  }
  if (basisDerReal) {
    PetscAssertPointer(basisDerReal, 4);
    *basisDerReal = prob->basisDerReal;
  }
  if (testReal) {
    PetscAssertPointer(testReal, 5);
    *testReal = prob->testReal;
  }
  if (testDerReal) {
    PetscAssertPointer(testDerReal, 6);
    *testDerReal = prob->testDerReal;
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  PetscDSAddBoundary - Add a boundary condition to the model.

  Collective

  Input Parameters:
+ ds       - The `PetscDS` object
. type     - The type of condition, e.g. `DM_BC_ESSENTIAL`/`DM_BC_ESSENTIAL_FIELD` (Dirichlet), or `DM_BC_NATURAL` (Neumann)
. name     - The name for the boundary condition
. label    - The label defining constrained points
. Nv       - The number of `DMLabel` values for constrained points
. values   - An array of label values for constrained points
. field    - The field to constrain
. Nc       - The number of constrained field components (0 will constrain all fields)
. comps    - An array of constrained component numbers
. bcFunc   - A pointwise function giving boundary values
. bcFunc_t - A pointwise function giving the time derivative of the boundary values, or `NULL`
- ctx      - An optional user context for `bcFunc`

  Output Parameter:
. bd - The boundary number

  Options Database Keys:
+ -bc_<boundary name> <num>      - Overrides the boundary ids
- -bc_<boundary name>_comp <num> - Overrides the boundary components

  Level: developer

  Note:
  Both `bcFunc` and `bcFunc_t` will depend on the boundary condition type. If the type if `DM_BC_ESSENTIAL`, then the calling sequence is\:
.vb
  void bcFunc(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar bcval[])
.ve

  If the type is `DM_BC_ESSENTIAL_FIELD` or other _FIELD value, then the calling sequence is\:
.vb
  void bcFunc(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 time, const PetscReal x[], PetscScalar bcval[])
.ve
+ dim          - the coordinate dimension
. Nf           - the number of fields
. uOff         - the offset into u[] and u_t[] for each field
. uOff_x       - the offset into u_x[] for each field
. u            - each field evaluated at the current point
. u_t          - the time derivative of each field evaluated at the current point
. u_x          - the gradient of each field evaluated at the current point
. aOff         - the offset into a[] and a_t[] for each auxiliary field
. aOff_x       - the offset into a_x[] for each auxiliary field
. a            - each auxiliary field evaluated at the current point
. a_t          - the time derivative of each auxiliary field evaluated at the current point
. a_x          - the gradient of auxiliary each field evaluated at the current point
. t            - current time
. x            - coordinates of the current point
. numConstants - number of constant parameters
. constants    - constant parameters
- bcval        - output values at the current point

  Notes:
  The pointwise functions are used to provide boundary values for essential boundary
  conditions. In FEM, they are acting upon by dual basis functionals to generate FEM
  coefficients which are fixed. Natural boundary conditions signal to PETSc that boundary
  integrals should be performed, using the kernels from `PetscDSSetBdResidual()`.

.seealso: `PetscDS`, `PetscWeakForm`, `DMLabel`, `DMBoundaryConditionType`, `PetscDSAddBoundaryByName()`, `PetscDSGetBoundary()`, `PetscDSSetResidual()`, `PetscDSSetBdResidual()`
@*/
PetscErrorCode PetscDSAddBoundary(PetscDS ds, DMBoundaryConditionType type, const char name[], DMLabel label, PetscInt Nv, const PetscInt values[], PetscInt field, PetscInt Nc, const PetscInt comps[], PetscVoidFn *bcFunc, PetscVoidFn *bcFunc_t, PetscCtx ctx, PetscInt *bd)
{
  DSBoundary  head = ds->boundary, b;
  PetscInt    n    = 0;
  const char *lname;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  PetscValidLogicalCollectiveEnum(ds, type, 2);
  PetscAssertPointer(name, 3);
  PetscValidHeaderSpecific(label, DMLABEL_CLASSID, 4);
  PetscValidLogicalCollectiveInt(ds, Nv, 5);
  PetscValidLogicalCollectiveInt(ds, field, 7);
  PetscValidLogicalCollectiveInt(ds, Nc, 8);
  PetscCheck(field >= 0 && field < ds->Nf, PetscObjectComm((PetscObject)ds), PETSC_ERR_ARG_OUTOFRANGE, "Field %" PetscInt_FMT " is not in [0, %" PetscInt_FMT ")", field, ds->Nf);
  if (Nc > 0) {
    PetscInt *fcomps;
    PetscInt  c;

    PetscCall(PetscDSGetComponents(ds, &fcomps));
    PetscCheck(Nc <= fcomps[field], PetscObjectComm((PetscObject)ds), PETSC_ERR_ARG_OUTOFRANGE, "Number of constrained components %" PetscInt_FMT " > %" PetscInt_FMT " components for field %" PetscInt_FMT, Nc, fcomps[field], field);
    for (c = 0; c < Nc; ++c) {
      PetscCheck(comps[c] >= 0 && comps[c] < fcomps[field], PetscObjectComm((PetscObject)ds), PETSC_ERR_ARG_OUTOFRANGE, "Constrained component[%" PetscInt_FMT "] %" PetscInt_FMT " not in [0, %" PetscInt_FMT ") components for field %" PetscInt_FMT, c, comps[c], fcomps[field], field);
    }
  }
  PetscCall(PetscNew(&b));
  PetscCall(PetscStrallocpy(name, (char **)&b->name));
  PetscCall(PetscWeakFormCreate(PETSC_COMM_SELF, &b->wf));
  PetscCall(PetscWeakFormSetNumFields(b->wf, ds->Nf));
  PetscCall(PetscMalloc1(Nv, &b->values));
  if (Nv) PetscCall(PetscArraycpy(b->values, values, Nv));
  PetscCall(PetscMalloc1(Nc, &b->comps));
  if (Nc) PetscCall(PetscArraycpy(b->comps, comps, Nc));
  PetscCall(PetscObjectGetName((PetscObject)label, &lname));
  PetscCall(PetscStrallocpy(lname, (char **)&b->lname));
  b->type   = type;
  b->label  = label;
  b->Nv     = Nv;
  b->field  = field;
  b->Nc     = Nc;
  b->func   = bcFunc;
  b->func_t = bcFunc_t;
  b->ctx    = ctx;
  b->next   = NULL;
  /* Append to linked list so that we can preserve the order */
  if (!head) ds->boundary = b;
  while (head) {
    if (!head->next) {
      head->next = b;
      head       = b;
    }
    head = head->next;
    ++n;
  }
  if (bd) {
    PetscAssertPointer(bd, 13);
    *bd = n;
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

// PetscClangLinter pragma ignore: -fdoc-section-header-unknown
/*@C
  PetscDSAddBoundaryByName - Add a boundary condition to the model.

  Collective

  Input Parameters:
+ ds       - The `PetscDS` object
. type     - The type of condition, e.g. `DM_BC_ESSENTIAL`/`DM_BC_ESSENTIAL_FIELD` (Dirichlet), or `DM_BC_NATURAL` (Neumann)
. name     - The boundary condition name
. lname    - The name of the label defining constrained points
. Nv       - The number of `DMLabel` values for constrained points
. values   - An array of label values for constrained points
. field    - The field to constrain
. Nc       - The number of constrained field components (0 will constrain all fields)
. comps    - An array of constrained component numbers
. bcFunc   - A pointwise function giving boundary values
. bcFunc_t - A pointwise function giving the time derivative of the boundary values, or `NULL`
- ctx      - An optional user context for `bcFunc`

  Output Parameter:
. bd - The boundary number

  Options Database Keys:
+ -bc_<boundary name> <num>      - Overrides the boundary ids
- -bc_<boundary name>_comp <num> - Overrides the boundary components

  Calling Sequence of `bcFunc` and `bcFunc_t`:
  If the type is `DM_BC_ESSENTIAL`
.vb
  void bcFunc(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar bcval[])
.ve
  If the type is `DM_BC_ESSENTIAL_FIELD` or other _FIELD value,
.vb
  void bcFunc(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 time, const PetscReal x[], PetscScalar bcval[])
.ve
+ dim          - the coordinate dimension
. Nf           - the number of fields
. uOff         - the offset into `u`[] and `u_t`[] for each field
. uOff_x       - the offset into `u_x`[] for each field
. u            - each field evaluated at the current point
. u_t          - the time derivative of each field evaluated at the current point
. u_x          - the gradient of each field evaluated at the current point
. aOff         - the offset into `a`[] and `a_t`[] for each auxiliary field
. aOff_x       - the offset into `a_x`[] for each auxiliary field
. a            - each auxiliary field evaluated at the current point
. a_t          - the time derivative of each auxiliary field evaluated at the current point
. a_x          - the gradient of auxiliary each field evaluated at the current point
. t            - current time
. x            - coordinates of the current point
. numConstants - number of constant parameters
. constants    - constant parameters
- bcval        - output values at the current point

  Level: developer

  Notes:
  The pointwise functions are used to provide boundary values for essential boundary
  conditions. In FEM, they are acting upon by dual basis functionals to generate FEM
  coefficients which are fixed. Natural boundary conditions signal to PETSc that boundary
  integrals should be performed, using the kernels from `PetscDSSetBdResidual()`.

  This function should only be used with `DMFOREST` currently, since labels cannot be defined before the underlying `DMPLEX` is built.

.seealso: `PetscDS`, `PetscWeakForm`, `DMLabel`, `DMBoundaryConditionType`, `PetscDSAddBoundary()`, `PetscDSGetBoundary()`, `PetscDSSetResidual()`, `PetscDSSetBdResidual()`
@*/
PetscErrorCode PetscDSAddBoundaryByName(PetscDS ds, DMBoundaryConditionType type, const char name[], const char lname[], PetscInt Nv, const PetscInt values[], PetscInt field, PetscInt Nc, const PetscInt comps[], PetscVoidFn *bcFunc, PetscVoidFn *bcFunc_t, PetscCtx ctx, PetscInt *bd)
{
  DSBoundary head = ds->boundary, b;
  PetscInt   n    = 0;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  PetscValidLogicalCollectiveEnum(ds, type, 2);
  PetscAssertPointer(name, 3);
  PetscAssertPointer(lname, 4);
  PetscValidLogicalCollectiveInt(ds, Nv, 5);
  PetscValidLogicalCollectiveInt(ds, field, 7);
  PetscValidLogicalCollectiveInt(ds, Nc, 8);
  PetscCall(PetscNew(&b));
  PetscCall(PetscStrallocpy(name, (char **)&b->name));
  PetscCall(PetscWeakFormCreate(PETSC_COMM_SELF, &b->wf));
  PetscCall(PetscWeakFormSetNumFields(b->wf, ds->Nf));
  PetscCall(PetscMalloc1(Nv, &b->values));
  if (Nv) PetscCall(PetscArraycpy(b->values, values, Nv));
  PetscCall(PetscMalloc1(Nc, &b->comps));
  if (Nc) PetscCall(PetscArraycpy(b->comps, comps, Nc));
  PetscCall(PetscStrallocpy(lname, (char **)&b->lname));
  b->type   = type;
  b->label  = NULL;
  b->Nv     = Nv;
  b->field  = field;
  b->Nc     = Nc;
  b->func   = bcFunc;
  b->func_t = bcFunc_t;
  b->ctx    = ctx;
  b->next   = NULL;
  /* Append to linked list so that we can preserve the order */
  if (!head) ds->boundary = b;
  while (head) {
    if (!head->next) {
      head->next = b;
      head       = b;
    }
    head = head->next;
    ++n;
  }
  if (bd) {
    PetscAssertPointer(bd, 13);
    *bd = n;
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  PetscDSUpdateBoundary - Change a boundary condition for the model.

  Input Parameters:
+ ds       - The `PetscDS` object
. bd       - The boundary condition number
. type     - The type of condition, e.g. `DM_BC_ESSENTIAL`/`DM_BC_ESSENTIAL_FIELD` (Dirichlet), or `DM_BC_NATURAL` (Neumann)
. name     - The boundary condition name
. label    - The label defining constrained points
. Nv       - The number of `DMLabel` ids for constrained points
. values   - An array of ids for constrained points
. field    - The field to constrain
. Nc       - The number of constrained field components
. comps    - An array of constrained component numbers
. bcFunc   - A pointwise function giving boundary values
. bcFunc_t - A pointwise function giving the time derivative of the boundary values, or `NULL`
- ctx      - An optional user context for `bcFunc`

  Level: developer

  Notes:
  The pointwise functions are used to provide boundary values for essential boundary
  conditions. In FEM, they are acting upon by dual basis functionals to generate FEM
  coefficients which are fixed. Natural boundary conditions signal to PETSc that boundary
  integrals should be performed, using the kernels from `PetscDSSetBdResidual()`.

  The boundary condition number is the order in which it was registered. The user can get the number of boundary conditions from `PetscDSGetNumBoundary()`.
  See `PetscDSAddBoundary()` for a description of the calling sequences for the callbacks.

.seealso: `PetscDS`, `PetscWeakForm`, `DMBoundaryConditionType`, `PetscDSAddBoundary()`, `PetscDSGetBoundary()`, `PetscDSGetNumBoundary()`, `DMLabel`
@*/
PetscErrorCode PetscDSUpdateBoundary(PetscDS ds, PetscInt bd, DMBoundaryConditionType type, const char name[], DMLabel label, PetscInt Nv, const PetscInt values[], PetscInt field, PetscInt Nc, const PetscInt comps[], PetscVoidFn *bcFunc, PetscVoidFn *bcFunc_t, PetscCtx ctx)
{
  DSBoundary b = ds->boundary;
  PetscInt   n = 0;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  while (b) {
    if (n == bd) break;
    b = b->next;
    ++n;
  }
  PetscCheck(b, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Boundary %" PetscInt_FMT " is not in [0, %" PetscInt_FMT ")", bd, n);
  if (name) {
    PetscCall(PetscFree(b->name));
    PetscCall(PetscStrallocpy(name, (char **)&b->name));
  }
  b->type = type;
  if (label) {
    const char *name;

    b->label = label;
    PetscCall(PetscFree(b->lname));
    PetscCall(PetscObjectGetName((PetscObject)label, &name));
    PetscCall(PetscStrallocpy(name, (char **)&b->lname));
  }
  if (Nv >= 0) {
    b->Nv = Nv;
    PetscCall(PetscFree(b->values));
    PetscCall(PetscMalloc1(Nv, &b->values));
    if (Nv) PetscCall(PetscArraycpy(b->values, values, Nv));
  }
  if (field >= 0) b->field = field;
  if (Nc >= 0) {
    b->Nc = Nc;
    PetscCall(PetscFree(b->comps));
    PetscCall(PetscMalloc1(Nc, &b->comps));
    if (Nc) PetscCall(PetscArraycpy(b->comps, comps, Nc));
  }
  if (bcFunc) b->func = bcFunc;
  if (bcFunc_t) b->func_t = bcFunc_t;
  if (ctx) b->ctx = ctx;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSGetNumBoundary - Get the number of registered boundary conditions

  Input Parameter:
. ds - The `PetscDS` object

  Output Parameter:
. numBd - The number of boundary conditions

  Level: intermediate

.seealso: `PetscDS`, `PetscDSAddBoundary()`, `PetscDSGetBoundary()`
@*/
PetscErrorCode PetscDSGetNumBoundary(PetscDS ds, PetscInt *numBd)
{
  DSBoundary b = ds->boundary;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  PetscAssertPointer(numBd, 2);
  *numBd = 0;
  while (b) {
    ++(*numBd);
    b = b->next;
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  PetscDSGetBoundary - Gets a boundary condition from the model

  Input Parameters:
+ ds - The `PetscDS` object
- bd - The boundary condition number

  Output Parameters:
+ wf     - The `PetscWeakForm` holding the pointwise functions
. type   - The type of condition, e.g. `DM_BC_ESSENTIAL`/`DM_BC_ESSENTIAL_FIELD` (Dirichlet), or `DM_BC_NATURAL` (Neumann)
. name   - The boundary condition name
. label  - The label defining constrained points
. Nv     - The number of `DMLabel` ids for constrained points
. values - An array of ids for constrained points
. field  - The field to constrain
. Nc     - The number of constrained field components
. comps  - An array of constrained component numbers
. func   - A pointwise function giving boundary values
. func_t - A pointwise function giving the time derivative of the boundary values
- ctx    - An optional user context for `bcFunc`

  Options Database Keys:
+ -bc_<boundary name> <num>      - Overrides the boundary ids
- -bc_<boundary name>_comp <num> - Overrides the boundary components

  Level: developer

.seealso: `PetscDS`, `PetscWeakForm`, `DMBoundaryConditionType`, `PetscDSAddBoundary()`, `DMLabel`
@*/
PetscErrorCode PetscDSGetBoundary(PetscDS ds, PetscInt bd, PetscWeakForm *wf, DMBoundaryConditionType *type, const char *name[], DMLabel *label, PetscInt *Nv, const PetscInt *values[], PetscInt *field, PetscInt *Nc, const PetscInt *comps[], PetscVoidFn **func, PetscVoidFn **func_t, void **ctx)
{
  DSBoundary b = ds->boundary;
  PetscInt   n = 0;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  while (b) {
    if (n == bd) break;
    b = b->next;
    ++n;
  }
  PetscCheck(b, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Boundary %" PetscInt_FMT " is not in [0, %" PetscInt_FMT ")", bd, n);
  if (wf) {
    PetscAssertPointer(wf, 3);
    *wf = b->wf;
  }
  if (type) {
    PetscAssertPointer(type, 4);
    *type = b->type;
  }
  if (name) {
    PetscAssertPointer(name, 5);
    *name = b->name;
  }
  if (label) {
    PetscAssertPointer(label, 6);
    *label = b->label;
  }
  if (Nv) {
    PetscAssertPointer(Nv, 7);
    *Nv = b->Nv;
  }
  if (values) {
    PetscAssertPointer(values, 8);
    *values = b->values;
  }
  if (field) {
    PetscAssertPointer(field, 9);
    *field = b->field;
  }
  if (Nc) {
    PetscAssertPointer(Nc, 10);
    *Nc = b->Nc;
  }
  if (comps) {
    PetscAssertPointer(comps, 11);
    *comps = b->comps;
  }
  if (func) {
    PetscAssertPointer(func, 12);
    *func = b->func;
  }
  if (func_t) {
    PetscAssertPointer(func_t, 13);
    *func_t = b->func_t;
  }
  if (ctx) {
    PetscAssertPointer(ctx, 14);
    *ctx = b->ctx;
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSUpdateBoundaryLabels - Update `DMLabel` in each boundary condition using the label name and the input `DM`

  Not Collective

  Input Parameters:
+ ds - The source `PetscDS` object
- dm - The `DM` holding labels

  Level: intermediate

.seealso: `PetscDS`, `DMBoundary`, `DM`, `PetscDSCopyBoundary()`, `PetscDSCreate()`, `DMGetLabel()`
@*/
PetscErrorCode PetscDSUpdateBoundaryLabels(PetscDS ds, DM dm)
{
  DSBoundary b;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  PetscValidHeaderSpecific(dm, DM_CLASSID, 2);
  for (b = ds->boundary; b; b = b->next) {
    if (b->lname) PetscCall(DMGetLabel(dm, b->lname, &b->label));
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode DSBoundaryDuplicate_Internal(DSBoundary b, DSBoundary *bNew)
{
  PetscFunctionBegin;
  PetscCall(PetscNew(bNew));
  PetscCall(PetscWeakFormCreate(PETSC_COMM_SELF, &(*bNew)->wf));
  PetscCall(PetscWeakFormCopy(b->wf, (*bNew)->wf));
  PetscCall(PetscStrallocpy(b->name, (char **)&((*bNew)->name)));
  PetscCall(PetscStrallocpy(b->lname, (char **)&((*bNew)->lname)));
  (*bNew)->type  = b->type;
  (*bNew)->label = b->label;
  (*bNew)->Nv    = b->Nv;
  PetscCall(PetscMalloc1(b->Nv, &(*bNew)->values));
  PetscCall(PetscArraycpy((*bNew)->values, b->values, b->Nv));
  (*bNew)->field = b->field;
  (*bNew)->Nc    = b->Nc;
  PetscCall(PetscMalloc1(b->Nc, &(*bNew)->comps));
  PetscCall(PetscArraycpy((*bNew)->comps, b->comps, b->Nc));
  (*bNew)->func   = b->func;
  (*bNew)->func_t = b->func_t;
  (*bNew)->ctx    = b->ctx;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSCopyBoundary - Copy all boundary condition objects to the new `PetscDS`

  Not Collective

  Input Parameters:
+ ds        - The source `PetscDS` object
. numFields - The number of selected fields, or `PETSC_DEFAULT` for all fields
- fields    - The selected fields, or `NULL` for all fields

  Output Parameter:
. newds - The target `PetscDS`, now with a copy of the boundary conditions

  Level: intermediate

.seealso: `PetscDS`, `DMBoundary`, `PetscDSCopyEquations()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`, `PetscDSSetRiemannSolver()`, `PetscDSSetBdResidual()`, `PetscDSSetBdJacobian()`, `PetscDSCreate()`
@*/
PetscErrorCode PetscDSCopyBoundary(PetscDS ds, PetscInt numFields, const PetscInt fields[], PetscDS newds)
{
  DSBoundary b, *lastnext;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  PetscValidHeaderSpecific(newds, PETSCDS_CLASSID, 4);
  if (ds == newds) PetscFunctionReturn(PETSC_SUCCESS);
  PetscCall(PetscDSDestroyBoundary(newds));
  lastnext = &newds->boundary;
  for (b = ds->boundary; b; b = b->next) {
    DSBoundary bNew;
    PetscInt   fieldNew = -1;

    if (numFields > 0 && fields) {
      PetscInt f;

      for (f = 0; f < numFields; ++f)
        if (b->field == fields[f]) break;
      if (f == numFields) continue;
      fieldNew = f;
    }
    PetscCall(DSBoundaryDuplicate_Internal(b, &bNew));
    bNew->field = fieldNew < 0 ? b->field : fieldNew;
    *lastnext   = bNew;
    lastnext    = &bNew->next;
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSDestroyBoundary - Remove all `DMBoundary` objects from the `PetscDS`

  Not Collective

  Input Parameter:
. ds - The `PetscDS` object

  Level: intermediate

.seealso: `PetscDS`, `DMBoundary`, `PetscDSCopyBoundary()`, `PetscDSCopyEquations()`
@*/
PetscErrorCode PetscDSDestroyBoundary(PetscDS ds)
{
  DSBoundary next = ds->boundary;

  PetscFunctionBegin;
  while (next) {
    DSBoundary b = next;

    next = b->next;
    PetscCall(PetscWeakFormDestroy(&b->wf));
    PetscCall(PetscFree(b->name));
    PetscCall(PetscFree(b->lname));
    PetscCall(PetscFree(b->values));
    PetscCall(PetscFree(b->comps));
    PetscCall(PetscFree(b));
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSSelectDiscretizations - Copy discretizations to the new `PetscDS` with different field layout

  Not Collective

  Input Parameters:
+ prob      - The `PetscDS` object
. numFields - Number of new fields
. fields    - Old field number for each new field
. minDegree - Minimum degree for a discretization, or `PETSC_DETERMINE` for no limit
- maxDegree - Maximum degree for a discretization, or `PETSC_DETERMINE` for no limit

  Output Parameter:
. newprob - The `PetscDS` copy

  Level: intermediate

.seealso: `PetscDS`, `PetscDSSelectEquations()`, `PetscDSCopyBoundary()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`, `PetscDSSetRiemannSolver()`, `PetscDSSetBdResidual()`, `PetscDSSetBdJacobian()`, `PetscDSCreate()`
@*/
PetscErrorCode PetscDSSelectDiscretizations(PetscDS prob, PetscInt numFields, const PetscInt fields[], PetscInt minDegree, PetscInt maxDegree, PetscDS newprob)
{
  PetscInt Nf, Nfn, fn;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
  if (fields) PetscAssertPointer(fields, 3);
  PetscValidHeaderSpecific(newprob, PETSCDS_CLASSID, 6);
  PetscCall(PetscDSGetNumFields(prob, &Nf));
  PetscCall(PetscDSGetNumFields(newprob, &Nfn));
  numFields = numFields < 0 ? Nf : numFields;
  for (fn = 0; fn < numFields; ++fn) {
    const PetscInt f = fields ? fields[fn] : fn;
    PetscObject    disc;
    PetscClassId   id;

    if (f >= Nf) continue;
    PetscCall(PetscDSGetDiscretization(prob, f, &disc));
    PetscCallContinue(PetscObjectGetClassId(disc, &id));
    if (id == PETSCFE_CLASSID) {
      PetscFE fe;

      PetscCall(PetscFELimitDegree((PetscFE)disc, minDegree, maxDegree, &fe));
      PetscCall(PetscDSSetDiscretization(newprob, fn, (PetscObject)fe));
      PetscCall(PetscFEDestroy(&fe));
    } else {
      PetscCall(PetscDSSetDiscretization(newprob, fn, disc));
    }
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSSelectEquations - Copy pointwise function pointers to the new `PetscDS` with different field layout

  Not Collective

  Input Parameters:
+ prob      - The `PetscDS` object
. numFields - Number of new fields
- fields    - Old field number for each new field

  Output Parameter:
. newprob - The `PetscDS` copy

  Level: intermediate

.seealso: `PetscDS`, `PetscDSSelectDiscretizations()`, `PetscDSCopyBoundary()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`, `PetscDSSetRiemannSolver()`, `PetscDSSetBdResidual()`, `PetscDSSetBdJacobian()`, `PetscDSCreate()`
@*/
PetscErrorCode PetscDSSelectEquations(PetscDS prob, PetscInt numFields, const PetscInt fields[], PetscDS newprob)
{
  PetscInt Nf, Nfn, fn, gn;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
  if (fields) PetscAssertPointer(fields, 3);
  PetscValidHeaderSpecific(newprob, PETSCDS_CLASSID, 4);
  PetscCall(PetscDSGetNumFields(prob, &Nf));
  PetscCall(PetscDSGetNumFields(newprob, &Nfn));
  PetscCheck(numFields <= Nfn, PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_SIZ, "Number of fields %" PetscInt_FMT " to transfer must not be greater than the total number of fields %" PetscInt_FMT, numFields, Nfn);
  for (fn = 0; fn < numFields; ++fn) {
    const PetscInt  f = fields ? fields[fn] : fn;
    PetscPointFn   *obj;
    PetscPointFn   *f0, *f1;
    PetscBdPointFn *f0Bd, *f1Bd;
    PetscRiemannFn *r;

    if (f >= Nf) continue;
    PetscCall(PetscDSGetObjective(prob, f, &obj));
    PetscCall(PetscDSGetResidual(prob, f, &f0, &f1));
    PetscCall(PetscDSGetBdResidual(prob, f, &f0Bd, &f1Bd));
    PetscCall(PetscDSGetRiemannSolver(prob, f, &r));
    PetscCall(PetscDSSetObjective(newprob, fn, obj));
    PetscCall(PetscDSSetResidual(newprob, fn, f0, f1));
    PetscCall(PetscDSSetBdResidual(newprob, fn, f0Bd, f1Bd));
    PetscCall(PetscDSSetRiemannSolver(newprob, fn, r));
    for (gn = 0; gn < numFields; ++gn) {
      const PetscInt     g = fields ? fields[gn] : gn;
      PetscPointJacFn   *g0, *g1, *g2, *g3;
      PetscPointJacFn   *g0p, *g1p, *g2p, *g3p;
      PetscBdPointJacFn *g0Bd, *g1Bd, *g2Bd, *g3Bd;

      if (g >= Nf) continue;
      PetscCall(PetscDSGetJacobian(prob, f, g, &g0, &g1, &g2, &g3));
      PetscCall(PetscDSGetJacobianPreconditioner(prob, f, g, &g0p, &g1p, &g2p, &g3p));
      PetscCall(PetscDSGetBdJacobian(prob, f, g, &g0Bd, &g1Bd, &g2Bd, &g3Bd));
      PetscCall(PetscDSSetJacobian(newprob, fn, gn, g0, g1, g2, g3));
      PetscCall(PetscDSSetJacobianPreconditioner(newprob, fn, gn, g0p, g1p, g2p, g3p));
      PetscCall(PetscDSSetBdJacobian(newprob, fn, gn, g0Bd, g1Bd, g2Bd, g3Bd));
    }
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSCopyEquations - Copy all pointwise function pointers to another `PetscDS`

  Not Collective

  Input Parameter:
. prob - The `PetscDS` object

  Output Parameter:
. newprob - The `PetscDS` copy

  Level: intermediate

.seealso: `PetscDS`, `PetscDSCopyBoundary()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`, `PetscDSSetRiemannSolver()`, `PetscDSSetBdResidual()`, `PetscDSSetBdJacobian()`, `PetscDSCreate()`
@*/
PetscErrorCode PetscDSCopyEquations(PetscDS prob, PetscDS newprob)
{
  PetscWeakForm wf, newwf;
  PetscInt      Nf, Ng;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
  PetscValidHeaderSpecific(newprob, PETSCDS_CLASSID, 2);
  PetscCall(PetscDSGetNumFields(prob, &Nf));
  PetscCall(PetscDSGetNumFields(newprob, &Ng));
  PetscCheck(Nf == Ng, PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_SIZ, "Number of fields must match %" PetscInt_FMT " != %" PetscInt_FMT, Nf, Ng);
  PetscCall(PetscDSGetWeakForm(prob, &wf));
  PetscCall(PetscDSGetWeakForm(newprob, &newwf));
  PetscCall(PetscWeakFormCopy(wf, newwf));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSCopyConstants - Copy all constants set with `PetscDSSetConstants()` to another `PetscDS`

  Not Collective

  Input Parameter:
. prob - The `PetscDS` object

  Output Parameter:
. newprob - The `PetscDS` copy

  Level: intermediate

.seealso: `PetscDS`, `PetscDSCopyBoundary()`, `PetscDSCopyEquations()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`, `PetscDSSetRiemannSolver()`, `PetscDSSetBdResidual()`, `PetscDSSetBdJacobian()`, `PetscDSCreate()`
@*/
PetscErrorCode PetscDSCopyConstants(PetscDS prob, PetscDS newprob)
{
  PetscInt           Nc;
  const PetscScalar *constants;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
  PetscValidHeaderSpecific(newprob, PETSCDS_CLASSID, 2);
  PetscCall(PetscDSGetConstants(prob, &Nc, &constants));
  PetscCall(PetscDSSetConstants(newprob, Nc, (PetscScalar *)constants));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSCopyExactSolutions - Copy all exact solutions set with `PetscDSSetExactSolution()` and `PetscDSSetExactSolutionTimeDerivative()` to another `PetscDS`

  Not Collective

  Input Parameter:
. ds - The `PetscDS` object

  Output Parameter:
. newds - The `PetscDS` copy

  Level: intermediate

.seealso: `PetscDS`, `PetscDSCopyBoundary()`, `PetscDSCopyEquations()`, `PetscDSCopyBounds()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`, `PetscDSSetRiemannSolver()`, `PetscDSSetBdResidual()`, `PetscDSSetBdJacobian()`, `PetscDSCreate()`
@*/
PetscErrorCode PetscDSCopyExactSolutions(PetscDS ds, PetscDS newds)
{
  PetscSimplePointFn *sol;
  void               *ctx;
  PetscInt            Nf, f;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  PetscValidHeaderSpecific(newds, PETSCDS_CLASSID, 2);
  PetscCall(PetscDSGetNumFields(ds, &Nf));
  for (f = 0; f < Nf; ++f) {
    PetscCall(PetscDSGetExactSolution(ds, f, &sol, &ctx));
    PetscCall(PetscDSSetExactSolution(newds, f, sol, ctx));
    PetscCall(PetscDSGetExactSolutionTimeDerivative(ds, f, &sol, &ctx));
    PetscCall(PetscDSSetExactSolutionTimeDerivative(newds, f, sol, ctx));
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PetscDSCopyBounds - Copy lower and upper solution bounds set with `PetscDSSetLowerBound()` and `PetscDSSetLowerBound()` to another `PetscDS`

  Not Collective

  Input Parameter:
. ds - The `PetscDS` object

  Output Parameter:
. newds - The `PetscDS` copy

  Level: intermediate

.seealso: `PetscDS`, `PetscDSCopyBoundary()`, `PetscDSCopyEquations()`, `PetscDSCopyExactSolutions()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`, `PetscDSSetRiemannSolver()`, `PetscDSSetBdResidual()`, `PetscDSSetBdJacobian()`, `PetscDSCreate()`
@*/
PetscErrorCode PetscDSCopyBounds(PetscDS ds, PetscDS newds)
{
  PetscSimplePointFn *bound;
  void               *ctx;
  PetscInt            Nf, f;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  PetscValidHeaderSpecific(newds, PETSCDS_CLASSID, 2);
  PetscCall(PetscDSGetNumFields(ds, &Nf));
  for (f = 0; f < Nf; ++f) {
    PetscCall(PetscDSGetLowerBound(ds, f, &bound, &ctx));
    PetscCall(PetscDSSetLowerBound(newds, f, bound, ctx));
    PetscCall(PetscDSGetUpperBound(ds, f, &bound, &ctx));
    PetscCall(PetscDSSetUpperBound(newds, f, bound, ctx));
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

PetscErrorCode PetscDSCopy(PetscDS ds, PetscInt minDegree, PetscInt maxDegree, DM dmNew, PetscDS dsNew)
{
  DSBoundary b;
  PetscInt   cdim, Nf, f, d;
  PetscBool  isCohesive;
  void      *ctx;

  PetscFunctionBegin;
  PetscCall(PetscDSCopyConstants(ds, dsNew));
  PetscCall(PetscDSCopyExactSolutions(ds, dsNew));
  PetscCall(PetscDSCopyBounds(ds, dsNew));
  PetscCall(PetscDSSelectDiscretizations(ds, PETSC_DETERMINE, NULL, minDegree, maxDegree, dsNew));
  PetscCall(PetscDSCopyEquations(ds, dsNew));
  PetscCall(PetscDSGetNumFields(ds, &Nf));
  for (f = 0; f < Nf; ++f) {
    PetscCall(PetscDSGetContext(ds, f, &ctx));
    PetscCall(PetscDSSetContext(dsNew, f, ctx));
    PetscCall(PetscDSGetCohesive(ds, f, &isCohesive));
    PetscCall(PetscDSSetCohesive(dsNew, f, isCohesive));
    PetscCall(PetscDSGetJetDegree(ds, f, &d));
    PetscCall(PetscDSSetJetDegree(dsNew, f, d));
  }
  if (Nf) {
    PetscCall(PetscDSGetCoordinateDimension(ds, &cdim));
    PetscCall(PetscDSSetCoordinateDimension(dsNew, cdim));
  }
  PetscCall(PetscDSCopyBoundary(ds, PETSC_DETERMINE, NULL, dsNew));
  for (b = dsNew->boundary; b; b = b->next) {
    PetscCall(DMGetLabel(dmNew, b->lname, &b->label));
    /* Do not check if label exists here, since p4est calls this for the reference tree which does not have the labels */
    //PetscCheck(b->label,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Label %s missing in new DM", name);
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

PetscErrorCode PetscDSGetHeightSubspace(PetscDS prob, PetscInt height, PetscDS *subprob)
{
  PetscInt dim, Nf, f;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(prob, PETSCDS_CLASSID, 1);
  PetscAssertPointer(subprob, 3);
  if (height == 0) {
    *subprob = prob;
    PetscFunctionReturn(PETSC_SUCCESS);
  }
  PetscCall(PetscDSGetNumFields(prob, &Nf));
  PetscCall(PetscDSGetSpatialDimension(prob, &dim));
  PetscCheck(height <= dim, PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_OUTOFRANGE, "DS can only handle height in [0, %" PetscInt_FMT "], not %" PetscInt_FMT, dim, height);
  if (!prob->subprobs) PetscCall(PetscCalloc1(dim, &prob->subprobs));
  if (!prob->subprobs[height - 1]) {
    PetscInt cdim;

    PetscCall(PetscDSCreate(PetscObjectComm((PetscObject)prob), &prob->subprobs[height - 1]));
    PetscCall(PetscDSGetCoordinateDimension(prob, &cdim));
    PetscCall(PetscDSSetCoordinateDimension(prob->subprobs[height - 1], cdim));
    for (f = 0; f < Nf; ++f) {
      PetscFE      subfe;
      PetscObject  obj;
      PetscClassId id;

      PetscCall(PetscDSGetDiscretization(prob, f, &obj));
      PetscCall(PetscObjectGetClassId(obj, &id));
      PetscCheck(id == PETSCFE_CLASSID, PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_WRONG, "Unsupported discretization type for field %" PetscInt_FMT, f);
      PetscCall(PetscFEGetHeightSubspace((PetscFE)obj, height, &subfe));
      PetscCall(PetscDSSetDiscretization(prob->subprobs[height - 1], f, (PetscObject)subfe));
    }
  }
  *subprob = prob->subprobs[height - 1];
  PetscFunctionReturn(PETSC_SUCCESS);
}

PetscErrorCode PetscDSPermuteQuadPoint(PetscDS ds, PetscInt ornt, PetscInt field, PetscInt q, PetscInt *qperm)
{
  IS              permIS;
  PetscQuadrature quad;
  DMPolytopeType  ct;
  const PetscInt *perm;
  PetscInt        Na, Nq;

  PetscFunctionBeginHot;
  PetscCall(PetscFEGetQuadrature((PetscFE)ds->disc[field], &quad));
  PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &Nq, NULL, NULL));
  PetscCall(PetscQuadratureGetCellType(quad, &ct));
  PetscCheck(q >= 0 && q < Nq, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Quadrature point %" PetscInt_FMT " is not in [0, %" PetscInt_FMT ")", q, Nq);
  Na = DMPolytopeTypeGetNumArrangements(ct) / 2;
  PetscCheck(ornt >= -Na && ornt < Na, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Orientation %" PetscInt_FMT " of %s is not in [%" PetscInt_FMT ", %" PetscInt_FMT ")", ornt, DMPolytopeTypes[ct], -Na, Na);
  if (!ds->quadPerm[(PetscInt)ct]) PetscCall(PetscQuadratureComputePermutations(quad, NULL, &ds->quadPerm[(PetscInt)ct]));
  permIS = ds->quadPerm[(PetscInt)ct][ornt + Na];
  PetscCall(ISGetIndices(permIS, &perm));
  *qperm = perm[q];
  PetscCall(ISRestoreIndices(permIS, &perm));
  PetscFunctionReturn(PETSC_SUCCESS);
}

PetscErrorCode PetscDSGetDiscType_Internal(PetscDS ds, PetscInt f, PetscDiscType *disctype)
{
  PetscObject  obj;
  PetscClassId id;
  PetscInt     Nf;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  PetscAssertPointer(disctype, 3);
  *disctype = PETSC_DISC_NONE;
  PetscCall(PetscDSGetNumFields(ds, &Nf));
  PetscCheck(f < Nf, PetscObjectComm((PetscObject)ds), PETSC_ERR_ARG_SIZ, "Field %" PetscInt_FMT " must be in [0, %" PetscInt_FMT ")", f, Nf);
  PetscCall(PetscDSGetDiscretization(ds, f, &obj));
  if (obj) {
    PetscCall(PetscObjectGetClassId(obj, &id));
    if (id == PETSCFE_CLASSID) *disctype = PETSC_DISC_FE;
    else *disctype = PETSC_DISC_FV;
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PetscDSDestroy_Basic(PetscDS ds)
{
  PetscFunctionBegin;
  PetscCall(PetscFree(ds->data));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PetscDSInitialize_Basic(PetscDS ds)
{
  PetscFunctionBegin;
  ds->ops->setfromoptions = NULL;
  ds->ops->setup          = NULL;
  ds->ops->view           = NULL;
  ds->ops->destroy        = PetscDSDestroy_Basic;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*MC
  PETSCDSBASIC = "basic" - A discrete system with pointwise residual and boundary residual functions

  Level: intermediate

.seealso: `PetscDSType`, `PetscDSCreate()`, `PetscDSSetType()`
M*/

PETSC_EXTERN PetscErrorCode PetscDSCreate_Basic(PetscDS ds)
{
  PetscDS_Basic *b;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(ds, PETSCDS_CLASSID, 1);
  PetscCall(PetscNew(&b));
  ds->data = b;

  PetscCall(PetscDSInitialize_Basic(ds));
  PetscFunctionReturn(PETSC_SUCCESS);
}