# Commits

committed 30d5754

finished up adding serial DMPlex test with a convergance test on a 5-point Laplacian finite volume test problem

# src/ts/examples/tutorials/ex32.c

-static char help[] = "Second Order Finite Volume Example of Laplacian div \alpha grad u = f.\n";
-/*F
-
-  Multigrid Poisson solver:
-     Laplace(u) = f
-
-  A numerical test on page 64 of
-    A MULTIGRID TUTORIAL by Briggs, Henson & McCormick'
-  u := (x^2-x^4)*(y^4-y^2)*(z^4-z^2) with alpha(x) = -1
-
-  Course grid mesh is read in from an ExodusII file, usually generated by Cubit.
-F*/
-#include <petscts.h>
-#include <petscdmplex.h>
-#include <petscsf.h>
-#include <petscblaslapack.h>
-#if defined(PETSC_HAVE_EXODUSII)
-#include <exodusII.h>
-#else
-#error This example requires ExodusII support. Reconfigure using --download-exodusii
-#endif
-
-#define DIM 2                   /* Geometric dimension */
-#define ALEN(a) (sizeof(a)/sizeof((a)[0]))
-
-/* Represents continuum physical equations. */
-typedef struct _n_Physics *Physics;
-
-/* Physical model includes boundary conditions, initial conditions, and functionals of interest. It is
- * discretization-independent, but its members depend on the scenario being solved. */
-typedef struct _n_Model *Model;
-
-/* 'User' implements a discretization of a continuous model. */
-typedef struct _n_User *User;
-
-typedef PetscErrorCode (*SolutionFunction)(Model,PetscReal,const PetscReal*,PetscScalar*,void*);
-typedef PetscErrorCode (*BoundaryFunction)(Model,PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*);
-typedef PetscErrorCode (*FunctionalFunction)(Model,PetscReal,const PetscReal*,const PetscScalar*,PetscReal*,void*);
-static PetscErrorCode ModelBoundaryRegister(Model,const char*,BoundaryFunction,void*,PetscInt,const PetscInt*);
-static PetscErrorCode ModelSolutionSetDefault(Model,SolutionFunction,void*);
-static PetscErrorCode ModelSourceSetDefault(Model,SolutionFunction,void*);
-static PetscErrorCode ModelAlphaSetDefault(Model,SolutionFunction,void*);
-static PetscErrorCode ModelFunctionalRegister(Model,const char*,PetscInt*,FunctionalFunction,void*);
-static PetscErrorCode OutputVTK(DM,const char*,PetscViewer*);
-
-struct FieldDescription {
-  const char *name;
-  PetscInt dof;
-};
-
-typedef struct _n_BoundaryLink *BoundaryLink;
-struct _n_BoundaryLink {
-  char             *name;
-  BoundaryFunction func;
-  void             *ctx;
-  PetscInt         numids;
-  PetscInt         *ids;
-  BoundaryLink     next;
-};
-
-typedef struct _n_FunctionalLink *FunctionalLink;
-struct _n_FunctionalLink {
-  char               *name;
-  FunctionalFunction func;
-  void               *ctx;
-  PetscInt           offset;
-  FunctionalLink     next;
-};
-
-struct _n_Physics {
-  PetscInt        dof;          /* number of degrees of freedom per cell */
-  void            *data;
-  PetscInt        nfields;
-  const struct FieldDescription *field_desc;
-};
-
-struct _n_Model {
-  MPI_Comm         comm;        /* Does not do collective communicaton, but some error conditions can be collective */
-  Physics          physics;
-  BoundaryLink     boundary;
-  FunctionalLink   functionalRegistry;
-  PetscInt         maxComputed;
-  PetscInt         numMonitored;
-  FunctionalLink   *functionalMonitored;
-  PetscInt         numCall;
-  FunctionalLink   *functionalCall;
-  SolutionFunction solution;
-  SolutionFunction source;
-  SolutionFunction alpha;
-  void             *solutionctx;
-};
-
-struct _n_User {
-  PetscInt       numGhostCells;
-  PetscInt       cEndInterior; /* First boundary ghost cell */
-  Vec            cellgeom, facegeom;
-  PetscInt       vtkInterval;   /* For monitor */
-  Model          model;
-  PetscReal      errorInf[32], error2[32];
-  PetscInt       N[32],ilev;
-};
-
-typedef struct {
-  PetscScalar normal[DIM];              /* Area-scaled normals */
-  PetscScalar centroid[DIM];            /* Location of centroid (quadrature point) */
-} FaceGeom;
-
-typedef struct {
-  PetscScalar centroid[DIM];
-  PetscScalar volume;
-} CellGeom;
-
-PETSC_STATIC_INLINE PetscScalar DotD(const PetscScalar *x, const PetscScalar *y) {PetscScalar sum = 0.0; PetscInt d; for (d = 0; d < DIM; ++d) sum += x[d]*y[d]; return sum;}
-PETSC_STATIC_INLINE PetscReal Norm2D(const PetscScalar *x) {return PetscSqrtReal(PetscAbsScalar(DotD(x,x)));}
-PETSC_STATIC_INLINE void axD(const PetscScalar a,PetscScalar *x){PetscInt i; for (i=0; i<DIM; i++) x[i] *= a;}
-PETSC_STATIC_INLINE void waxD(const PetscScalar a,const PetscScalar *x, PetscScalar *w){PetscInt i; for (i=0; i<DIM; i++) w[i] = x[i]*a;}
-PETSC_STATIC_INLINE void NormalizeD(PetscScalar *x) {PetscReal a = 1./Norm2D(x); PetscInt d; for (d = 0; d < DIM; ++d) x[d] *= a;}
-PETSC_STATIC_INLINE void WaxpyD(PetscScalar a, const PetscScalar *x, const PetscScalar *y, PetscScalar *w) {PetscInt d; for (d = 0; d < DIM; ++d) w[d] = a*x[d] + y[d];}
-
-/******************* Laplacian ********************/
-typedef struct {
-  struct {
-    PetscInt Error;
-  } functional;
-} Physics_Lap;
-
-static const struct FieldDescription PhysicsFields_Lap[] = {{"U",1},{NULL,0}};
-
-#undef __FUNCT__
-#define __FUNCT__ "PhysicsBoundary_DiriHomo"
-static PetscErrorCode PhysicsBoundary_DiriHomo(Model mod, PetscReal time, const PetscReal *fCentroid, const PetscReal *normal, const PetscScalar *xInterior, PetscScalar *xGhost, void *ctx)
-{
-  /* Physics      phys = (Physics)ctx; */
-  /* Physics_Lap *lap = (Physics_Lap*)phys->data; */
-  PetscFunctionBeginUser;
-  xGhost[0] = -xInterior[0];
-  PetscFunctionReturn(0);
-}
-
-/* U: (x^2-x^4)*(y^4-y^2)*(z^4-z^2) */
-#undef __FUNCT__
-#define __FUNCT__ "PhysicsSolution_Lap_x4_x2"
-static PetscErrorCode PhysicsSolution_Lap_x4_x2(Model mod,PetscReal time,const PetscReal *x,PetscScalar *u,void *ctx)
-{
-  /* Physics        phys    = (Physics)ctx; */
-  /* Physics_Lap *lap = (Physics_Lap*)phys->data; */
-  PetscScalar r;
-  PetscScalar e[DIM];
-  PetscInt d;
-
-  PetscFunctionBeginUser;
-  for (d = 0; d < DIM; ++d) e[d] = x[d]*x[d];
-  for (d = 0; d < DIM; ++d) e[d] *= 1 - e[d];
-  for (r = 1.0, d = 0; d < DIM; ++d) r *= e[d];
-  u[0] = r;
-  PetscFunctionReturn(0);
-}
-
-/* RHS: Lap((x^2-x^4)*(y^4-y^2)*(z^4-z^2)) */
-#undef __FUNCT__
-#define __FUNCT__ "PhysicsSource_Lap_x4_x2"
-static PetscErrorCode PhysicsSource_Lap_x4_x2(Model mod,PetscReal time,const PetscReal *x,PetscScalar *rhs,void *ctx)
-{
-  /* Physics        phys    = (Physics)ctx; */
-  /* Physics_Lap *lap = (Physics_Lap*)phys->data; */
-  PetscScalar r;
-  PetscScalar x2[DIM]; /* = x*x; */
-  PetscScalar a[DIM]; /* = x2*(1-x2) */
-  PetscScalar b[DIM]; /* = 2-12*x2; */
-  PetscInt d;
-
-  PetscFunctionBeginUser;
-  for (d = 0; d < DIM; ++d) x2[d] = x[d]*x[d];
-  for (d = 0; d < DIM; ++d) a[d] = x2[d]*(1-x2[d]);
-  for (d = 0; d < DIM; ++d) b[d] = 2-12*x2[d];
-
-  /* cross products of coordinates prevents the use of D_TERM */
-  if (DIM==2)
-    r = b[0]*a[1] + a[0]*b[1];
-  else if (DIM==3)
-    r = b[0]*a[1]*a[2] + a[0]*b[1]*a[2] + a[0]*a[1]*b[2];
-  else  if (DIM==1)
-   r = b[0];
-  else
-    SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for solution this DIM %d",DIM);
-  rhs[0] = -r; /* keep solution positive */
-  PetscFunctionReturn(0);
-}
-#undef __FUNCT__
-#define __FUNCT__ "PhysicsAlpha_n1"
-static PetscErrorCode PhysicsAlpha_n1(Model mod,PetscReal time,const PetscReal *x,PetscScalar *alpha,void *ctx)
-{
-  /* Physics        phys    = (Physics)ctx; */
-  /* Physics_Lap *lap = (Physics_Lap*)phys->data; */
-  PetscFunctionBeginUser;
-  alpha[0] = -1.0;
-  PetscFunctionReturn(0);
-}
-#undef __FUNCT__
-#define __FUNCT__ "PhysicsFunctional_Error"
-static PetscErrorCode PhysicsFunctional_Error(Model mod,PetscReal time,const PetscScalar *x,const PetscScalar *u,PetscReal *f,void *ctx)
-{
-  Physics        phys    = (Physics)ctx;
-  Physics_Lap   *lap = (Physics_Lap*)phys->data;
-  PetscScalar    yexact[1];
-  PetscErrorCode ierr;
-
-  PetscFunctionBeginUser;
-  ierr = (*mod->solution)(mod,time,x,yexact,phys);CHKERRQ(ierr);
-  /* f[lap->functional.Error] = PetscAbsScalar(u[0]-yexact[0]); */
-  f[lap->functional.Error] = PetscRealPart(u[0]-yexact[0]); /* lets do real error for plotting */
-  PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "PhysicsCreate_Lap"
-static PetscErrorCode PhysicsCreate_Lap(Model mod,Physics phys)
-{
-  Physics_Lap *lap;
-  PetscErrorCode ierr;
-
-  PetscFunctionBeginUser;
-  phys->field_desc = PhysicsFields_Lap;
-  ierr = PetscNew(Physics_Lap,&phys->data);CHKERRQ(ierr);
-  lap = (Physics_Lap*)phys->data;
-
-  {
-    const PetscInt diriids[] = {100,101,102,103};
-    /* Register "canned" boundary conditions and defaults for where to apply. */
-    ierr = ModelBoundaryRegister(mod,"diri",PhysicsBoundary_DiriHomo,phys,ALEN(diriids),diriids);CHKERRQ(ierr);
-    /* Initial/transient solution with default boundary conditions */
-    ierr = ModelSolutionSetDefault(mod,PhysicsSolution_Lap_x4_x2,phys);CHKERRQ(ierr);
-    ierr = ModelSourceSetDefault(mod,PhysicsSource_Lap_x4_x2,phys);CHKERRQ(ierr);
-    ierr = ModelAlphaSetDefault(mod,PhysicsAlpha_n1,phys);CHKERRQ(ierr);
-    /* Register "canned" functionals */
-    ierr = ModelFunctionalRegister(mod,"error",&lap->functional.Error,PhysicsFunctional_Error,phys);CHKERRQ(ierr);
-  }
-  PetscFunctionReturn(0);
-}
-
-/* #undef __FUNCT__ */
-/* #define __FUNCT__ "IsExteriorOrGhostFace" */
-/* static PetscErrorCode IsExteriorOrGhostFace(DM dm,PetscInt face,PetscBool *isghost) */
-/* { */
-/*   PetscErrorCode ierr; */
-/*   PetscInt       ghost,boundary; */
-
-/*   PetscFunctionBeginUser; */
-/*   *isghost = PETSC_FALSE; */
-/*   ierr = DMPlexGetLabelValue(dm, "ghost", face, &ghost);CHKERRQ(ierr); */
-/*   ierr = DMPlexGetLabelValue(dm, "Face Sets", face, &boundary);CHKERRQ(ierr); */
-/*   if (ghost >= 0 || boundary >= 0) *isghost = PETSC_TRUE; */
-/*   PetscFunctionReturn(0); */
-/* } */
-
-#undef __FUNCT__
-#define __FUNCT__ "ConstructGeometry"
-/* Set up face data and cell data */
-PetscErrorCode ConstructGeometry(DM dm, Vec *facegeom, Vec *cellgeom, User user)
-{
-  DM             dmFace, dmCell;
-  DMLabel        ghostLabel;
-  PetscSection   sectionFace, sectionCell;
-  PetscSection   coordSection;
-  Vec            coordinates;
-  PetscScalar    *fgeom, *cgeom;
-  PetscInt       dim, cStart, cEnd, c, fStart, fEnd, f;
-  PetscErrorCode ierr;
-
-  PetscFunctionBeginUser;
-  ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
-  if (dim != DIM) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No support for dim %D != DIM %D",dim,DIM);
-  ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
-  ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
-
-  /* Make cell centroids and volumes */
-  ierr = DMClone(dm, &dmCell);CHKERRQ(ierr);
-  ierr = DMSetCoordinateSection(dmCell, coordSection);CHKERRQ(ierr);
-  ierr = DMSetCoordinatesLocal(dmCell, coordinates);CHKERRQ(ierr);
-  ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &sectionCell);CHKERRQ(ierr);
-  ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
-  ierr = PetscSectionSetChart(sectionCell, cStart, cEnd);CHKERRQ(ierr);
-  for (c = cStart; c < cEnd; ++c) {
-    ierr = PetscSectionSetDof(sectionCell, c, sizeof(CellGeom)/sizeof(PetscScalar));CHKERRQ(ierr);
-  }
-  ierr = PetscSectionSetUp(sectionCell);CHKERRQ(ierr);
-  ierr = DMSetDefaultSection(dmCell, sectionCell);CHKERRQ(ierr);
-  ierr = PetscSectionDestroy(&sectionCell);CHKERRQ(ierr); /* relinquish our reference */
-
-  ierr = DMCreateLocalVector(dmCell, cellgeom);CHKERRQ(ierr);
-  ierr = VecGetArray(*cellgeom, &cgeom);CHKERRQ(ierr);
-  for (c = cStart; c < user->cEndInterior; ++c) {
-    CellGeom *cg;
-
-    ierr = DMPlexPointLocalRef(dmCell, c, cgeom, &cg);CHKERRQ(ierr);
-    ierr = PetscMemzero(cg,sizeof(*cg));CHKERRQ(ierr);
-    ierr = DMPlexComputeCellGeometryFVM(dmCell, c, &cg->volume, cg->centroid, NULL);CHKERRQ(ierr);
-  }
-  /* Compute face normals, centroid(used?) */
-  ierr = DMClone(dm, &dmFace);CHKERRQ(ierr);
-  ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &sectionFace);CHKERRQ(ierr);
-  ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr);
-  ierr = PetscSectionSetChart(sectionFace, fStart, fEnd);CHKERRQ(ierr);
-
-  for (f = fStart; f < fEnd; ++f) {
-    ierr = PetscSectionSetDof(sectionFace, f, sizeof(FaceGeom)/sizeof(PetscScalar));CHKERRQ(ierr);
-  }
-  ierr = PetscSectionSetUp(sectionFace);CHKERRQ(ierr);
-  ierr = DMSetDefaultSection(dmFace, sectionFace);CHKERRQ(ierr);
-  ierr = PetscSectionDestroy(&sectionFace);CHKERRQ(ierr);
-  ierr = DMCreateLocalVector(dmFace, facegeom);CHKERRQ(ierr);
-  ierr = VecGetArray(*facegeom, &fgeom);CHKERRQ(ierr);
-  ierr = DMPlexGetLabel(dm, "ghost", &ghostLabel);CHKERRQ(ierr);
-  for (f = fStart; f < fEnd; ++f) {
-    FaceGeom *fg;
-    PetscInt  ghost;
-
-    ierr = DMLabelGetValue(ghostLabel, f, &ghost);CHKERRQ(ierr);
-    if (ghost >= 0) continue;
-    ierr = DMPlexPointLocalRef(dmFace, f, fgeom, &fg);CHKERRQ(ierr);
-    ierr = DMPlexComputeCellGeometryFVM(dm, f, NULL, fg->centroid, fg->normal);CHKERRQ(ierr);
-    /* Flip face orientation if necessary to match ordering in support, and Update minimum radius */
-    {
-      CellGeom       *cL, *cR;
-      const PetscInt *cells;
-      PetscReal      *lcentroid, *rcentroid;
-      PetscScalar     v[3];
-      PetscInt        d;
-
-      ierr = DMPlexGetSupport(dm, f, &cells);CHKERRQ(ierr);
-      ierr = DMPlexPointLocalRead(dmCell, cells[0], cgeom, &cL);CHKERRQ(ierr);
-      ierr = DMPlexPointLocalRead(dmCell, cells[1], cgeom, &cR);CHKERRQ(ierr);
-      lcentroid = cells[0] >= user->cEndInterior ? fg->centroid : cL->centroid;
-      rcentroid = cells[1] >= user->cEndInterior ? fg->centroid : cR->centroid;
-      WaxpyD(-1, lcentroid, rcentroid, v);
-      if (DotD(fg->normal, v) < 0) {
-        for (d = 0; d < dim; ++d) fg->normal[d] = -fg->normal[d];
-      }
-      if (DotD(fg->normal, v) <= 0) {
-#if DIM == 2
-        SETERRQ5(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Direction for face %d could not be fixed, normal (%g,%g) v (%g,%g)", f, fg->normal[0], fg->normal[1], v[0], v[1]);
-#elif DIM == 3
-        SETERRQ7(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Direction for face %d could not be fixed, normal (%g,%g,%g) v (%g,%g,%g)", f, fg->normal[0], fg->normal[1], fg->normal[2], v[0], v[1], v[2]);
-#else
-#  error DIM not supported
-#endif
-      }
-    }
-  }
-  /* Compute centroids of ghost cells */
-  for (c = user->cEndInterior; c < cEnd; ++c) {
-    FaceGeom       *fg;
-    const PetscInt *cone,    *support;
-    PetscInt        coneSize, supportSize, s;
-
-    ierr = DMPlexGetConeSize(dmCell, c, &coneSize);CHKERRQ(ierr);
-    if (coneSize != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ghost cell %d has cone size %d != 1", c, coneSize);
-    ierr = DMPlexGetCone(dmCell, c, &cone);CHKERRQ(ierr);
-    ierr = DMPlexGetSupportSize(dmCell, cone[0], &supportSize);CHKERRQ(ierr);
-    if (supportSize != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %d has support size %d != 1", cone[0], supportSize);
-    ierr = DMPlexGetSupport(dmCell, cone[0], &support);CHKERRQ(ierr);
-    ierr = DMPlexPointLocalRef(dmFace, cone[0], fgeom, &fg);CHKERRQ(ierr);
-    for (s = 0; s < 2; ++s) {
-      /* Reflect ghost centroid across plane of face */
-      if (support[s] == c) {
-        const CellGeom *ci;
-        CellGeom       *cg;
-        PetscScalar     c2f[3], a;
-
-        ierr = DMPlexPointLocalRead(dmCell, support[(s+1)%2], cgeom, &ci);CHKERRQ(ierr);
-        WaxpyD(-1, ci->centroid, fg->centroid, c2f); /* cell to face centroid */
-        a    = DotD(c2f, fg->normal)/DotD(fg->normal, fg->normal);
-        ierr = DMPlexPointLocalRef(dmCell, support[s], cgeom, &cg);CHKERRQ(ierr);
-        WaxpyD(2*a, fg->normal, ci->centroid, cg->centroid);
-        cg->volume = ci->volume;
-      }
-    }
-  }
-  ierr = VecRestoreArray(*facegeom, &fgeom);CHKERRQ(ierr);
-  ierr = VecRestoreArray(*cellgeom, &cgeom);CHKERRQ(ierr);
-  ierr = DMDestroy(&dmCell);CHKERRQ(ierr);
-  ierr = DMDestroy(&dmFace);CHKERRQ(ierr);
-  PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "SetUpLocalSpace"
-PetscErrorCode SetUpLocalSpace(DM dm, User user)
-{
-  PetscSection   stateSection;
-  Physics        phys;
-  PetscInt       dof = user->model->physics->dof, *cind, d, stateSize, cStart, cEnd, c, i;
-  PetscErrorCode ierr;
-
-  PetscFunctionBeginUser;
-  ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
-  ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &stateSection);CHKERRQ(ierr);
-  phys = user->model->physics;
-  ierr = PetscSectionSetNumFields(stateSection,phys->nfields);CHKERRQ(ierr);
-  for (i=0; i<phys->nfields; i++) {
-    ierr = PetscSectionSetFieldName(stateSection,i,phys->field_desc[i].name);CHKERRQ(ierr);
-    ierr = PetscSectionSetFieldComponents(stateSection,i,phys->field_desc[i].dof);CHKERRQ(ierr);
-  }
-  ierr = PetscSectionSetChart(stateSection, cStart, cEnd);CHKERRQ(ierr);
-
-  for (c = cStart; c < cEnd; ++c) {
-    for (i=0; i<phys->nfields; i++) {
-      ierr = PetscSectionSetFieldDof(stateSection,c,i,phys->field_desc[i].dof);CHKERRQ(ierr);
-    }
-    ierr = PetscSectionSetDof(stateSection, c, dof);CHKERRQ(ierr);
-  }
-  for (c = user->cEndInterior; c < cEnd; ++c) {
-    ierr = PetscSectionSetConstraintDof(stateSection, c, dof);CHKERRQ(ierr);
-  }
-  ierr = PetscSectionSetUp(stateSection);CHKERRQ(ierr);
-  ierr = PetscMalloc(dof * sizeof(PetscInt), &cind);CHKERRQ(ierr);
-  for (d = 0; d < dof; ++d) cind[d] = d;
-  for (c = user->cEndInterior; c < cEnd; ++c) {
-    ierr = PetscSectionSetConstraintIndices(stateSection, c, cind);CHKERRQ(ierr);
-  }
-  ierr = PetscFree(cind);CHKERRQ(ierr);
-  ierr = PetscSectionGetStorageSize(stateSection, &stateSize);CHKERRQ(ierr);
-  ierr = DMSetDefaultSection(dm,stateSection);CHKERRQ(ierr);
-  ierr = PetscSectionDestroy(&stateSection);CHKERRQ(ierr);
-  PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "SetUpBoundaries"
-PetscErrorCode SetUpBoundaries(MPI_Comm comm, User user)
-{
-  Model          mod = user->model;
-  PetscErrorCode ierr;
-  BoundaryLink   b;
-
-  PetscFunctionBeginUser;
-  ierr = PetscOptionsBegin(comm,NULL,"Boundary condition options","");CHKERRQ(ierr);
-  for (b = mod->boundary; b; b=b->next) {
-    char      optname[512];
-    PetscInt  ids[512],len = 512;
-    PetscBool flg;
-    ierr = PetscSNPrintf(optname,sizeof optname,"-bc_%s",b->name);CHKERRQ(ierr);
-    ierr = PetscMemzero(ids,sizeof(ids));CHKERRQ(ierr);
-    ierr = PetscOptionsIntArray(optname,"List of boundary IDs","",ids,&len,&flg);CHKERRQ(ierr);
-    if (flg) {
-      /* TODO: check all IDs to make sure they exist in the mesh */
-      ierr      = PetscFree(b->ids);CHKERRQ(ierr);
-      b->numids = len;
-      ierr      = PetscMalloc(len*sizeof(PetscInt),&b->ids);CHKERRQ(ierr);
-      ierr      = PetscMemcpy(b->ids,ids,len*sizeof(PetscInt));CHKERRQ(ierr);
-    }
-  }
-  ierr = PetscOptionsEnd();CHKERRQ(ierr);
-  PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "ModelBoundaryRegister"
-/* The ids are just defaults, can be overridden at command line. I expect to set defaults based on names in the future. */
-static PetscErrorCode ModelBoundaryRegister(Model mod,const char *name,BoundaryFunction bcFunc,void *ctx,PetscInt numids,const PetscInt *ids)
-{
-  PetscErrorCode ierr;
-  BoundaryLink   link;
-
-  PetscFunctionBeginUser;
-  ierr          = PetscNew(struct _n_BoundaryLink,&link);CHKERRQ(ierr);
-  ierr          = PetscStrallocpy(name,&link->name);CHKERRQ(ierr);
-  link->numids  = numids;
-  ierr          = PetscMalloc(numids*sizeof(PetscInt),&link->ids);CHKERRQ(ierr);
-  ierr          = PetscMemcpy(link->ids,ids,numids*sizeof(PetscInt));CHKERRQ(ierr);
-  link->func    = bcFunc;
-  link->ctx     = ctx;
-  link->next    = mod->boundary;
-  mod->boundary = link;
-  PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "BoundaryLinkDestroy"
-static PetscErrorCode BoundaryLinkDestroy(BoundaryLink *link)
-{
-  PetscErrorCode ierr;
-  BoundaryLink   l,next;
-
-  PetscFunctionBeginUser;
-  if (!link) PetscFunctionReturn(0);
-  l     = *link;
-  *link = NULL;
-  for (; l; l=next) {
-    next = l->next;
-    ierr = PetscFree(l->ids);CHKERRQ(ierr);
-    ierr = PetscFree(l->name);CHKERRQ(ierr);
-    ierr = PetscFree(l);CHKERRQ(ierr);
-  }
-  PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "ModelBoundaryFind"
-static PetscErrorCode ModelBoundaryFind(Model mod,PetscInt id,BoundaryFunction *bcFunc,void **ctx)
-{
-  BoundaryLink link;
-  PetscInt     i;
-
-  PetscFunctionBeginUser;
-  *bcFunc = NULL;
-  for (link=mod->boundary; link; link=link->next) {
-    for (i=0; i<link->numids; i++) {
-      if (link->ids[i] == id) {
-        *bcFunc = link->func;
-        *ctx    = link->ctx;
-        PetscFunctionReturn(0);
-      }
-    }
-  }
-  SETERRQ1(mod->comm,PETSC_ERR_USER,"Boundary ID %D not associated with any registered boundary condition",id);
-  PetscFunctionReturn(0);
-}
-#undef __FUNCT__
-#define __FUNCT__ "ModelSolutionSetDefault"
-/* Behavior will be different for multi-physics or when using non-default boundary conditions */
-static PetscErrorCode ModelSolutionSetDefault(Model mod,SolutionFunction func,void *ctx)
-{
-  PetscFunctionBeginUser;
-  mod->solution    = func;
-  mod->solutionctx = ctx;
-  PetscFunctionReturn(0);
-}
-#undef __FUNCT__
-#define __FUNCT__ "ModelSourceSetDefault"
-/* Behavior will be different for multi-physics or when using non-default boundary conditions */
-static PetscErrorCode ModelSourceSetDefault(Model mod,SolutionFunction func,void *ctx)
-{
-  PetscFunctionBeginUser;
-  mod->source    = func;
-  /* mod->solutionctx = ctx; */
-  PetscFunctionReturn(0);
-}
-#undef __FUNCT__
-#define __FUNCT__ "ModelAlphaSetDefault"
-/* Behavior will be different for multi-physics or when using non-default boundary conditions */
-static PetscErrorCode ModelAlphaSetDefault(Model mod,SolutionFunction func,void *ctx)
-{
-  PetscFunctionBeginUser;
-  mod->alpha    = func;
-  /* mod->solutionctx = ctx; */
-  PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "ModelFunctionalRegister"
-static PetscErrorCode ModelFunctionalRegister(Model mod,const char *name,PetscInt *offset,FunctionalFunction func,void *ctx)
-{
-  PetscErrorCode ierr;
-  FunctionalLink link,*ptr;
-  PetscInt       lastoffset = -1;
-
-  PetscFunctionBeginUser;
-  for (ptr=&mod->functionalRegistry; *ptr; ptr = &(*ptr)->next) lastoffset = (*ptr)->offset;
-  ierr         = PetscNew(struct _n_FunctionalLink,&link);CHKERRQ(ierr);
-  ierr         = PetscStrallocpy(name,&link->name);CHKERRQ(ierr);
-  link->offset = lastoffset + 1;
-  link->func   = func;
-  link->ctx    = ctx;
-  link->next   = NULL;
-  *ptr         = link;
-  *offset      = link->offset;
-  PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "ModelFunctionalSetFromOptions"
-static PetscErrorCode ModelFunctionalSetFromOptions(Model mod)
-{
-  PetscErrorCode ierr;
-  PetscInt       i,j;
-  FunctionalLink link;
-  char           *names[256];
-
-  PetscFunctionBeginUser;
-  mod->numMonitored = ALEN(names);
-  ierr = PetscOptionsStringArray("-monitor","list of functionals to monitor","",names,&mod->numMonitored,NULL);CHKERRQ(ierr);
-  /* Create list of functionals that will be computed somehow */
-  ierr = PetscMalloc(mod->numMonitored*sizeof(FunctionalLink),&mod->functionalMonitored);CHKERRQ(ierr);
-  /* Create index of calls that we will have to make to compute these functionals (over-allocation in general). */
-  ierr = PetscMalloc(mod->numMonitored*sizeof(FunctionalLink),&mod->functionalCall);CHKERRQ(ierr);
-  mod->numCall = 0;
-
-  for (i=0; i<mod->numMonitored; i++) {
-    for (link=mod->functionalRegistry; link; link=link->next) {
-      PetscBool match;
-      ierr = PetscStrcasecmp(names[i],link->name,&match);CHKERRQ(ierr);
-      if (match) break;
-    }
-    if (!link) SETERRQ1(mod->comm,PETSC_ERR_USER,"No known functional '%s'",names[i]);
-    mod->functionalMonitored[i] = link;
-    for (j=0; j<i; j++) {
-      if (mod->functionalCall[j]->func == link->func && mod->functionalCall[j]->ctx == link->ctx) goto next_name;
-    }
-    mod->functionalCall[mod->numCall++] = link; /* Just points to the first link using the result. There may be more results. */
-next_name:
-    ierr = PetscFree(names[i]);CHKERRQ(ierr);
-  }
-
-  /* Find out the maximum index of any functional computed by a function we will be calling (even if we are not using it) */
-  mod->maxComputed = -1;
-  for (link=mod->functionalRegistry; link; link=link->next) {
-    for (i=0; i<mod->numCall; i++) {
-      FunctionalLink call = mod->functionalCall[i];
-      if (link->func == call->func && link->ctx == call->ctx) {
-        mod->maxComputed = PetscMax(mod->maxComputed,link->offset);
-      }
-    }
-  }
-  PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "FunctionalLinkDestroy"
-static PetscErrorCode FunctionalLinkDestroy(FunctionalLink *link)
-{
-  PetscErrorCode ierr;
-  FunctionalLink l,next;
-
-  PetscFunctionBeginUser;
-  if (!link) PetscFunctionReturn(0);
-  l     = *link;
-  *link = NULL;
-  for (; l; l=next) {
-    next = l->next;
-    ierr = PetscFree(l->name);CHKERRQ(ierr);
-    ierr = PetscFree(l);CHKERRQ(ierr);
-  }
-  PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "SetInitialCondition"
-PetscErrorCode SetInitialCondition(DM dm, Vec X, User user)
-{
-  PetscScalar       *x;
-  PetscInt          cStart, cEnd, cEndInterior = user->cEndInterior, c;
-  PetscErrorCode    ierr;
-
-  PetscFunctionBeginUser;
-  ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
-  ierr = VecGetArray(X, &x);CHKERRQ(ierr);
-  for (c = cStart; c < cEndInterior; ++c) {
-    PetscScalar    *xc;
-    ierr = DMPlexPointGlobalRef(dm,c,x,&xc);CHKERRQ(ierr);
-    if (xc) *xc = 0.; /* process ghost */
-  }
-  ierr = VecRestoreArray(X, &x);CHKERRQ(ierr);
-  PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "OutputVTK"
-static PetscErrorCode OutputVTK(DM dm, const char *filename, PetscViewer *viewer)
-{
-  PetscErrorCode ierr;
-
-  PetscFunctionBeginUser;
-  ierr = PetscViewerCreate(PetscObjectComm((PetscObject)dm), viewer);CHKERRQ(ierr);
-  ierr = PetscViewerSetType(*viewer, PETSCVIEWERVTK);CHKERRQ(ierr);
-  ierr = PetscViewerFileSetName(*viewer, filename);CHKERRQ(ierr);
-  PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "MonitorVTK"
-static PetscErrorCode MonitorVTK(TS ts,PetscInt stepnum,PetscReal time,Vec X,void *ctx)
-{
-  User           user = (User)ctx;
-  DM             dm;
-  PetscViewer    viewer;
-  char           filename[PETSC_MAX_PATH_LEN],*ftable = NULL;
-  PetscReal      xnorm;
-  Vec            E=0;
-  PetscErrorCode ierr;
-
-  PetscFunctionBeginUser;
-  ierr = PetscObjectSetName((PetscObject) X, "solution");CHKERRQ(ierr);
-  ierr = VecGetDM(X,&dm);CHKERRQ(ierr);
-  ierr = VecNorm(X,NORM_INFINITY,&xnorm);CHKERRQ(ierr);
-  ierr = VecDuplicate(X,&E);CHKERRQ(ierr);
-  if (stepnum >= 0) {           /* No summary for final time */
-    Model             mod = user->model;
-    PetscInt          c,cStart,cEnd,fcount,i;
-    size_t            ftableused,ftablealloc;
-    const PetscScalar *cellgeom_v,*x_v;
-    PetscScalar       *func_v;
-    DM                dmCell;
-    PetscReal         *fmin,*fmax,*fintegral,*ftmp;
-    fcount = mod->maxComputed+1;
-    ierr   = PetscMalloc4(fcount,PetscReal,&fmin,fcount,PetscReal,&fmax,fcount,PetscReal,&fintegral,fcount,PetscReal,&ftmp);CHKERRQ(ierr);
-    for (i=0; i<fcount; i++) {
-      fmin[i]      = PETSC_MAX_REAL;
-      fmax[i]      = PETSC_MIN_REAL;
-      fintegral[i] = 0;
-    }
-    ierr = DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);CHKERRQ(ierr);
-    ierr = VecGetDM(user->cellgeom,&dmCell);CHKERRQ(ierr);
-    ierr = VecGetArrayRead(user->cellgeom,&cellgeom_v);CHKERRQ(ierr);
-    ierr = VecGetArrayRead(X,&x_v);CHKERRQ(ierr);
-    ierr = VecGetArray(E,&func_v);CHKERRQ(ierr);
-    for (c=cStart; c<user->cEndInterior; c++) {
-      const CellGeom    *cg;
-      const PetscScalar *cx;
-      PetscScalar *cfunc;
-      ierr = DMPlexPointLocalRead(dmCell,c,cellgeom_v,&cg);CHKERRQ(ierr);
-      ierr = DMPlexPointGlobalRead(dm,c,x_v,&cx);CHKERRQ(ierr);
-      ierr = DMPlexPointGlobalRef(dm,c,func_v,&cfunc);CHKERRQ(ierr);
-      if (!cx) continue;        /* not a global cell */
-      for (i=0; i<mod->numCall; i++) {
-        FunctionalLink flink = mod->functionalCall[i];
-        ierr = (*flink->func)(mod,time,cg->centroid,cx,&ftmp[i],flink->ctx);CHKERRQ(ierr);
-        *cfunc = ftmp[i]; /* save last functional */
-        if (i==mod->numCall-1) ftmp[i] = PetscAbs(ftmp[i]); /* allready a real, this assumes we want an abs - error */
-      }
-      for (i=0; i<fcount; i++) {
-        fmin[i]       = PetscMin(fmin[i],ftmp[i]);
-        fmax[i]       = PetscMax(fmax[i],ftmp[i]);
-        fintegral[i] += cg->volume * ftmp[i];
-      }
-    }
-    ierr = VecRestoreArrayRead(user->cellgeom,&cellgeom_v);CHKERRQ(ierr);
-    ierr = VecRestoreArrayRead(X,&x_v);CHKERRQ(ierr);
-    ierr = VecRestoreArray(E,&func_v);CHKERRQ(ierr);
-    ierr = MPI_Allreduce(MPI_IN_PLACE,fmin,fcount,MPIU_REAL,MPI_MIN,PetscObjectComm((PetscObject)ts));CHKERRQ(ierr);
-    ierr = MPI_Allreduce(MPI_IN_PLACE,fmax,fcount,MPIU_REAL,MPI_MAX,PetscObjectComm((PetscObject)ts));CHKERRQ(ierr);
-    ierr = MPI_Allreduce(MPI_IN_PLACE,fintegral,fcount,MPIU_REAL,MPI_SUM,PetscObjectComm((PetscObject)ts));CHKERRQ(ierr);
-    if (fcount) { user->errorInf[user->ilev] = fmax[mod->numCall-1]; user->error2[user->ilev] = fintegral[mod->numCall-1]; } /* keep for convergence test */
-    ftablealloc = fcount * 100;
-    ftableused  = 0;
-    ierr        = PetscMalloc(ftablealloc,&ftable);CHKERRQ(ierr);
-    for (i=0; i<mod->numMonitored; i++) {
-      size_t         countused;
-      char           buffer[256],*p;
-      FunctionalLink flink = mod->functionalMonitored[i];
-      PetscInt       id    = flink->offset;
-      if (i % 3) {
-        ierr = PetscMemcpy(buffer,"  ",2);CHKERRQ(ierr);
-        p    = buffer + 2;
-      } else if (i) {
-        char newline[] = "\n";
-        ierr = PetscMemcpy(buffer,newline,sizeof newline-1);CHKERRQ(ierr);
-        p    = buffer + sizeof newline - 1;
-      } else {
-        p = buffer;
-      }
-      ierr = PetscSNPrintfCount(p,sizeof buffer-(p-buffer),"%12s functional: [%10.7G,%10.7G] integral=%10.7G",&countused,flink->name,fmin[id],fmax[id],fintegral[id]);CHKERRQ(ierr);
-      countused += p - buffer;
-      if (countused > ftablealloc-ftableused-1) { /* reallocate */
-        char *ftablenew;
-        ftablealloc = 2*ftablealloc + countused;
-        ierr = PetscMalloc(ftablealloc,&ftablenew);CHKERRQ(ierr);
-        ierr = PetscMemcpy(ftablenew,ftable,ftableused);CHKERRQ(ierr);
-        ierr = PetscFree(ftable);CHKERRQ(ierr);
-        ftable = ftablenew;
-      }
-      ierr = PetscMemcpy(ftable+ftableused,buffer,countused);CHKERRQ(ierr);
-      ftableused += countused;
-      ftable[ftableused] = 0;
-    }
-    ierr = PetscFree4(fmin,fmax,fintegral,ftmp);CHKERRQ(ierr);
-
-    ierr = PetscPrintf(PetscObjectComm((PetscObject)ts),"% 3D  time %8.4G  |x| %8.4G  %s\n",stepnum,time,xnorm,ftable ? ftable : "");CHKERRQ(ierr);
-    ierr = PetscFree(ftable);CHKERRQ(ierr);
-  }
-  if (user->vtkInterval < 1) {
-    ierr = VecDestroy(&E);CHKERRQ(ierr);
-    PetscFunctionReturn(0);
-  }
-  if (user->model->functionalCall && ((stepnum == -1) ^ (stepnum % user->vtkInterval == 0))) {
-    Model          mod   = user->model;
-    FunctionalLink flink = mod->functionalCall[mod->numCall-1];
-    if (stepnum == -1) {        /* Final time is not multiple of normal time interval, write it anyway */
-      ierr = TSGetTimeStepNumber(ts,&stepnum);CHKERRQ(ierr);
-    }
-    ierr = PetscSNPrintf(filename,sizeof filename,"ex32-%03D-u.vtu",stepnum);CHKERRQ(ierr);
-    ierr = OutputVTK(dm,filename,&viewer);CHKERRQ(ierr);
-    ierr = VecView(X,viewer);CHKERRQ(ierr);
-    ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
-    ierr = PetscObjectSetName((PetscObject) E, flink->name);CHKERRQ(ierr);
-    ierr = PetscSNPrintf(filename,sizeof filename,"ex32-%03D-%s.vtu",stepnum,flink->name);CHKERRQ(ierr);
-    ierr = OutputVTK(dm,filename,&viewer);CHKERRQ(ierr);
-    ierr = VecView(E,viewer);CHKERRQ(ierr);
-    ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
-  }
-  ierr = VecDestroy(&E);CHKERRQ(ierr);
-  PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "CreatePartitionVec"
-PetscErrorCode CreatePartitionVec(DM dm, DM *dmCell, Vec *partition)
-{
-  PetscSF        sfPoint;
-  PetscSection   coordSection;
-  Vec            coordinates;
-  PetscSection   sectionCell;
-  PetscScalar    *part;
-  PetscInt       cStart, cEnd, c;
-  PetscMPIInt    rank;
-  PetscErrorCode ierr;
-
-  PetscFunctionBeginUser;
-  ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
-  ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
-  ierr = DMClone(dm, dmCell);CHKERRQ(ierr);
-  ierr = DMGetPointSF(dm, &sfPoint);CHKERRQ(ierr);
-  ierr = DMSetPointSF(*dmCell, sfPoint);CHKERRQ(ierr);
-  ierr = DMSetCoordinateSection(*dmCell, coordSection);CHKERRQ(ierr);
-  ierr = DMSetCoordinatesLocal(*dmCell, coordinates);CHKERRQ(ierr);
-  ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr);
-  ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &sectionCell);CHKERRQ(ierr);
-  ierr = DMPlexGetHeightStratum(*dmCell, 0, &cStart, &cEnd);CHKERRQ(ierr);
-  ierr = PetscSectionSetChart(sectionCell, cStart, cEnd);CHKERRQ(ierr);
-  for (c = cStart; c < cEnd; ++c) {
-    ierr = PetscSectionSetDof(sectionCell, c, 1);CHKERRQ(ierr);
-  }
-  ierr = PetscSectionSetUp(sectionCell);CHKERRQ(ierr);
-  ierr = DMSetDefaultSection(*dmCell, sectionCell);CHKERRQ(ierr);
-  ierr = PetscSectionDestroy(&sectionCell);CHKERRQ(ierr);
-  ierr = DMCreateLocalVector(*dmCell, partition);CHKERRQ(ierr);
-  ierr = PetscObjectSetName((PetscObject)*partition, "partition");CHKERRQ(ierr);
-  ierr = VecGetArray(*partition, &part);CHKERRQ(ierr);
-  for (c = cStart; c < cEnd; ++c) {
-    PetscScalar *p;
-
-    ierr = DMPlexPointLocalRef(*dmCell, c, part, &p);CHKERRQ(ierr);
-    p[0] = rank;
-  }
-  ierr = VecRestoreArray(*partition, &part);CHKERRQ(ierr);
-  PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "ApplyBC"
-static PetscErrorCode ApplyBC(DM dm, PetscReal time, Vec locX, User user)
-{
-  Model             mod = user->model;
-  const char        *name = "Face Sets";
-  DM                dmFace;
-  IS                idIS;
-  const PetscInt    *ids;
-  PetscScalar       *x;
-  const PetscScalar *facegeom;
-  PetscInt          numFS, fs;
-  PetscErrorCode    ierr;
-
-  PetscFunctionBeginUser;
-  ierr = VecGetDM(user->facegeom,&dmFace);CHKERRQ(ierr);
-  ierr = DMPlexGetLabelIdIS(dm, name, &idIS);CHKERRQ(ierr);
-  if (!idIS) PetscFunctionReturn(0);
-  ierr = ISGetLocalSize(idIS, &numFS);CHKERRQ(ierr);
-  ierr = ISGetIndices(idIS, &ids);CHKERRQ(ierr);
-  ierr = VecGetArrayRead(user->facegeom, &facegeom);CHKERRQ(ierr);
-  ierr = VecGetArray(locX, &x);CHKERRQ(ierr);
-  for (fs = 0; fs < numFS; ++fs) {
-    BoundaryFunction bcFunc;
-    void             *bcCtx;
-    IS               faceIS;
-    const PetscInt   *faces;
-    PetscInt         numFaces, f;
-    ierr = ModelBoundaryFind(mod,ids[fs],&bcFunc,&bcCtx);CHKERRQ(ierr);
-    ierr = DMPlexGetStratumIS(dm, name, ids[fs], &faceIS);CHKERRQ(ierr);
-    ierr = ISGetLocalSize(faceIS, &numFaces);CHKERRQ(ierr);
-    ierr = ISGetIndices(faceIS, &faces);CHKERRQ(ierr);
-    for (f = 0; f < numFaces; ++f) {
-      const PetscInt    face = faces[f], *cells;
-      const PetscScalar *xI;
-      PetscScalar       *xG;
-      const FaceGeom    *fg;
-      ierr = DMPlexPointLocalRead(dmFace, face, facegeom, &fg);CHKERRQ(ierr);
-      ierr = DMPlexGetSupport(dm, face, &cells);CHKERRQ(ierr);
-      ierr = DMPlexPointLocalRead(dm, cells[0], x, &xI);CHKERRQ(ierr);
-      ierr = DMPlexPointLocalRef(dm, cells[1], x, &xG);CHKERRQ(ierr);
-      ierr = (*bcFunc)(mod, time, fg->centroid, fg->normal, xI, xG, bcCtx);CHKERRQ(ierr);
-    }
-    ierr = ISRestoreIndices(faceIS, &faces);CHKERRQ(ierr);
-    ierr = ISDestroy(&faceIS);CHKERRQ(ierr);
-  }
-  ierr = VecRestoreArray(locX, &x);CHKERRQ(ierr);
-  ierr = VecRestoreArrayRead(user->facegeom,&facegeom);CHKERRQ(ierr);
-  ierr = ISRestoreIndices(idIS, &ids);CHKERRQ(ierr);
-  ierr = ISDestroy(&idIS);CHKERRQ(ierr);
-  PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "ApplyLaplacianLocal"
-static PetscErrorCode ApplyLaplacianLocal(DM dm,DM dmFace,DM dmCell,PetscReal time,Vec locX,Vec F,User user)
-{
-  /* Physics           phys = user->model->physics; */
-  DMLabel           ghostLabel;
-  PetscErrorCode    ierr;
-  const PetscScalar *facegeom, *cellgeom, *x;
-  PetscScalar       *f;
-  PetscInt          fStart, fEnd, face;
-
-  PetscFunctionBeginUser;
-  ierr = VecGetArrayRead(user->facegeom,&facegeom);CHKERRQ(ierr);
-  ierr = VecGetArrayRead(user->cellgeom,&cellgeom);CHKERRQ(ierr);
-  ierr = VecGetArrayRead(locX,&x);CHKERRQ(ierr);
-  ierr = VecGetArray(F,&f);CHKERRQ(ierr);
-  ierr = DMPlexGetLabel(dm, "ghost", &ghostLabel);CHKERRQ(ierr);
-  ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr);
-  for (face = fStart; face < fEnd; ++face) {
-    const PetscInt    *cells;
-    PetscInt          ghost;
-    PetscScalar       *fL,*fR,v[DIM],dx,area,flux,alpha;
-    const FaceGeom    *fg;
-    const CellGeom    *cgL,*cgR;
-    const PetscScalar *xL,*xR;
-
-    ierr = DMLabelGetValue(ghostLabel, face, &ghost);CHKERRQ(ierr);
-    if (ghost >= 0) continue;
-    ierr = DMPlexGetSupport(dm, face, &cells);CHKERRQ(ierr);
-    ierr = DMPlexPointLocalRead(dmFace,face,facegeom,&fg);CHKERRQ(ierr);
-    ierr = DMPlexPointLocalRead(dmCell,cells[0],cellgeom,&cgL);CHKERRQ(ierr);
-    ierr = DMPlexPointLocalRead(dmCell,cells[1],cellgeom,&cgR);CHKERRQ(ierr);
-    ierr = DMPlexPointLocalRead(dm,cells[0],x,&xL);CHKERRQ(ierr);
-    ierr = DMPlexPointLocalRead(dm,cells[1],x,&xR);CHKERRQ(ierr);
-    ierr = DMPlexPointGlobalRef(dm,cells[0],f,&fL);CHKERRQ(ierr);
-    ierr = DMPlexPointGlobalRef(dm,cells[1],f,&fR);CHKERRQ(ierr);
-    ierr = (*user->model->alpha)(user->model,time,fg->centroid,&alpha,user->model->physics);CHKERRQ(ierr); /* store this? */
-    WaxpyD(-1,cgL->centroid,cgR->centroid,v); /* dx vector */
-    dx = Norm2D(v);
-    area = Norm2D(fg->normal);
-    flux = alpha*(xR[0]-xL[0])/dx;
-    if (fL) fL[0] -= area*flux / cgL->volume;
-    if (fR) fR[0] += area*flux / cgR->volume;
-  }
-  ierr = VecRestoreArrayRead(user->facegeom,&facegeom);CHKERRQ(ierr);
-  ierr = VecRestoreArrayRead(user->cellgeom,&cellgeom);CHKERRQ(ierr);
-  ierr = VecRestoreArrayRead(locX,&x);CHKERRQ(ierr);
-  ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
-  PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "ComputeFunctionLap"
-/*
- Finite difference operator for Jacobian.  The operator in SNES form.
- */
-PetscErrorCode ComputeFunctionLap(SNES snes,Vec U,Vec FU,void *ctx)
-{
-  User           user = (User)ctx;
-  DM             dm, dmFace, dmCell;
-  Vec            locU;
-  PetscErrorCode ierr;
-
-  PetscFunctionBeginUser;
-  ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr);
-  ierr = VecGetDM(user->facegeom,&dmFace);CHKERRQ(ierr);
-  ierr = VecGetDM(user->cellgeom,&dmCell);CHKERRQ(ierr);
-  ierr = DMGetLocalVector(dm,&locU);CHKERRQ(ierr);
-  ierr = DMGlobalToLocalBegin(dm, U, INSERT_VALUES, locU);CHKERRQ(ierr); /* get ghost values - echange */
-  ierr = DMGlobalToLocalEnd(dm, U, INSERT_VALUES, locU);CHKERRQ(ierr);
-
-  ierr = ApplyBC(dm, 0., locU, user);CHKERRQ(ierr); /* set ghost values - boundary */
-
-  ierr = VecZeroEntries(FU);CHKERRQ(ierr);
-  ierr = ApplyLaplacianLocal(dm,dmFace,dmCell,0.,locU,FU,user);CHKERRQ(ierr);
-
-  ierr = DMRestoreLocalVector(dm,&locU);CHKERRQ(ierr);
-
-  PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "RHSFunctionLap"
-/*
- Laplacian + Source term for TS RHS: Au - b
- */
-PetscErrorCode RHSFunctionLap(TS ts,PetscReal t,Vec globalin,Vec globalout,void *ctx)
-{
-  PetscErrorCode ierr;
-  Mat            A;
-  User           user = (User)ctx;
-  Model          mod = user->model;
-  PetscInt       cStart, cEnd, c;
-  DM             dm,dmCell;
-  const PetscScalar *cgeom_v;
-  PetscScalar *x_v,source[1];
-
-  PetscFunctionBeginUser;
-  ierr = TSGetRHSJacobian(ts,&A,NULL,NULL,&ctx);CHKERRQ(ierr);
-  ierr = MatMult(A,globalin,globalout);CHKERRQ(ierr);
-  /* subtract off RHS to solve Ax=b */
-  ierr = TSGetDM(ts,&dm);CHKERRQ(ierr);
-  ierr = VecGetDM(user->cellgeom, &dmCell);CHKERRQ(ierr);
-  ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
-  ierr = VecGetArrayRead(user->cellgeom, &cgeom_v);CHKERRQ(ierr);
-  ierr = VecGetArray(globalout, &x_v);CHKERRQ(ierr);
-  for (c = cStart; c < user->cEndInterior; ++c) {
-    CellGeom *cg;
-    PetscScalar *xc;
-    ierr = DMPlexPointLocalRead(dmCell,c,cgeom_v,&cg);CHKERRQ(ierr);
-    ierr = (*mod->source)(mod,t,cg->centroid,source,mod->physics);CHKERRQ(ierr);
-    ierr = DMPlexPointGlobalRef(dm,c,x_v,&xc);CHKERRQ(ierr);
-    if (xc) *xc -= source[0];
-  }
-  ierr = VecRestoreArrayRead(user->cellgeom, &cgeom_v);CHKERRQ(ierr);
-  ierr = VecRestoreArray(globalout, &x_v);CHKERRQ(ierr);
-  VecScale(globalout, -1.); /* why is this needed to get SNES residual correct? */
-  PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "main"
-int main(int argc, char **argv)
-{
-  MPI_Comm          comm;
-  User              user;
-  Model             mod;
-  Physics           phys;
-  DM                dm, dmDist;
-  PetscReal         ftime,dt;
-  PetscInt          dim, overlap, nsteps, nrefine,i,ilev;
-  int               CPU_word_size = 0, IO_word_size = 0, exoid;
-  float             version;
-  TS                ts;
-  TSConvergedReason reason;
-  Vec               X;
-  Mat               A;
-  PetscViewer       viewer;
-  PetscMPIInt       rank;
-  char              filename[PETSC_MAX_PATH_LEN] = "./square-4x4.exo";
-  PetscBool         vtkCellGeom;
-  PetscErrorCode    ierr;
-
-  ierr = PetscInitialize(&argc, &argv, (char*) 0, help);CHKERRQ(ierr);
-  comm = PETSC_COMM_WORLD;
-  ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
-
-  ierr = PetscNew(struct _n_User,&user);CHKERRQ(ierr);
-  ierr = PetscNew(struct _n_Model,&user->model);CHKERRQ(ierr);
-  ierr = PetscNew(struct _n_Physics,&user->model->physics);CHKERRQ(ierr);
-  mod  = user->model;
-  phys = mod->physics;
-  mod->comm = comm;
-
-  ierr = PetscOptionsBegin(comm,NULL,"Unstructured Finite Volume Options","");CHKERRQ(ierr);
-  {
-    ierr              = PetscOptionsString("-f","Exodus.II filename to read","",filename,filename,sizeof(filename),NULL);CHKERRQ(ierr);
-    user->vtkInterval = 1;
-    ierr = PetscOptionsInt("-ufv_vtk_interval","VTK output interval (0 to disable)","",user->vtkInterval,&user->vtkInterval,NULL);CHKERRQ(ierr);
-    overlap = 1;
-    ierr = PetscOptionsInt("-ufv_mesh_overlap","Number of cells to overlap partitions","",overlap,&overlap,NULL);CHKERRQ(ierr);
-    vtkCellGeom = PETSC_FALSE;
-    nrefine=0;
-    ierr = PetscOptionsInt("-nrefine","Number of refinment steps","",nrefine,&nrefine,NULL);CHKERRQ(ierr);
-    if (nrefine>=32) nrefine = 31;
-    ierr = PetscOptionsBool("-ufv_vtk_cellgeom","Write cell geometry (for debugging)","",vtkCellGeom,&vtkCellGeom,NULL);CHKERRQ(ierr);
-    ierr = PetscMemzero(phys,sizeof(struct _n_Physics));CHKERRQ(ierr);
-    ierr = PhysicsCreate_Lap(mod,phys);CHKERRQ(ierr);
-    /* Count number of fields and dofs */
-    for (phys->nfields=0,phys->dof=0; phys->field_desc[phys->nfields].name; phys->nfields++) phys->dof += phys->field_desc[phys->nfields].dof;
-
-    if (phys->dof != 1) SETERRQ1(comm,PETSC_ERR_ARG_WRONGSTATE,"Physics dof=%d",phys->dof);
-    ierr = ModelFunctionalSetFromOptions(mod);CHKERRQ(ierr);
-  }
-  ierr = PetscOptionsEnd();CHKERRQ(ierr);
-
-  ierr = SetUpBoundaries(comm, user);CHKERRQ(ierr);
-
-  /* outer congerance loop */
-  for (ilev=0;ilev<=nrefine;ilev++) {
-    user->ilev = ilev;
-    /* create base mesh */
-    if (!rank) {
-      exoid = ex_open(filename, EX_READ, &CPU_word_size, &IO_word_size, &version);
-      if (exoid <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"ex_open(\"%s\",...) did not return a valid file ID",filename);
-    } else exoid = -1;                 /* Not used */
-    ierr = DMPlexCreateExodus(comm, exoid, PETSC_TRUE, &dm);CHKERRQ(ierr);
-    if (!rank) {ierr = ex_close(exoid);CHKERRQ(ierr);}
-
-    /* refine mesh */
-    /* ierr = DMPlexSetRefinementUniform(dm, PETSC_TRUE); */
-    for (i=0;i<ilev;i++) {
-      DM refinedMesh     = NULL;
-      /* Refine mesh using ... */
-      ierr = DMRefine(dm, comm, &refinedMesh);CHKERRQ(ierr);
-      if (refinedMesh) {
-        DMLabel  faceSets;
-        PetscInt fStart, fEnd;
-        ierr = DMDestroy(&dm);CHKERRQ(ierr);
-        ierr = DMPlexGetHeightStratum(refinedMesh, 1, &fStart, &fEnd);CHKERRQ(ierr);
-        ierr = DMPlexGetLabel(refinedMesh, "Face Sets", &faceSets);CHKERRQ(ierr);
-        if (faceSets) {ierr = DMLabelFilter(faceSets, fStart, fEnd);CHKERRQ(ierr);}
-        dm  = refinedMesh;
-      }
-    }
-
-    /* Distribute mesh */
-    ierr = DMPlexDistribute(dm, "metis", overlap, NULL, &dmDist);CHKERRQ(ierr);
-    if (dmDist) {
-      ierr = DMDestroy(&dm);CHKERRQ(ierr);
-      dm   = dmDist;
-    }
-    ierr = DMSetFromOptions(dm);CHKERRQ(ierr);
-
-    {
-      DM gdm;
-      ierr = DMPlexGetHeightStratum(dm, 0, NULL, &user->cEndInterior);CHKERRQ(ierr);
-      ierr = DMPlexConstructGhostCells(dm, NULL, &user->numGhostCells, &gdm);CHKERRQ(ierr);
-      ierr = DMDestroy(&dm);CHKERRQ(ierr);
-      dm   = gdm;
-    }
-    /* DMView(dm,PETSC_VIEWER_STDOUT_WORLD); */
-    ierr = ConstructGeometry(dm, &user->facegeom, &user->cellgeom, user);CHKERRQ(ierr);
-    if (0) {ierr = VecView(user->cellgeom, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);}
-    ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
-    ierr = DMPlexSetPreallocationCenterDimension(dm, 0);CHKERRQ(ierr);
-
-    /* Set up DM with section describing local vector and configure local vector. */
-    ierr = SetUpLocalSpace(dm, user);CHKERRQ(ierr);
-
-    ierr = DMCreateGlobalVector(dm, &X);CHKERRQ(ierr);
-    ierr = SetInitialCondition(dm, X, user);CHKERRQ(ierr);
-
-    if (vtkCellGeom) {
-      DM  dmCell;
-      Vec partition;
-
-      ierr = OutputVTK(dm, "ex32-cellgeom.vtk", &viewer);CHKERRQ(ierr);
-      ierr = VecView(user->cellgeom, viewer);CHKERRQ(ierr);
-      ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
-      ierr = CreatePartitionVec(dm, &dmCell, &partition);CHKERRQ(ierr);
-      ierr = OutputVTK(dmCell, "ex32-partition.vtk", &viewer);CHKERRQ(ierr);
-      ierr = VecView(partition, viewer);CHKERRQ(ierr);
-      ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
-      ierr = VecDestroy(&partition);CHKERRQ(ierr);
-      ierr = DMDestroy(&dmCell);CHKERRQ(ierr);
-    }
-
-    ierr = DMCreateMatrix(dm, &A);CHKERRQ(ierr);
-    { /* finite difference Jacobian A */
-      SNES snes;
-      MatStructure flag;
-
-      ierr = SNESCreate(comm, &snes);CHKERRQ(ierr);
-      ierr = SNESSetDM(snes,dm);CHKERRQ(ierr);
-      ierr = SNESSetFunction(snes,NULL,ComputeFunctionLap,user);CHKERRQ(ierr);
-      ierr = SNESSetJacobian(snes,A,A,SNESComputeJacobianDefaultColor,0);CHKERRQ(ierr);
-
-      ierr = SNESComputeJacobian(snes,X,&A,&A,&flag);CHKERRQ(ierr);
-      /* ierr = MatView(A,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); */
-      ierr = SNESDestroy(&snes);CHKERRQ(ierr);
-    }
-    ierr = MatGetSize(A, &user->N[ilev], NULL);CHKERRQ(ierr);
-
-    /* create TS */
-    ierr = TSCreate(comm, &ts);CHKERRQ(ierr);
-    ierr = TSSetProblemType(ts,TS_LINEAR);CHKERRQ(ierr);
-    ierr = TSSetType(ts, TSBEULER);CHKERRQ(ierr);
-    ierr = TSSetDM(ts, dm);CHKERRQ(ierr);
-    ierr = TSMonitorSet(ts,MonitorVTK,user,NULL);CHKERRQ(ierr);
-    ierr = TSSetRHSJacobian(ts,A,A,TSComputeRHSJacobianConstant,&user);CHKERRQ(ierr);
-    ierr = TSSetRHSFunction(ts,NULL,RHSFunctionLap,user);CHKERRQ(ierr);
-    ierr = TSSetDuration(ts,1,2.0);CHKERRQ(ierr);
-    dt   = 1.e20;
-    ierr = TSSetInitialTimeStep(ts,0.0,dt);CHKERRQ(ierr);
-    ierr = TSSetSolution(ts,X);CHKERRQ(ierr);
-    ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
-
-    /* solve */
-    ierr = TSSolve(ts,X);CHKERRQ(ierr);
-    ierr = TSGetSolveTime(ts,&ftime);CHKERRQ(ierr);
-    ierr = TSGetTimeStepNumber(ts,&nsteps);CHKERRQ(ierr);
-    ierr = TSGetConvergedReason(ts,&reason);CHKERRQ(ierr);
-    PetscPrintf(PETSC_COMM_WORLD,"%s at time %G after %D steps. Errors: |e|_inf=%e |e|_2=%e, N=%d\n",TSConvergedReasons[reason],ftime,nsteps,user->errorInf[ilev],user->error2[ilev],user->N[ilev]);
-
-    /* clean up */
-    ierr = TSDestroy(&ts);CHKERRQ(ierr);
-    ierr = DMDestroy(&dm);CHKERRQ(ierr);
-    ierr = VecDestroy(&user->cellgeom);CHKERRQ(ierr);
-    ierr = VecDestroy(&user->facegeom);CHKERRQ(ierr);
-    ierr = VecDestroy(&X);CHKERRQ(ierr);
-    ierr = MatDestroy(&A);CHKERRQ(ierr);
-  }
-
-  /* print out convergence rates */
-  if (nrefine>0 && mod->numCall>0) {
-    PetscPrintf(PETSC_COMM_WORLD,"\tN | Convergance rates |e|_inf\t\t|e|_2  (%d)\n",mod->numCall);
-    for (ilev=1;ilev<=nrefine;ilev++) {
-      const PetscReal log2r = 1.0/log(2.0);
-      PetscReal ratio = user->errorInf[ilev-1]/user->errorInf[ilev];
-      const PetscReal convergeRate_Inf = log(ratio)*log2r;
-      ratio = user->error2[ilev-1]/user->error2[ilev];
-      const PetscReal convergeRate_2 = log(ratio)*log2r;
-      PetscPrintf(PETSC_COMM_WORLD,"\t%d\t\t\t%e\t%e\n",user->N[ilev],convergeRate_Inf,convergeRate_2);
-    }
-  }
-
-  ierr = FunctionalLinkDestroy(&user->model->functionalRegistry);CHKERRQ(ierr);
-  ierr = BoundaryLinkDestroy(&user->model->boundary);CHKERRQ(ierr);
-  ierr = PetscFree(user->model->functionalMonitored);CHKERRQ(ierr);
-  ierr = PetscFree(user->model->functionalCall);CHKERRQ(ierr);
-  ierr = PetscFree(user->model->physics->data);CHKERRQ(ierr);
-  ierr = PetscFree(user->model->physics);CHKERRQ(ierr);
-  ierr = PetscFree(user->model);CHKERRQ(ierr);
-  ierr = PetscFree(user);CHKERRQ(ierr);
-
-  /* { */
-  /*   // setup interpolation */
-  /*   const PetscInt nref=3; */
-  /*   PetscReal targetx[nref+1],sourcex[nref+1],RCoef[nref*nref]; */
-  /*   int i,j,k; */
-  /*   for (i=0,j=-(nref-1)/2;i<=nref;i++,j++) targetx[i] = j; */
-  /*   //  normal interpolant */
-  /*   //  sources:     0             1            2           3 */
-  /*   //  targets (*)  |      o      |      *     |     o     | */
-  /*   //  coord       -1             0.           1           2 */
-  /*   for (i=0,j=-(nref-1)/2;i<=nref;i++,j++) sourcex[i] = j; */
-  /*   printf("source: "); */
-  /*   for (i=0;i<=nref;i++) printf(" %e",sourcex[i]); printf("\ntargetx:"); */
-  /*   for (i=0;i<=nref;i++) printf(" %e",targetx[i]); printf("\n"); */
-  /*   // PetscDTReconstructPoly(degree,nsource,*sourcex,ntarget,*targetx,*R) */
-  /*   ierr = PetscDTReconstructPoly(nref-1,nref,sourcex,nref,targetx,RCoef);CHKERRQ(ierr); */
-  /*   for (i=0,k=0;i<nref;i++) { */
-  /*     printf("RCoef:"); */
-  /*     for (j=0;j<nref;j++,k++) printf(" %e",RCoef[k]); */
-  /*     printf("\n"); */
-  /*   } */
-  /* } */
-
-  ierr = PetscFinalize();
-  return(0);
-  }

# src/ts/examples/tutorials/ex33.c

+static char help[] = "Second Order Finite Volume Example of Laplacian div \alpha grad u = f.\n";
+/*F
+
+  Multigrid Poisson solver:
+     Laplace(u) = f
+
+  A numerical test on page 64 of
+    A MULTIGRID TUTORIAL by Briggs, Henson & McCormick'
+  u := (x^2-x^4)*(y^4-y^2)*(z^4-z^2) with alpha(x) = -1
+
+  Course grid mesh is read in from an ExodusII file, usually generated by Cubit.
+F*/
+#include <petscts.h>
+#include <petscdmplex.h>
+#include <petscsf.h>
+#include <petscblaslapack.h>
+#if defined(PETSC_HAVE_EXODUSII)
+#include <exodusII.h>
+#else
+#error This example requires ExodusII support. Reconfigure using --download-exodusii
+#endif
+
+#define DIM 2                   /* Geometric dimension */
+#define ALEN(a) (sizeof(a)/sizeof((a)[0]))
+
+/* Represents continuum physical equations. */
+typedef struct _n_Physics *Physics;
+
+/* Physical model includes boundary conditions, initial conditions, and functionals of interest. It is
+ * discretization-independent, but its members depend on the scenario being solved. */
+typedef struct _n_Model *Model;
+
+/* 'User' implements a discretization of a continuous model. */
+typedef struct _n_User *User;
+
+typedef PetscErrorCode (*SolutionFunction)(Model,PetscReal,const PetscReal*,PetscScalar*,void*);
+typedef PetscErrorCode (*BoundaryFunction)(Model,PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*);
+typedef PetscErrorCode (*FunctionalFunction)(Model,PetscReal,const PetscReal*,const PetscScalar*,PetscReal*,void*);
+static PetscErrorCode ModelBoundaryRegister(Model,const char*,BoundaryFunction,void*,PetscInt,const PetscInt*);
+static PetscErrorCode ModelSolutionSetDefault(Model,SolutionFunction,void*);
+static PetscErrorCode ModelSourceSetDefault(Model,SolutionFunction,void*);
+static PetscErrorCode ModelAlphaSetDefault(Model,SolutionFunction,void*);
+static PetscErrorCode ModelFunctionalRegister(Model,const char*,PetscInt*,FunctionalFunction,void*);
+static PetscErrorCode OutputVTK(DM,const char*,PetscViewer*);
+
+struct FieldDescription {
+  const char *name;
+  PetscInt dof;
+};
+
+typedef struct _n_BoundaryLink *BoundaryLink;
+struct _n_BoundaryLink {
+  char             *name;
+  BoundaryFunction func;
+  void             *ctx;
+  PetscInt         numids;
+  PetscInt         *ids;
+  BoundaryLink     next;
+};
+
+typedef struct _n_FunctionalLink *FunctionalLink;
+struct _n_FunctionalLink {
+  char               *name;
+  FunctionalFunction func;
+  void               *ctx;
+  PetscInt           offset;
+  FunctionalLink     next;
+};
+
+struct _n_Physics {
+  PetscInt        dof;          /* number of degrees of freedom per cell */
+  void            *data;
+  PetscInt        nfields;
+  const struct FieldDescription *field_desc;
+};
+
+struct _n_Model {
+  MPI_Comm         comm;        /* Does not do collective communicaton, but some error conditions can be collective */
+  Physics          physics;
+  BoundaryLink     boundary;
+  FunctionalLink   functionalRegistry;
+  PetscInt         maxComputed;
+  PetscInt         numMonitored;
+  FunctionalLink   *functionalMonitored;
+  PetscInt         numCall;
+  FunctionalLink   *functionalCall;
+  SolutionFunction solution;
+  SolutionFunction source;
+  SolutionFunction alpha;
+  void             *solutionctx;
+};
+
+struct _n_User {
+  PetscInt       numGhostCells;
+  PetscInt       cEndInterior; /* First boundary ghost cell */
+  Vec            cellgeom, facegeom;
+  PetscInt       vtkInterval;   /* For monitor */
+  Model          model;
+  PetscReal      errorInf[32], error2[32];
+  PetscInt       N[32],ilev;
+};
+
+typedef struct {
+  PetscScalar normal[DIM];              /* Area-scaled normals */
+  PetscScalar centroid[DIM];            /* Location of centroid (quadrature point) */
+} FaceGeom;
+
+typedef struct {
+  PetscScalar centroid[DIM];
+  PetscScalar volume;
+} CellGeom;
+
+PETSC_STATIC_INLINE PetscScalar DotD(const PetscScalar *x, const PetscScalar *y) {PetscScalar sum = 0.0; PetscInt d; for (d = 0; d < DIM; ++d) sum += x[d]*y[d]; return sum;}
+PETSC_STATIC_INLINE PetscReal Norm2D(const PetscScalar *x) {return PetscSqrtReal(PetscAbsScalar(DotD(x,x)));}
+PETSC_STATIC_INLINE void axD(const PetscScalar a,PetscScalar *x){PetscInt i; for (i=0; i<DIM; i++) x[i] *= a;}
+PETSC_STATIC_INLINE void waxD(const PetscScalar a,const PetscScalar *x, PetscScalar *w){PetscInt i; for (i=0; i<DIM; i++) w[i] = x[i]*a;}
+PETSC_STATIC_INLINE void NormalizeD(PetscScalar *x) {PetscReal a = 1./Norm2D(x); PetscInt d; for (d = 0; d < DIM; ++d) x[d] *= a;}
+PETSC_STATIC_INLINE void WaxpyD(PetscScalar a, const PetscScalar *x, const PetscScalar *y, PetscScalar *w) {PetscInt d; for (d = 0; d < DIM; ++d) w[d] = a*x[d] + y[d];}
+
+/******************* Laplacian ********************/
+typedef struct {
+  struct {
+    PetscInt Error;
+  } functional;
+} Physics_Lap;
+
+static const struct FieldDescription PhysicsFields_Lap[] = {{"U",1},{NULL,0}};
+
+#undef __FUNCT__
+#define __FUNCT__ "PhysicsBoundary_DiriHomo"
+static PetscErrorCode PhysicsBoundary_DiriHomo(Model mod, PetscReal time, const PetscReal *fCentroid, const PetscReal *normal, const PetscScalar *xInterior, PetscScalar *xGhost, void *ctx)
+{
+  /* Physics      phys = (Physics)ctx; */
+  /* Physics_Lap *lap = (Physics_Lap*)phys->data; */
+  PetscFunctionBeginUser;
+  xGhost[0] = -xInterior[0];
+  PetscFunctionReturn(0);
+}
+
+/* U: (x^2-x^4)*(y^4-y^2)*(z^4-z^2) */
+#undef __FUNCT__
+#define __FUNCT__ "PhysicsSolution_Lap_x4_x2"
+static PetscErrorCode PhysicsSolution_Lap_x4_x2(Model mod,PetscReal time,const PetscReal *x,PetscScalar *u,void *ctx)
+{
+  /* Physics        phys    = (Physics)ctx; */
+  /* Physics_Lap *lap = (Physics_Lap*)phys->data; */
+  PetscScalar r;
+  PetscScalar e[DIM];
+  PetscInt d;
+
+  PetscFunctionBeginUser;
+  for (d = 0; d < DIM; ++d) e[d] = x[d]*x[d];
+  for (d = 0; d < DIM; ++d) e[d] *= 1 - e[d];
+  for (r = 1.0, d = 0; d < DIM; ++d) r *= e[d];
+  u[0] = r;
+  PetscFunctionReturn(0);
+}
+
+/* RHS: Lap((x^2-x^4)*(y^4-y^2)*(z^4-z^2)) */
+#undef __FUNCT__
+#define __FUNCT__ "PhysicsSource_Lap_x4_x2"
+static PetscErrorCode PhysicsSource_Lap_x4_x2(Model mod,PetscReal time,const PetscReal *x,PetscScalar *rhs,void *ctx)
+{
+  /* Physics        phys    = (Physics)ctx; */
+  /* Physics_Lap *lap = (Physics_Lap*)phys->data; */
+  PetscScalar r;
+  PetscScalar x2[DIM]; /* = x*x; */
+  PetscScalar a[DIM]; /* = x2*(1-x2) */
+  PetscScalar b[DIM]; /* = 2-12*x2; */
+  PetscInt d;
+
+  PetscFunctionBeginUser;
+  for (d = 0; d < DIM; ++d) x2[d] = x[d]*x[d];
+  for (d = 0; d < DIM; ++d) a[d] = x2[d]*(1-x2[d]);
+  for (d = 0; d < DIM; ++d) b[d] = 2-12*x2[d];
+
+  /* cross products of coordinates prevents the use of D_TERM */
+  if (DIM==2)
+    r = b[0]*a[1] + a[0]*b[1];
+  else if (DIM==3)
+    r = b[0]*a[1]*a[2] + a[0]*b[1]*a[2] + a[0]*a[1]*b[2];
+  else  if (DIM==1)
+   r = b[0];
+  else
+    SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for solution this DIM %d",DIM);
+  rhs[0] = -r; /* keep solution positive */
+  PetscFunctionReturn(0);
+}
+#undef __FUNCT__
+#define __FUNCT__ "PhysicsAlpha_n1"
+static PetscErrorCode PhysicsAlpha_n1(Model mod,PetscReal time,const PetscReal *x,PetscScalar *alpha,void *ctx)
+{
+  /* Physics        phys    = (Physics)ctx; */
+  /* Physics_Lap *lap = (Physics_Lap*)phys->data; */
+  PetscFunctionBeginUser;
+  alpha[0] = -1.0;
+  PetscFunctionReturn(0);
+}
+#undef __FUNCT__
+#define __FUNCT__ "PhysicsFunctional_Error"
+static PetscErrorCode PhysicsFunctional_Error(Model mod,PetscReal time,const PetscScalar *x,const PetscScalar *u,PetscReal *f,void *ctx)
+{
+  Physics        phys    = (Physics)ctx;
+  Physics_Lap   *lap = (Physics_Lap*)phys->data;
+  PetscScalar    yexact[1];
+  PetscErrorCode ierr;
+
+  PetscFunctionBeginUser;
+  ierr = (*mod->solution)(mod,time,x,yexact,phys);CHKERRQ(ierr);
+  /* f[lap->functional.Error] = PetscAbsScalar(u[0]-yexact[0]); */
+  f[lap->functional.Error] = PetscRealPart(u[0]-yexact[0]); /* lets do real error for plotting */
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "PhysicsCreate_Lap"
+static PetscErrorCode PhysicsCreate_Lap(Model mod,Physics phys)
+{
+  Physics_Lap *lap;
+  PetscErrorCode ierr;
+
+  PetscFunctionBeginUser;
+  phys->field_desc = PhysicsFields_Lap;
+  ierr = PetscNew(Physics_Lap,&phys->data);CHKERRQ(ierr);
+  lap = (Physics_Lap*)phys->data;
+
+  {
+    const PetscInt diriids[] = {100,101,102,103};
+    /* Register "canned" boundary conditions and defaults for where to apply. */
+    ierr = ModelBoundaryRegister(mod,"diri",PhysicsBoundary_DiriHomo,phys,ALEN(diriids),diriids);CHKERRQ(ierr);
+    /* Initial/transient solution with default boundary conditions */
+    ierr = ModelSolutionSetDefault(mod,PhysicsSolution_Lap_x4_x2,phys);CHKERRQ(ierr);
+    ierr = ModelSourceSetDefault(mod,PhysicsSource_Lap_x4_x2,phys);CHKERRQ(ierr);
+    ierr = ModelAlphaSetDefault(mod,PhysicsAlpha_n1,phys);CHKERRQ(ierr);
+    /* Register "canned" functionals */
+    ierr = ModelFunctionalRegister(mod,"error",&lap->functional.Error,PhysicsFunctional_Error,phys);CHKERRQ(ierr);
+  }
+  PetscFunctionReturn(0);
+}
+
+/* #undef __FUNCT__ */
+/* #define __FUNCT__ "IsExteriorOrGhostFace" */
+/* static PetscErrorCode IsExteriorOrGhostFace(DM dm,PetscInt face,PetscBool *isghost) */
+/* { */
+/*   PetscErrorCode ierr; */
+/*   PetscInt       ghost,boundary; */
+
+/*   PetscFunctionBeginUser; */
+/*   *isghost = PETSC_FALSE; */
+/*   ierr = DMPlexGetLabelValue(dm, "ghost", face, &ghost);CHKERRQ(ierr); */
+/*   ierr = DMPlexGetLabelValue(dm, "Face Sets", face, &boundary);CHKERRQ(ierr); */
+/*   if (ghost >= 0 || boundary >= 0) *isghost = PETSC_TRUE; */
+/*   PetscFunctionReturn(0); */
+/* } */
+
+#undef __FUNCT__
+#define __FUNCT__ "ConstructGeometry"
+/* Set up face data and cell data */
+PetscErrorCode ConstructGeometry(DM dm, Vec *facegeom, Vec *cellgeom, User user)
+{
+  DM             dmFace, dmCell;
+  DMLabel        ghostLabel;
+  PetscSection   sectionFace, sectionCell;
+  PetscSection   coordSection;
+  Vec            coordinates;
+  PetscScalar    *fgeom, *cgeom;
+  PetscInt       dim, cStart, cEnd, c, fStart, fEnd, f;
+  PetscErrorCode ierr;
+
+  PetscFunctionBeginUser;
+  ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
+  if (dim != DIM) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No support for dim %D != DIM %D",dim,DIM);
+  ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
+  ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
+
+  /* Make cell centroids and volumes */
+  ierr = DMClone(dm, &dmCell);CHKERRQ(ierr);
+  ierr = DMSetCoordinateSection(dmCell, coordSection);CHKERRQ(ierr);
+  ierr = DMSetCoordinatesLocal(dmCell, coordinates);CHKERRQ(ierr);
+  ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &sectionCell);CHKERRQ(ierr);
+  ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
+  ierr = PetscSectionSetChart(sectionCell, cStart, cEnd);CHKERRQ(ierr);
+  for (c = cStart; c < cEnd; ++c) {
+    ierr = PetscSectionSetDof(sectionCell, c, sizeof(CellGeom)/sizeof(PetscScalar));CHKERRQ(ierr);
+  }
+  ierr = PetscSectionSetUp(sectionCell);CHKERRQ(ierr);
+  ierr = DMSetDefaultSection(dmCell, sectionCell);CHKERRQ(ierr);
+  ierr = PetscSectionDestroy(&sectionCell);CHKERRQ(ierr); /* relinquish our reference */
+
+  ierr = DMCreateLocalVector(dmCell, cellgeom);CHKERRQ(ierr);
+  ierr = VecGetArray(*cellgeom, &cgeom);CHKERRQ(ierr);
+  for (c = cStart; c < user->cEndInterior; ++c) {
+    CellGeom *cg;
+
+    ierr = DMPlexPointLocalRef(dmCell, c, cgeom, &cg);CHKERRQ(ierr);
+    ierr = PetscMemzero(cg,sizeof(*cg));CHKERRQ(ierr);
+    ierr = DMPlexComputeCellGeometryFVM(dmCell, c, &cg->volume, cg->centroid, NULL);CHKERRQ(ierr);
+  }
+  /* Compute face normals, centroid(used?) */
+  ierr = DMClone(dm, &dmFace);CHKERRQ(ierr);
+  ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &sectionFace);CHKERRQ(ierr);
+  ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr);
+  ierr = PetscSectionSetChart(sectionFace, fStart, fEnd);CHKERRQ(ierr);
+
+  for (f = fStart; f < fEnd; ++f) {
+    ierr = PetscSectionSetDof(sectionFace, f, sizeof(FaceGeom)/sizeof(PetscScalar));CHKERRQ(ierr);
+  }
+  ierr = PetscSectionSetUp(sectionFace);CHKERRQ(ierr);
+  ierr = DMSetDefaultSection(dmFace, sectionFace);CHKERRQ(ierr);
+  ierr = PetscSectionDestroy(&sectionFace);CHKERRQ(ierr);
+  ierr = DMCreateLocalVector(dmFace, facegeom);CHKERRQ(ierr);
+  ierr = VecGetArray(*facegeom, &fgeom);CHKERRQ(ierr);
+  ierr = DMPlexGetLabel(dm, "ghost", &ghostLabel);CHKERRQ(ierr);
+  for (f = fStart; f < fEnd; ++f) {
+    FaceGeom *fg;
+    PetscInt  ghost;
+
+    ierr = DMLabelGetValue(ghostLabel, f, &ghost);CHKERRQ(ierr);
+    if (ghost >= 0) continue;
+    ierr = DMPlexPointLocalRef(dmFace, f, fgeom, &fg);CHKERRQ(ierr);
+    ierr = DMPlexComputeCellGeometryFVM(dm, f, NULL, fg->centroid, fg->normal);CHKERRQ(ierr);
+    /* Flip face orientation if necessary to match ordering in support, and Update minimum radius */
+    {
+      CellGeom       *cL, *cR;
+      const PetscInt *cells;
+      PetscReal      *lcentroid, *rcentroid;
+      PetscScalar     v[3];
+      PetscInt        d;
+
+      ierr = DMPlexGetSupport(dm, f, &cells);CHKERRQ(ierr);
+      ierr = DMPlexPointLocalRead(dmCell, cells[0], cgeom, &cL);CHKERRQ(ierr);
+      ierr = DMPlexPointLocalRead(dmCell, cells[1], cgeom, &cR);CHKERRQ(ierr);
+      lcentroid = cells[0] >= user->cEndInterior ? fg->centroid : cL->centroid;
+      rcentroid = cells[1] >= user->cEndInterior ? fg->centroid : cR->centroid;
+      WaxpyD(-1, lcentroid, rcentroid, v);
+      if (DotD(fg->normal, v) < 0) {
+        for (d = 0; d < dim; ++d) fg->normal[d] = -fg->normal[d];
+      }
+      if (DotD(fg->normal, v) <= 0) {
+#if DIM == 2
+        SETERRQ5(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Direction for face %d could not be fixed, normal (%g,%g) v (%g,%g)", f, fg->normal[0], fg->normal[1], v[0], v[1]);
+#elif DIM == 3
+        SETERRQ7(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Direction for face %d could not be fixed, normal (%g,%g,%g) v (%g,%g,%g)", f, fg->normal[0], fg->normal[1], fg->normal[2], v[0], v[1], v[2]);
+#else
+#  error DIM not supported
+#endif
+      }
+    }
+  }
+  /* Compute centroids of ghost cells */
+  for (c = user->cEndInterior; c < cEnd; ++c) {
+    FaceGeom       *fg;
+    const PetscInt *cone,    *support;
+    PetscInt        coneSize, supportSize, s;
+
+    ierr = DMPlexGetConeSize(dmCell, c, &coneSize);CHKERRQ(ierr);
+    if (coneSize != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ghost cell %d has cone size %d != 1", c, coneSize);
+    ierr = DMPlexGetCone(dmCell, c, &cone);CHKERRQ(ierr);
+    ierr = DMPlexGetSupportSize(dmCell, cone[0], &supportSize);CHKERRQ(ierr);
+    if (supportSize != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %d has support size %d != 1", cone[0], supportSize);
+    ierr = DMPlexGetSupport(dmCell, cone[0], &support);CHKERRQ(ierr);
+    ierr = DMPlexPointLocalRef(dmFace, cone[0], fgeom, &fg);CHKERRQ(ierr);
+    for (s = 0; s < 2; ++s) {
+      /* Reflect ghost centroid across plane of face */
+      if (support[s] == c) {
+        const CellGeom *ci;
+        CellGeom       *cg;
+        PetscScalar     c2f[3], a;
+
+        ierr = DMPlexPointLocalRead(dmCell, support[(s+1)%2], cgeom, &ci);CHKERRQ(ierr);
+        WaxpyD(-1, ci->centroid, fg->centroid, c2f); /* cell to face centroid */
+        a    = DotD(c2f, fg->normal)/DotD(fg->normal, fg->normal);
+        ierr = DMPlexPointLocalRef(dmCell, support[s], cgeom, &cg);CHKERRQ(ierr);
+        WaxpyD(2*a, fg->normal, ci->centroid, cg->centroid);
+        cg->volume = ci->volume;
+      }
+    }
+  }
+  ierr = VecRestoreArray(*facegeom, &fgeom);CHKERRQ(ierr);
+  ierr = VecRestoreArray(*cellgeom, &cgeom);CHKERRQ(ierr);
+  ierr = DMDestroy(&dmCell);CHKERRQ(ierr);
+  ierr = DMDestroy(&dmFace);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "SetUpLocalSpace"
+PetscErrorCode SetUpLocalSpace(DM dm, User user)
+{
+  PetscSection   stateSection;
+  Physics        phys;
+  PetscInt       dof = user->model->physics->dof, *cind, d, stateSize, cStart, cEnd, c, i;
+  PetscErrorCode ierr;
+
+  PetscFunctionBeginUser;
+  ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
+  ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &stateSection);CHKERRQ(ierr);
+  phys = user->model->physics;
+  ierr = PetscSectionSetNumFields(stateSection,phys->nfields);CHKERRQ(ierr);
+  for (i=0; i<phys->nfields; i++) {
+    ierr = PetscSectionSetFieldName(stateSection,i,phys->field_desc[i].name);CHKERRQ(ierr);
+    ierr = PetscSectionSetFieldComponents(stateSection,i,phys->field_desc[i].dof);CHKERRQ(ierr);
+  }
+  ierr = PetscSectionSetChart(stateSection, cStart, cEnd);CHKERRQ(ierr);
+
+  for (c = cStart; c < cEnd; ++c) {
+    for (i=0; i<phys->nfields; i++) {
+      ierr = PetscSectionSetFieldDof(stateSection,c,i,phys->field_desc[i].dof);CHKERRQ(ierr);
+    }
+    ierr = PetscSectionSetDof(stateSection, c, dof);CHKERRQ(ierr);
+  }
+  for (c = user->cEndInterior; c < cEnd; ++c) {
+    ierr = PetscSectionSetConstraintDof(stateSection, c, dof);CHKERRQ(ierr);
+  }
+  ierr = PetscSectionSetUp(stateSection);CHKERRQ(ierr);
+  ierr = PetscMalloc(dof * sizeof(PetscInt), &cind);CHKERRQ(ierr);
+  for (d = 0; d < dof; ++d) cind[d] = d;
+  for (c = user->cEndInterior; c < cEnd; ++c) {
+    ierr = PetscSectionSetConstraintIndices(stateSection, c, cind);CHKERRQ(ierr);
+  }
+  ierr = PetscFree(cind);CHKERRQ(ierr);
+  ierr = PetscSectionGetStorageSize(stateSection, &stateSize);CHKERRQ(ierr);
+  ierr = DMSetDefaultSection(dm,stateSection);CHKERRQ(ierr);
+  ierr = PetscSectionDestroy(&stateSection);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "SetUpBoundaries"
+PetscErrorCode SetUpBoundaries(MPI_Comm comm, User user)
+{
+  Model          mod = user->model;
+  PetscErrorCode ierr;
+  BoundaryLink   b;
+
+  PetscFunctionBeginUser;
+  ierr = PetscOptionsBegin(comm,NULL,"Boundary condition options","");CHKERRQ(ierr);
+  for (b = mod->boundary; b; b=b->next) {
+    char      optname[512];
+    PetscInt  ids[512],len = 512;
+    PetscBool flg;
+    ierr = PetscSNPrintf(optname,sizeof optname,"-bc_%s",b->name);CHKERRQ(ierr);
+    ierr = PetscMemzero(ids,sizeof(ids));CHKERRQ(ierr);
+    ierr = PetscOptionsIntArray(optname,"List of boundary IDs","",ids,&len,&flg);CHKERRQ(ierr);
+    if (flg) {
+      /* TODO: check all IDs to make sure they exist in the mesh */
+      ierr      = PetscFree(b->ids);CHKERRQ(ierr);
+      b->numids = len;
+      ierr      = PetscMalloc(len*sizeof(PetscInt),&b->ids);CHKERRQ(ierr);
+      ierr      = PetscMemcpy(b->ids,ids,len*sizeof(PetscInt));CHKERRQ(ierr);
+    }
+  }
+  ierr = PetscOptionsEnd();CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "ModelBoundaryRegister"
+/* The ids are just defaults, can be overridden at command line. I expect to set defaults based on names in the future. */
+static PetscErrorCode ModelBoundaryRegister(Model mod,const char *name,BoundaryFunction bcFunc,void *ctx,PetscInt numids,const PetscInt *ids)
+{
+  PetscErrorCode ierr;
+  BoundaryLink   link;
+
+  PetscFunctionBeginUser;
+  ierr          = PetscNew(struct _n_BoundaryLink,&link);CHKERRQ(ierr);
+  ierr          = PetscStrallocpy(name,&link->name);CHKERRQ(ierr);
+  link->numids  = numids;
+  ierr          = PetscMalloc(numids*sizeof(PetscInt),&link->ids);CHKERRQ(ierr);
+  ierr          = PetscMemcpy(link->ids,ids,numids*sizeof(PetscInt));CHKERRQ(ierr);
+  link->func    = bcFunc;
+  link->ctx     = ctx;
+  link->next    = mod->boundary;
+  mod->boundary = link;
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "BoundaryLinkDestroy"
+static PetscErrorCode BoundaryLinkDestroy(BoundaryLink *link)
+{
+  PetscErrorCode ierr;
+  BoundaryLink   l,next;
+
+  PetscFunctionBeginUser;
+  if (!link) PetscFunctionReturn(0);
+  l     = *link;
+  *link = NULL;
+  for (; l; l=next) {
+    next = l->next;
+    ierr = PetscFree(l->ids);CHKERRQ(ierr);
+    ierr = PetscFree(l->name);CHKERRQ(ierr);
+    ierr = PetscFree(l);CHKERRQ(ierr);
+  }
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "ModelBoundaryFind"
+static PetscErrorCode ModelBoundaryFind(Model mod,PetscInt id,BoundaryFunction *bcFunc,void **ctx)
+{
+  BoundaryLink link;
+  PetscInt     i;
+
+  PetscFunctionBeginUser;
+  *bcFunc = NULL;
+  for (link=mod->boundary; link; link=link->next) {
+    for (i=0; i<link->numids; i++) {
+      if (link->ids[i] == id) {
+        *bcFunc = link->func;
+        *ctx    = link->ctx;
+        PetscFunctionReturn(0);
+      }
+    }
+  }
+  SETERRQ1(mod->comm,PETSC_ERR_USER,"Boundary ID %D not associated with any registered boundary condition",id);
+  PetscFunctionReturn(0);
+}
+#undef __FUNCT__
+#define __FUNCT__ "ModelSolutionSetDefault"
+/* Behavior will be different for multi-physics or when using non-default boundary conditions */
+static PetscErrorCode ModelSolutionSetDefault(Model mod,SolutionFunction func,void *ctx)
+{
+  PetscFunctionBeginUser;
+  mod->solution    = func;
+  mod->solutionctx = ctx;
+  PetscFunctionReturn(0);
+}
+#undef __FUNCT__
+#define __FUNCT__ "ModelSourceSetDefault"
+/* Behavior will be different for multi-physics or when using non-default boundary conditions */
+static PetscErrorCode ModelSourceSetDefault(Model mod,SolutionFunction func,void *ctx)
+{
+  PetscFunctionBeginUser;
+  mod->source    = func;
+  /* mod->solutionctx = ctx; */
+  PetscFunctionReturn(0);
+}
+#undef __FUNCT__
+#define __FUNCT__ "ModelAlphaSetDefault"
+/* Behavior will be different for multi-physics or when using non-default boundary conditions */
+static PetscErrorCode ModelAlphaSetDefault(Model mod,SolutionFunction func,void *ctx)
+{
+  PetscFunctionBeginUser;
+  mod->alpha    = func;
+  /* mod->solutionctx = ctx; */
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "ModelFunctionalRegister"
+static PetscErrorCode ModelFunctionalRegister(Model mod,const char *name,PetscInt *offset,FunctionalFunction func,void *ctx)
+{
+  PetscErrorCode ierr;
+  FunctionalLink link,*ptr;
+  PetscInt       lastoffset = -1;
+
+  PetscFunctionBeginUser;
+  for (ptr=&mod->functionalRegistry; *ptr; ptr = &(*ptr)->next) lastoffset = (*ptr)->offset;
+  ierr         = PetscNew(struct _n_FunctionalLink,&link);CHKERRQ(ierr);
+  ierr         = PetscStrallocpy(name,&link->name);CHKERRQ(ierr);
+  link->offset = lastoffset + 1;
+  link->func   = func;
+  link->ctx    = ctx;
+  link->next   = NULL;
+  *ptr         = link;
+  *offset      = link->offset;
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "ModelFunctionalSetFromOptions"
+static PetscErrorCode ModelFunctionalSetFromOptions(Model mod)
+{
+  PetscErrorCode ierr;
+  PetscInt       i,j;
+  FunctionalLink link;
+  char           *names[256];
+
+  PetscFunctionBeginUser;
+  mod->numMonitored = ALEN(names);
+  ierr = PetscOptionsStringArray("-monitor","list of functionals to monitor","",names,&mod->numMonitored,NULL);CHKERRQ(ierr);
+  /* Create list of functionals that will be computed somehow */
+  ierr = PetscMalloc(mod->numMonitored*sizeof(FunctionalLink),&mod->functionalMonitored);CHKERRQ(ierr);
+  /* Create index of calls that we will have to make to compute these functionals (over-allocation in general). */
+  ierr = PetscMalloc(mod->numMonitored*sizeof(FunctionalLink),&mod->functionalCall);CHKERRQ(ierr);
+  mod->numCall = 0;
+
+  for (i=0; i<mod->numMonitored; i++) {
+    for (link=mod->functionalRegistry; link; link=link->next) {
+      PetscBool match;
+      ierr = PetscStrcasecmp(names[i],link->name,&match);CHKERRQ(ierr);
+      if (match) break;
+    }
+    if (!link) SETERRQ1(mod->comm,PETSC_ERR_USER,"No known functional '%s'",names[i]);
+    mod->functionalMonitored[i] = link;
+    for (j=0; j<i; j++) {
+      if (mod->functionalCall[j]->func == link->func && mod->functionalCall[j]->ctx == link->ctx) goto next_name;
+    }
+    mod->functionalCall[mod->numCall++] = link; /* Just points to the first link using the result. There may be more results. */
+next_name:
+    ierr = PetscFree(names[i]);CHKERRQ(ierr);
+  }
+
+  /* Find out the maximum index of any functional computed by a function we will be calling (even if we are not using it) */
+  mod->maxComputed = -1;
+  for (link=mod->functionalRegistry; link; link=link->next) {
+    for (i=0; i<mod->numCall; i++) {
+      FunctionalLink call = mod->functionalCall[i];
+      if (link->func == call->func && link->ctx == call->ctx) {
+        mod->maxComputed = PetscMax(mod->maxComputed,link->offset);
+      }
+    }
+  }
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "FunctionalLinkDestroy"
+static PetscErrorCode FunctionalLinkDestroy(FunctionalLink *link)
+{
+  PetscErrorCode ierr;
+  FunctionalLink l,next;
+
+  PetscFunctionBeginUser;
+  if (!link) PetscFunctionReturn(0);
+  l     = *link;
+  *link = NULL;
+  for (; l; l=next) {
+    next = l->next;
+    ierr = PetscFree(l->name);CHKERRQ(ierr);
+    ierr = PetscFree(l);CHKERRQ(ierr);
+  }
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "SetInitialCondition"
+PetscErrorCode SetInitialCondition(DM dm, Vec X, User user)
+{
+  PetscScalar       *x;
+  PetscInt          cStart, cEnd, cEndInterior = user->cEndInterior, c;
+  PetscErrorCode    ierr;
+
+  PetscFunctionBeginUser;
+  ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
+  ierr = VecGetArray(X, &x);CHKERRQ(ierr);
+  for (c = cStart; c < cEndInterior; ++c) {
+    PetscScalar    *xc;
+    ierr = DMPlexPointGlobalRef(dm,c,x,&xc);CHKERRQ(ierr);
+    if (xc) *xc = 0.; /* process ghost */
+  }
+  ierr = VecRestoreArray(X, &x);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "OutputVTK"
+static PetscErrorCode OutputVTK(DM dm, const char *filename, PetscViewer *viewer)
+{
+  PetscErrorCode ierr;
+
+  PetscFunctionBeginUser;
+  ierr = PetscViewerCreate(PetscObjectComm((PetscObject)dm), viewer);CHKERRQ(ierr);
+  ierr = PetscViewerSetType(*viewer, PETSCVIEWERVTK);CHKERRQ(ierr);
+  ierr = PetscViewerFileSetName(*viewer, filename);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "MonitorVTK"
+static PetscErrorCode MonitorVTK(TS ts,PetscInt stepnum,PetscReal time,Vec X,void *ctx)
+{
+  User           user = (User)ctx;
+  DM             dm;
+  PetscViewer    viewer;
+  char           filename[PETSC_MAX_PATH_LEN],*ftable = NULL;
+  PetscReal      xnorm;
+  Vec            E=0;
+  PetscErrorCode ierr;
+
+  PetscFunctionBeginUser;
+  ierr = PetscObjectSetName((PetscObject) X, "solution");CHKERRQ(ierr);
+  ierr = VecGetDM(X,&dm);CHKERRQ(ierr);
+  ierr = VecNorm(X,NORM_INFINITY,&xnorm);CHKERRQ(ierr);
+  ierr = VecDuplicate(X,&E);CHKERRQ(ierr);
+  if (stepnum >= 0) {           /* No summary for final time */
+    Model             mod = user->model;
+    PetscInt          c,cStart,cEnd,fcount,i;
+    size_t            ftableused,ftablealloc;
+    const PetscScalar *cellgeom_v,*x_v;
+    PetscScalar       *func_v;
+    DM                dmCell;
+    PetscReal         *fmin,*fmax,*fintegral,*ftmp;
+    fcount = mod->maxComputed+1;
+    ierr   = PetscMalloc4(fcount,PetscReal,&fmin,fcount,PetscReal,&fmax,fcount,PetscReal,&fintegral,fcount,PetscReal,&ftmp);CHKERRQ(ierr);
+    for (i=0; i<fcount; i++) {
+      fmin[i]      = PETSC_MAX_REAL;
+      fmax[i]      = PETSC_MIN_REAL;
+      fintegral[i] = 0;
+    }
+    ierr = DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);CHKERRQ(ierr);
+    ierr = VecGetDM(user->cellgeom,&dmCell);CHKERRQ(ierr);
+    ierr = VecGetArrayRead(user->cellgeom,&cellgeom_v);CHKERRQ(ierr);
+    ierr = VecGetArrayRead(X,&x_v);CHKERRQ(ierr);
+    ierr = VecGetArray(E,&func_v);CHKERRQ(ierr);
+    for (c=cStart; c<user->cEndInterior; c++) {
+      const CellGeom    *cg;
+      const PetscScalar *cx;
+      PetscScalar *cfunc;
+      ierr = DMPlexPointLocalRead(dmCell,c,cellgeom_v,&cg);CHKERRQ(ierr);
+      ierr = DMPlexPointGlobalRead(dm,c,x_v,&cx);CHKERRQ(ierr);
+      ierr = DMPlexPointGlobalRef(dm,c,func_v,&cfunc);CHKERRQ(ierr);
+      if (!cx) continue;        /* not a global cell */
+      for (i=0; i<mod->numCall; i++) {
+        FunctionalLink flink = mod->functionalCall[i];
+        ierr = (*flink->func)(mod,time,cg->centroid,cx,&ftmp[i],flink->ctx);CHKERRQ(ierr);
+        *cfunc = ftmp[i]; /* save last functional */
+        if (i==mod->numCall-1) ftmp[i] = PetscAbs(ftmp[i]); /* allready a real, this assumes we want an abs - error */
+      }
+      for (i=0; i<fcount; i++) {
+        fmin[i]       = PetscMin(fmin[i],ftmp[i]);
+        fmax[i]       = PetscMax(fmax[i],ftmp[i]);
+        fintegral[i] += cg->volume * ftmp[i];
+      }
+    }
+    ierr = VecRestoreArrayRead(user->cellgeom,&cellgeom_v);CHKERRQ(ierr);
+    ierr = VecRestoreArrayRead(X,&x_v);CHKERRQ(ierr);
+    ierr = VecRestoreArray(E,&func_v);CHKERRQ(ierr);
+    ierr = MPI_Allreduce(MPI_IN_PLACE,fmin,fcount,MPIU_REAL,MPI_MIN,PetscObjectComm((PetscObject)ts));CHKERRQ(ierr);
+    ierr = MPI_Allreduce(MPI_IN_PLACE,fmax,fcount,MPIU_REAL,MPI_MAX,PetscObjectComm((PetscObject)ts));CHKERRQ(ierr);
+    ierr = MPI_Allreduce(MPI_IN_PLACE,fintegral,fcount,MPIU_REAL,MPI_SUM,PetscObjectComm((PetscObject)ts));CHKERRQ(ierr);
+    if (fcount) { user->errorInf[user->ilev] = fmax[mod->numCall-1]; user->error2[user->ilev] = fintegral[mod->numCall-1]; } /* keep for convergence test */
+    ftablealloc = fcount * 100;
+    ftableused  = 0;
+    ierr        = PetscMalloc(ftablealloc,&ftable);CHKERRQ(ierr);
+    for (i=0; i<mod->numMonitored; i++) {
+      size_t         countused;
+      char           buffer[256],*p;
+      FunctionalLink flink = mod->functionalMonitored[i];
+      PetscInt       id    = flink->offset;
+      if (i % 3) {
+        ierr = PetscMemcpy(buffer,"  ",2);CHKERRQ(ierr);
+        p    = buffer + 2;
+      } else if (i) {
+        char newline[] = "\n";
+        ierr = PetscMemcpy(buffer,newline,sizeof newline-1);CHKERRQ(ierr);
+        p    = buffer + sizeof newline - 1;
+      } else {
+        p = buffer;
+      }
+      ierr = PetscSNPrintfCount(p,sizeof buffer-(p-buffer),"%12s functional: [%10.7G,%10.7G] integral=%10.7G",&countused,flink->name,fmin[id],fmax[id],fintegral[id]);CHKERRQ(ierr);
+      countused += p - buffer;
+      if (countused > ftablealloc-ftableused-1) { /* reallocate */
+        char *ftablenew;
+        ftablealloc = 2*ftablealloc + countused;
+        ierr = PetscMalloc(ftablealloc,&ftablenew);CHKERRQ(ierr);
+        ierr = PetscMemcpy(ftablenew,ftable,ftableused);CHKERRQ(ierr);
+        ierr = PetscFree(ftable);CHKERRQ(ierr);
+        ftable = ftablenew;
+      }
+      ierr = PetscMemcpy(ftable+ftableused,buffer,countused);CHKERRQ(ierr);
+      ftableused += countused;
+      ftable[ftableused] = 0;
+    }
+    ierr = PetscFree4(fmin,fmax,fintegral,ftmp);CHKERRQ(ierr);
+
+    ierr = PetscPrintf(PetscObjectComm((PetscObject)ts),"% 3D  time %8.4G  |x| %8.4G  %s\n",stepnum,time,xnorm,ftable ? ftable : "");CHKERRQ(ierr);
+    ierr = PetscFree(ftable);CHKERRQ(ierr);
+  }
+  if (user->vtkInterval < 1) {
+    ierr = VecDestroy(&E);CHKERRQ(ierr);
+    PetscFunctionReturn(0);
+  }
+  if (user->model->functionalCall && ((stepnum == -1) ^ (stepnum % user->vtkInterval == 0))) {
+    Model          mod   = user->model;
+    FunctionalLink flink = mod->functionalCall[mod->numCall-1];
+    if (stepnum == -1) {        /* Final time is not multiple of normal time interval, write it anyway */
+      ierr = TSGetTimeStepNumber(ts,&stepnum);CHKERRQ(ierr);
+    }
+    ierr = PetscSNPrintf(filename,sizeof filename,"ex32-%03D-u.vtu",stepnum);CHKERRQ(ierr);
+    ierr = OutputVTK(dm,filename,&viewer);CHKERRQ(ierr);
+    ierr = VecView(X,viewer);CHKERRQ(ierr);
+    ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
+    ierr = PetscObjectSetName((PetscObject) E, flink->name);CHKERRQ(ierr);
+    ierr = PetscSNPrintf(filename,sizeof filename,"ex32-%03D-%s.vtu",stepnum,flink->name);CHKERRQ(ierr);
+    ierr = OutputVTK(dm,filename,&viewer);CHKERRQ(ierr);
+    ierr = VecView(E,viewer);CHKERRQ(ierr);
+    ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
+  }