Commits

Tobin Isaac committed 701e0d5

dmplex, petscfe: proposed interface for nonconforming elements

Comments (0)

Files changed (7)

include/petsc-private/dmpleximpl.h

   PetscInt             dim;               /* Topological mesh dimension */
 
   /* Sieve */
-  PetscSection         coneSection;       /* Layout of cones (inedges for DAG) */
-  PetscInt             maxConeSize;       /* Cached for fast lookup */
-  PetscInt            *cones;             /* Cone for each point */
-  PetscInt            *coneOrientations;  /* Orientation of each cone point, means cone traveral should start on point 'o', and if negative start on -(o+1) and go in reverse */
-  PetscSection         supportSection;    /* Layout of cones (inedges for DAG) */
-  PetscInt             maxSupportSize;    /* Cached for fast lookup */
-  PetscInt            *supports;          /* Cone for each point */
-  PetscBool            refinementUniform; /* Flag for uniform cell refinement */
-  PetscReal            refinementLimit;   /* Maximum volume for refined cell */
-  PetscInt             hybridPointMax[8]; /* Allow segregation of some points, each dimension has a divider (used in VTK output and refinement) */
+  PetscSection         coneSection;        /* Layout of cones (inedges for DAG) */
+  PetscInt             maxConeSize;        /* Cached for fast lookup */
+  PetscInt            *cones;              /* Cone for each point */
+  PetscInt            *coneOrientations;   /* Orientation of each cone point, means cone traveral should start on point 'o', and if negative start on -(o+1) and go in reverse */
+  PetscInt            *coneInterpolations; /* Interpolation of cone data, for nonconforming topologies */
+  PetscSection         supportSection;     /* Layout of cones (inedges for DAG) */
+  PetscInt             maxSupportSize;     /* Cached for fast lookup */
+  PetscInt            *supports;           /* Cone for each point */
+  PetscBool            refinementUniform;  /* Flag for uniform cell refinement */
+  PetscReal            refinementLimit;    /* Maximum volume for refined cell */
+  PetscInt             hybridPointMax[8];  /* Allow segregation of some points, each dimension has a divider (used in VTK output and refinement) */
 
   PetscInt            *facesTmp;          /* Work space for faces operation */
 

include/petsc-private/petscfeimpl.h

                                       void (*)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]),
                                       void (*)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]),
                                       PetscScalar[]);
+  PetscErrorCode (*interpolateclosure)(PetscFE,PetscInt,const PetscInt[],PetscScalar[]);
+  PetscErrorCode (*interpolateclosuretranspose)(PetscFE,PetscInt,const PetscInt[],PetscScalar[]);
 };
 
 struct _p_PetscFE {

include/petscdmplex.h

 PETSC_EXTERN PetscErrorCode DMPlexInsertConeOrientation(DM, PetscInt, PetscInt, PetscInt);
 PETSC_EXTERN PetscErrorCode DMPlexGetConeOrientation(DM, PetscInt, const PetscInt *[]);
 PETSC_EXTERN PetscErrorCode DMPlexSetConeOrientation(DM, PetscInt, const PetscInt[]);
+PETSC_EXTERN PetscErrorCode DMPlexGetConeInterpolation(DM, PetscInt, const PetscInt *[]);
+PETSC_EXTERN PetscErrorCode DMPlexSetConeInterpolation(DM, PetscInt, const PetscInt[]);
 PETSC_EXTERN PetscErrorCode DMPlexGetSupportSize(DM, PetscInt, PetscInt *);
 PETSC_EXTERN PetscErrorCode DMPlexSetSupportSize(DM, PetscInt, PetscInt);
 PETSC_EXTERN PetscErrorCode DMPlexGetSupport(DM, PetscInt, const PetscInt *[]);
 
 PETSC_EXTERN PetscErrorCode DMPlexComputeCellGeometry(DM, PetscInt, PetscReal *, PetscReal *, PetscReal *, PetscReal *);
 PETSC_EXTERN PetscErrorCode DMPlexVecGetClosure(DM, PetscSection, Vec, PetscInt, PetscInt *, PetscScalar *[]);
+PETSC_EXTERN PetscErrorCode DMPlexVecGetClosureInterpolated(DM, PetscFE [], PetscSection, Vec, PetscInt, PetscInt *, PetscScalar *[]);
 PETSC_EXTERN PetscErrorCode DMPlexVecRestoreClosure(DM, PetscSection, Vec, PetscInt, PetscInt *, PetscScalar *[]);
 PETSC_EXTERN PetscErrorCode DMPlexVecSetClosure(DM, PetscSection, Vec, PetscInt, const PetscScalar[], InsertMode);
+PETSC_EXTERN PetscErrorCode DMPlexVecSetClosureInterpolated(DM, PetscFE [],PetscSection, Vec, PetscInt, const PetscScalar[], InsertMode);
 PETSC_EXTERN PetscErrorCode DMPlexMatSetClosure(DM, PetscSection, PetscSection, Mat, PetscInt, const PetscScalar[], InsertMode);
 PETSC_EXTERN PetscErrorCode DMPlexCreateClosureIndex(DM, PetscSection);
 
 PETSC_EXTERN PetscErrorCode PetscFEOpenCLSetRealType(PetscFE, PetscDataType);
 PETSC_EXTERN PetscErrorCode PetscFEOpenCLGetRealType(PetscFE, PetscDataType *);
 
+PETSC_EXTERN PetscErrorCode PetscFEInterpolateClosure(PetscFE, PetscInt, const PetscInt [], PetscScalar[]);
+PETSC_EXTERN PetscErrorCode PetscFEInterpolateClosureTranspose(PetscFE, PetscInt, const PetscInt [], PetscScalar[]);
+
 /* TODO: Should be moved inside DM */
 typedef struct {
   PetscFE         *fe;
   void (**f1BdFuncs)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]); /* The f_1 functions for each field */
 } PetscFEM;
 
+/*E
+     PetscFEInterpolationRange - a range that is the subset of a point (cell, facet, etc.) to which a PetscSpace
+                                 defined on the point can be interpolated
+
+   Level: advanced
+
+   A PetscSpace whose domain is a point in a mesh, such as a cell or a facet, may need a representation on a subset of
+   that point, during solution/residual transfer or because the mesh has nonconforming interfaces.
+
+   Interpolation ranges:
+.   FE_INTERPOLATE_FULL - The range of interpolation is the domain.
+.   FE_INTERPOLATE_HALFx - The range is a half of the domain.
+$
+$    If the domain is a simplex
+$     FE_INTERPOLATE_HALF{2x,2x+1} are the first and the second half of the simplex after bisection along edge x.
+$
+$    If the domain is a hypercube
+$     FE_INTERPOLATE_HALF{2x,2x+1} are the first and second half of the triangle after bisection in direction x.
+$
+.   FE_INTERPOLATE_QUARTERx - The range is a quarter of the domain.
+$
+$    If the domain is a triangle, the quarters are formed from connecting the midpoints of the edges:
+$     FE_INTERPOLATE_QUARTER{0,1,2} are the quarters touching the vertices,
+$     FE_INTERPOLATE_QUARTER3 is the quarter in the middle.
+$
+$    If the domain is a quadrilateral, the quarters are formed from bisecting in both directions:
+$     FE_INTERPOLATE_QUARTERx is the quadter touching vertex x.
+$
+$    If the domain is a hexahedron,
+$     FE_INTERPOLATE_QUARTER{4x,4x+1,4x+2,4x+3} are the quarters formed from bisection in every direction except
+$      direction x.
+$
+.   FE_INTERPOLATE_EIGHTHx - The range is an eighth of the domain.
+$
+$    If the domain is a tetrahedron, four of the eighths each connect a vertex to the midpoints of its adjacent edges,
+$     and four of the eighths quadrisect the remaining octahedron in the middle.
+$     FE_INTERPOLATE_EIGHTH{0,1,2,3} are the eighths touching the vertices,
+$     FE_INTERPOLATE_EIGHTH{4,5,6,7} are the eighths in the middle
+$
+$    If the domain is a hexahedron, the eights are formed from bisection in all directions,
+$     FE_INTERPOLATE_EIGHTHx is the eighth touching vertex x.
+$
+
+.seealso: PetscFEInterpolateClosure(), DMPlexGetConeInterpolation(), DMPlexSetConeInterpolation()
+E*/
+typedef enum{
+             FE_INTERPOLATE_FULL = 0,
+             FE_INTERPOLATE_HALF0,
+             FE_INTERPOLATE_HALF1,
+             FE_INTERPOLATE_HALF2,
+             FE_INTERPOLATE_HALF3,
+             FE_INTERPOLATE_HALF4,
+             FE_INTERPOLATE_HALF5,
+             FE_INTERPOLATE_HALF6,
+             FE_INTERPOLATE_HALF7,
+             FE_INTERPOLATE_HALF8,
+             FE_INTERPOLATE_HALF9,
+             FE_INTERPOLATE_HALF10,
+             FE_INTERPOLATE_HALF11,
+             FE_INTERPOLATE_QUARTER0,
+             FE_INTERPOLATE_QUARTER1,
+             FE_INTERPOLATE_QUARTER2,
+             FE_INTERPOLATE_QUARTER3,
+             FE_INTERPOLATE_QUARTER4,
+             FE_INTERPOLATE_QUARTER5,
+             FE_INTERPOLATE_QUARTER6,
+             FE_INTERPOLATE_QUARTER7,
+             FE_INTERPOLATE_QUARTER8,
+             FE_INTERPOLATE_QUARTER9,
+             FE_INTERPOLATE_QUARTER10,
+             FE_INTERPOLATE_QUARTER11,
+             FE_INTERPOLATE_EIGHTH0,
+             FE_INTERPOLATE_EIGHTH1,
+             FE_INTERPOLATE_EIGHTH2,
+             FE_INTERPOLATE_EIGHTH3,
+             FE_INTERPOLATE_EIGHTH4,
+             FE_INTERPOLATE_EIGHTH5,
+             FE_INTERPOLATE_EIGHTH6,
+             FE_INTERPOLATE_EIGHTH7,
+             FE_INTERPOLATE_NUM_PREDEFINED} PetscFEInterpolationRange;
+
+
 #endif

src/dm/dt/interface/dtfe.c

 }
 
 #undef __FUNCT__
+#define __FUNCT__ "PetscFEInterpolateClosure"
+/*C
+  PetscFEInterpolateClosure - Interpolate values from nonconforming basis functions to the basis functions of a finite element.
+
+  Not collective
+
+  Input Parameters:
++ fe           = The PetscFE object for the field being interpolated
+. numPoints    - The number of points from the reference element's closure
+. interp       - A PetscFEInterpolationRange for each point, which indicates the type of nonconforming basis involved
+- values       - the values to be interpolated
+
+  Output Parameter
+. values       - the values of the field after interpolation
+
+.seealso: PetscFEInterpolationRange, PetscFEInterpolateClosureTranspose()
+*/
+PetscErrorCode PetscFEInterpolateClosure (PetscFE fe, PetscInt numPoints, const PetscInt interp[], PetscScalar values[])
+{
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(fe, PETSCFE_CLASSID, 1);
+  PetscValidPointer(interp, 3);
+  PetscValidPointer(values, 4);
+  if (fe->ops->interpolateclosure) {ierr = (*fe->ops->interpolateclosure)(fe, numPoints, interp, values);CHKERRQ(ierr);}
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "PetscFEInterpolateClosureTranspose"
+/*C
+  PetscFEInterpolateClosureTranspose - Apply the transpose of interpolation from nonconforming basis functions to the dual basis functions of a finite element.
+
+  Not collective
+
+  Input Parameters:
++ fe           = The PetscFE object for the field being interpolated
+. numPoints    - The number of points from the reference element's closure
+. interp       - A PetscFEInterpolationRange for each point, which indicates the type of nonconforming basis involved
+- values       - the values to be interpolated
+
+  Output Parameter
+. values       - the values of the field after interpolation
+
+.seealso: PetscFEInterpolationRange, PetscFEInterpolateClosure()
+*/
+PetscErrorCode PetscFEInterpolateClosureTranspose (PetscFE fe, PetscInt numPoints, const PetscInt interp[], PetscScalar values[])
+{
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(fe, PETSCFE_CLASSID, 1);
+  PetscValidPointer(interp, 3);
+  PetscValidPointer(values, 4);
+  if (fe->ops->interpolateclosuretranspose) {ierr = (*fe->ops->interpolateclosuretranspose)(fe, numPoints, interp, values);CHKERRQ(ierr);}
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
 #define __FUNCT__ "PetscFECreateDefault"
 /*@
   PetscFECreateDefault - Create a PetscFE for basic FEM computation

src/dm/impls/plex/plex.c

   ierr = PetscSectionDestroy(&mesh->coneSection);CHKERRQ(ierr);
   ierr = PetscFree(mesh->cones);CHKERRQ(ierr);
   ierr = PetscFree(mesh->coneOrientations);CHKERRQ(ierr);
+  ierr = PetscFree(mesh->coneInterpolations);CHKERRQ(ierr);
   ierr = PetscSectionDestroy(&mesh->supportSection);CHKERRQ(ierr);
   ierr = PetscFree(mesh->supports);CHKERRQ(ierr);
   ierr = PetscFree(mesh->facesTmp);CHKERRQ(ierr);
 }
 
 #undef __FUNCT__
+#define __FUNCT__ "DMPlexCheckInterpolation"
+PetscErrorCode DMPlexCheckInterpolation(DM dm)
+{
+  DM_Plex       *mesh = (DM_Plex*) dm->data;
+  PetscInt       size;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  if (!mesh->coneInterpolations) {
+    ierr = PetscSectionGetStorageSize(mesh->coneSection, &size);CHKERRQ(ierr);
+    ierr = PetscCalloc1(size, &mesh->coneInterpolations);CHKERRQ(ierr);
+  }
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "DMPlexGetConeInterpolation"
+/*@C
+  DMPlexGetConeInterpolation - Return the interpolations on the in-edges for this point in the Sieve DAG
+
+  Not collective
+
+  Input Parameters:
++ mesh - The DMPlex
+- p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
+
+  Output Parameter:
+. coneInterpolation - An array of interpolations which are on the in-edges for point p. An interpolation is an integer
+                      giving the prescription for interpolation of fields on the point from fields on the cone points.
+                      If it is null, there are no interpolations, i.e. each value would be FE_INTERPOLATE_FULL,
+                      otherwise the interpretation depends on geometric object the point p represents.
+  Level: advanced
+
+  Fortran Notes:
+  Since it returns an array, this routine is only available in Fortran 90, and you must
+  include petsc.h90 in your code.
+
+  You must also call DMPlexRestoreConeInterpolation() after you finish using the returned array.
+
+.seealso: DMPlexCreate(), DMPlexGetCone(), DMPlexSetCone(), DMPlexSetChart()
+@*/
+PetscErrorCode DMPlexGetConeInterpolation(DM dm, PetscInt p, const PetscInt *coneInterpolation[])
+{
+  DM_Plex       *mesh = (DM_Plex*) dm->data;
+  PetscInt       off;
+  PetscErrorCode ierr;
+  PetscInt dof;
+
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
+  if (!mesh->coneInterpolations) {
+    *coneInterpolation = NULL;
+    PetscFunctionReturn(0);
+  }
+#if defined(PETSC_USE_DEBUG)
+  {
+    ierr = PetscSectionGetDof(mesh->coneSection, p, &dof);CHKERRQ(ierr);
+    if (dof) PetscValidPointer(coneInterpolation, 3);
+  }
+#endif
+  ierr = PetscSectionGetOffset(mesh->coneSection, p, &off);CHKERRQ(ierr);
+
+  *coneInterpolation = &mesh->coneInterpolations[off];
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "DMPlexSetConeInterpolation"
+/*@
+  DMPlexSetConeInterpolation - Set the interpolations on the in-edges for this point in the Sieve DAG
+
+  Not collective
+
+  Input Parameters:
++ mesh - The DMPlex
+. p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
+- coneInterpolation - An array of interpolations which are on the in-edges for point p. An interpolation is an integer
+                      giving the prescription for interpolation of fields on the point from fields on the cone points.
+                      If it is 0, there is no interpolation, i.e. identity, otherwise the interpretation depends on
+                      the cone size of the cone point.
+
+  Output Parameter:
+
+  Note:
+  This should be called after all calls to DMPlexSetConeSize() and DMSetUp().
+
+  Level: advanced
+
+.seealso: DMPlexCreate(), DMPlexGetConeInterpolation(), DMPlexSetCone(), DMPlexSetChart(), DMPlexSetConeSize(), DMSetUp()
+@*/
+PetscErrorCode DMPlexSetConeInterpolation(DM dm, PetscInt p, const PetscInt coneInterpolation[])
+{
+  DM_Plex       *mesh = (DM_Plex*) dm->data;
+  PetscInt       pStart, pEnd;
+  PetscInt       dof, off, c;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
+  ierr = PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);CHKERRQ(ierr);
+  ierr = PetscSectionGetDof(mesh->coneSection, p, &dof);CHKERRQ(ierr);
+  if (dof) PetscValidPointer(coneInterpolation, 3);
+  ierr = PetscSectionGetOffset(mesh->coneSection, p, &off);CHKERRQ(ierr);
+  if ((p < pStart) || (p >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Mesh point %D is not in the valid range [%D, %D)", p, pStart, pEnd);
+  ierr = DMPlexCheckInterpolation(dm);CHKERRQ(ierr);
+  for (c = 0; c < dof; ++c) {
+    PetscInt cdof, i = coneInterpolation[c];
+
+    ierr = PetscSectionGetDof(mesh->coneSection, mesh->cones[off+c], &cdof);CHKERRQ(ierr);
+    mesh->coneInterpolations[off+c] = i;
+  }
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
 #define __FUNCT__ "DMPlexInsertCone"
 PetscErrorCode DMPlexInsertCone(DM dm, PetscInt p, PetscInt conePos, PetscInt conePoint)
 {
   PetscFunctionReturn(0);
 }
 
+PetscErrorCode DMPlexGetTransitiveClosure_Internal(DM,PetscInt,PetscInt,PetscBool,PetscInt*,PetscInt**,PetscInt**);
+
 #undef __FUNCT__
 #define __FUNCT__ "DMPlexGetTransitiveClosure"
 /*@C
 @*/
 PetscErrorCode DMPlexGetTransitiveClosure(DM dm, PetscInt p, PetscBool useCone, PetscInt *numPoints, PetscInt *points[])
 {
-  DM_Plex        *mesh = (DM_Plex*) dm->data;
-  PetscInt       *closure, *fifo;
-  const PetscInt *tmp = NULL, *tmpO = NULL;
-  PetscInt        tmpSize, t;
-  PetscInt        depth       = 0, maxSize;
-  PetscInt        closureSize = 2, fifoSize = 0, fifoStart = 0;
-  PetscErrorCode  ierr;
+  return DMPlexGetTransitiveClosure_Internal(dm,p,0,useCone,numPoints,points,NULL);
+}
 
-  PetscFunctionBegin;
-  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
-  ierr    = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr);
-  /* This is only 1-level */
-  if (useCone) {
-    ierr = DMPlexGetConeSize(dm, p, &tmpSize);CHKERRQ(ierr);
-    ierr = DMPlexGetCone(dm, p, &tmp);CHKERRQ(ierr);
-    ierr = DMPlexGetConeOrientation(dm, p, &tmpO);CHKERRQ(ierr);
-  } else {
-    ierr = DMPlexGetSupportSize(dm, p, &tmpSize);CHKERRQ(ierr);
-    ierr = DMPlexGetSupport(dm, p, &tmp);CHKERRQ(ierr);
-  }
-  if (depth == 1) {
-    if (*points) {
-      closure = *points;
-    } else {
-      maxSize = 2*(PetscMax(mesh->maxConeSize, mesh->maxSupportSize)+1);
-      ierr = DMGetWorkArray(dm, maxSize, PETSC_INT, &closure);CHKERRQ(ierr);
-    }
-    closure[0] = p; closure[1] = 0;
-    for (t = 0; t < tmpSize; ++t, closureSize += 2) {
-      closure[closureSize]   = tmp[t];
-      closure[closureSize+1] = tmpO ? tmpO[t] : 0;
-    }
-    if (numPoints) *numPoints = closureSize/2;
-    if (points)    *points    = closure;
-    PetscFunctionReturn(0);
-  }
-  maxSize = 2*PetscMax(PetscMax(PetscPowInt(mesh->maxConeSize,depth+1),PetscPowInt(mesh->maxSupportSize,depth+1)),depth+1);
-  ierr    = DMGetWorkArray(dm, maxSize, PETSC_INT, &fifo);CHKERRQ(ierr);
-  if (*points) {
-    closure = *points;
-  } else {
-    ierr = DMGetWorkArray(dm, maxSize, PETSC_INT, &closure);CHKERRQ(ierr);
-  }
-  closure[0] = p; closure[1] = 0;
-  for (t = 0; t < tmpSize; ++t, closureSize += 2, fifoSize += 2) {
-    const PetscInt cp = tmp[t];
-    const PetscInt co = tmpO ? tmpO[t] : 0;
+#undef __FUNCT__
+#define __FUNCT__ "DMPlexGetTransitiveClosureInterpolation"
+/*@C
+  DMPlexGetTransitiveClosureInterpolation - Return the points on the transitive closure of the in-edges or out-edges for this point in the Sieve DAG,
+                                            with interpolation information in addition to orientation informaition
 
-    closure[closureSize]   = cp;
-    closure[closureSize+1] = co;
-    fifo[fifoSize]         = cp;
-    fifo[fifoSize+1]       = co;
-  }
-  /* Should kick out early when depth is reached, rather than checking all vertices for empty cones */
-  while (fifoSize - fifoStart) {
-    const PetscInt q   = fifo[fifoStart];
-    const PetscInt o   = fifo[fifoStart+1];
-    const PetscInt rev = o >= 0 ? 0 : 1;
-    const PetscInt off = rev ? -(o+1) : o;
+  Not collective
 
-    if (useCone) {
-      ierr = DMPlexGetConeSize(dm, q, &tmpSize);CHKERRQ(ierr);
-      ierr = DMPlexGetCone(dm, q, &tmp);CHKERRQ(ierr);
-      ierr = DMPlexGetConeOrientation(dm, q, &tmpO);CHKERRQ(ierr);
-    } else {
-      ierr = DMPlexGetSupportSize(dm, q, &tmpSize);CHKERRQ(ierr);
-      ierr = DMPlexGetSupport(dm, q, &tmp);CHKERRQ(ierr);
-      tmpO = NULL;
-    }
-    for (t = 0; t < tmpSize; ++t) {
-      const PetscInt i  = ((rev ? tmpSize-t : t) + off)%tmpSize;
-      const PetscInt cp = tmp[i];
-      /* Must propogate orientation: When we reverse orientation, we both reverse the direction of iteration and start at the other end of the chain. */
-      /* HACK: It is worse to get the size here, than to change the interpretation of -(*+1)
-       const PetscInt co = tmpO ? (rev ? -(tmpO[i]+1) : tmpO[i]) : 0; */
-      PetscInt       co = tmpO ? tmpO[i] : 0;
-      PetscInt       c;
+  Input Parameters:
++ mesh - The DMPlex
+. p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
+. useCone - PETSC_TRUE for in-edges,  otherwise use out-edges
+. points - If points is NULL on input, internal storage will be returned, otherwise the provided array is used
+- interp - If interp is NULL on input, internal storage will be returned, otherwise the provided array is used
 
-      if (rev) {
-        PetscInt childSize, coff;
-        ierr = DMPlexGetConeSize(dm, cp, &childSize);CHKERRQ(ierr);
-        coff = tmpO[i] < 0 ? -(tmpO[i]+1) : tmpO[i];
-        co   = childSize ? -(((coff+childSize-1)%childSize)+1) : 0;
-      }
-      /* Check for duplicate */
-      for (c = 0; c < closureSize; c += 2) {
-        if (closure[c] == cp) break;
-      }
-      if (c == closureSize) {
-        closure[closureSize]   = cp;
-        closure[closureSize+1] = co;
-        fifo[fifoSize]         = cp;
-        fifo[fifoSize+1]       = co;
-        closureSize           += 2;
-        fifoSize              += 2;
-      }
-    }
-    fifoStart += 2;
-  }
-  if (numPoints) *numPoints = closureSize/2;
-  if (points)    *points    = closure;
-  ierr = DMRestoreWorkArray(dm, maxSize, PETSC_INT, &fifo);CHKERRQ(ierr);
-  PetscFunctionReturn(0);
+  Output Parameters:
++ numPoints - The number of points in the closure, so points[] is of size 2*numPoints
+. points - The points and point orientations, interleaved as pairs [p0, o0, p1, o1, ...]
+- interp - The PetscFEInterpolationRange values of the points
+
+  Note:
+  If using internal storage (points is NULL on input), each call overwrites the last output.
+
+  TODO:
+  Explain how interpolations and orientations are combined
+
+  Fortran Notes:
+  Since it returns an array, this routine is only available in Fortran 90, and you must
+  include petsc.h90 in your code.
+
+  The numPoints argument is not present in the Fortran 90 binding since it is internal to the array.
+
+  Level: advanced
+
+.seealso: DMPlexRestoreTransitiveClosureInterpolation(), DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart(), DMPlexGetCone()
+@*/
+PetscErrorCode DMPlexGetTransitiveClosureInterpolation(DM dm, PetscInt p, PetscBool useCone, PetscInt *numPoints, PetscInt *points[], PetscInt *interp[])
+{
+  return DMPlexGetTransitiveClosure_Internal(dm,p,0,useCone,numPoints,points,interp);
 }
 
-#undef __FUNCT__
-#define __FUNCT__ "DMPlexGetTransitiveClosure_Internal"
 /*@C
   DMPlexGetTransitiveClosure_Internal - Return the points on the transitive closure of the in-edges or out-edges for this point in the Sieve DAG with a specified initial orientation
 
 . p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
 . orientation - The orientation of the point
 . useCone - PETSC_TRUE for in-edges,  otherwise use out-edges
-- points - If points is NULL on input, internal storage will be returned, otherwise the provided array is used
+. points - If points is NULL on input, internal storage will be returned, otherwise the provided array is used
+- interp - If interp is NULL on input, internal storage will be returned, otherwise the provided array is used
 
   Output Parameters:
 + numPoints - The number of points in the closure, so points[] is of size 2*numPoints
-- points - The points and point orientations, interleaved as pairs [p0, o0, p1, o1, ...]
+. points - The points and point orientations, interleaved as pairs [p0, o0, p1, o1, ...],
+- interp - The PetscFEInterpolationRange values of the points
 
   Note:
   If using internal storage (points is NULL on input), each call overwrites the last output.
 
 .seealso: DMPlexRestoreTransitiveClosure(), DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart(), DMPlexGetCone()
 @*/
-PetscErrorCode DMPlexGetTransitiveClosure_Internal(DM dm, PetscInt p, PetscInt ornt, PetscBool useCone, PetscInt *numPoints, PetscInt *points[])
+PetscErrorCode DMPlexGetTransitiveClosure_Internal(DM dm, PetscInt p, PetscInt ornt, PetscBool useCone, PetscInt *numPoints, PetscInt *points[], PetscInt *interp[])
 {
   DM_Plex        *mesh = (DM_Plex*) dm->data;
-  PetscInt       *closure, *fifo;
-  const PetscInt *tmp = NULL, *tmpO = NULL;
+  PetscInt       *closure, *closureinterp, *fifo, *fifointerp;
+  const PetscInt *tmp = NULL, *tmpO = NULL, *tmp1 = NULL;
   PetscInt        tmpSize, t;
   PetscInt        depth       = 0, maxSize;
   PetscInt        closureSize = 2, fifoSize = 0, fifoStart = 0;
     ierr = DMPlexGetConeSize(dm, p, &tmpSize);CHKERRQ(ierr);
     ierr = DMPlexGetCone(dm, p, &tmp);CHKERRQ(ierr);
     ierr = DMPlexGetConeOrientation(dm, p, &tmpO);CHKERRQ(ierr);
+    if (interp) {
+      ierr = DMPlexGetConeInterpolation(dm, p, &tmp1);CHKERRQ(ierr);
+    }
   } else {
     ierr = DMPlexGetSupportSize(dm, p, &tmpSize);CHKERRQ(ierr);
     ierr = DMPlexGetSupport(dm, p, &tmp);CHKERRQ(ierr);
       maxSize = 2*(PetscMax(mesh->maxConeSize, mesh->maxSupportSize)+1);
       ierr = DMGetWorkArray(dm, maxSize, PETSC_INT, &closure);CHKERRQ(ierr);
     }
+    if (interp) {
+     if (*interp) {
+       closureinterp = *interp;
+     }
+     else {
+       maxSize = (PetscMax(mesh->maxConeSize, mesh->maxSupportSize)+1);
+       ierr = DMGetWorkArray(dm, maxSize, PETSC_INT, &closureinterp);CHKERRQ(ierr);
+     }
+    }
     closure[0] = p; closure[1] = ornt;
     for (t = 0; t < tmpSize; ++t, closureSize += 2) {
       const PetscInt i = ornt >= 0 ? (t+ornt)%tmpSize : (-(ornt+1) + tmpSize-t)%tmpSize;
       closure[closureSize]   = tmp[i];
       closure[closureSize+1] = tmpO ? tmpO[i] : 0;
+      if (interp) {
+        closureinterp[closureSize/2] = tmp1 ? tmp1[i] : FE_INTERPOLATE_FULL;
+      }
     }
     if (numPoints) *numPoints = closureSize/2;
     if (points)    *points    = closure;
+    if (interp)    *interp    = closureinterp;
     PetscFunctionReturn(0);
   }
   maxSize = 2*PetscMax(PetscMax(PetscPowInt(mesh->maxConeSize,depth+1),PetscPowInt(mesh->maxSupportSize,depth+1)),depth+1);
   } else {
     ierr = DMGetWorkArray(dm, maxSize, PETSC_INT, &closure);CHKERRQ(ierr);
   }
+  if (interp) {
+    ierr = DMGetWorkArray(dm, maxSize/2, PETSC_INT, &fifointerp);CHKERRQ(ierr);
+    if (*interp) {
+      closureinterp = *interp;
+    }
+    else {
+      ierr = DMGetWorkArray(dm, maxSize/2, PETSC_INT, &closureinterp);CHKERRQ(ierr);
+    }
+  }
   closure[0] = p; closure[1] = ornt;
   for (t = 0; t < tmpSize; ++t, closureSize += 2, fifoSize += 2) {
     const PetscInt i  = ornt >= 0 ? (t+ornt)%tmpSize : (-(ornt+1) + tmpSize-t)%tmpSize;
     const PetscInt cp = tmp[i];
     PetscInt       co = tmpO ? tmpO[i] : 0;
+    PetscInt       ci = tmp1 ? tmp1[i] : FE_INTERPOLATE_FULL;
 
     if (ornt < 0) {
       PetscInt childSize, coff;
     closure[closureSize+1] = co;
     fifo[fifoSize]         = cp;
     fifo[fifoSize+1]       = co;
+    if (interp) {
+      /* TODO: apply ornt to ci */
+      closureinterp[closureSize/2]  = ci;
+      fifointerp[fifoSize/2] = ci;
+    }
   }
   /* Should kick out early when depth is reached, rather than checking all vertices for empty cones */
   while (fifoSize - fifoStart) {
     const PetscInt o   = fifo[fifoStart+1];
     const PetscInt rev = o >= 0 ? 0 : 1;
     const PetscInt off = rev ? -(o+1) : o;
+    //const PetscInt itp = interp ? fifointerp[fifoStart/2] : FE_INTERPOLATE_FULL; // turn this back on when used below
 
     if (useCone) {
       ierr = DMPlexGetConeSize(dm, q, &tmpSize);CHKERRQ(ierr);
       ierr = DMPlexGetCone(dm, q, &tmp);CHKERRQ(ierr);
       ierr = DMPlexGetConeOrientation(dm, q, &tmpO);CHKERRQ(ierr);
+      if (interp) {
+        ierr = DMPlexGetConeInterpolation(dm, q, &tmp1);CHKERRQ(ierr);
+      }
     } else {
       ierr = DMPlexGetSupportSize(dm, q, &tmpSize);CHKERRQ(ierr);
       ierr = DMPlexGetSupport(dm, q, &tmp);CHKERRQ(ierr);
       /* HACK: It is worse to get the size here, than to change the interpretation of -(*+1)
        const PetscInt co = tmpO ? (rev ? -(tmpO[i]+1) : tmpO[i]) : 0; */
       PetscInt       co = tmpO ? tmpO[i] : 0;
+      PetscInt       ci = tmp1 ? tmp1[i] : FE_INTERPOLATE_FULL;
       PetscInt       c;
 
       if (rev) {
         fifo[fifoSize+1]       = co;
         closureSize           += 2;
         fifoSize              += 2;
+        if (interp) {
+          /* TODO: apply rev, off, and itp to ci */
+          closureinterp[closureSize/2] = ci;
+          fifointerp[fifoSize/2]       = ci;
+        }
       }
     }
     fifoStart += 2;
   }
   if (numPoints) *numPoints = closureSize/2;
   if (points)    *points    = closure;
+  if (interp)    *interp    = closureinterp;
   ierr = DMRestoreWorkArray(dm, maxSize, PETSC_INT, &fifo);CHKERRQ(ierr);
+  ierr = DMRestoreWorkArray(dm, maxSize, PETSC_INT, &fifointerp);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
 
 }
 
 #undef __FUNCT__
+#define __FUNCT__ "DMPlexRestoreTransitiveClosureInterpolation"
+/*@C
+  DMPlexRestoreTransitiveClosureInterpolation - Restore the array of points on the transitive closure of the in-edges or out-edges for this point in the Sieve DAG
+
+  Not collective
+
+  Input Parameters:
++ mesh - The DMPlex
+. p - The Sieve point, which must lie in the chart set with DMPlexSetChart()
+. useCone - PETSC_TRUE for in-edges,  otherwise use out-edges
+. numPoints - The number of points in the closure, so points[] is of size 2*numPoints, zeroed on exit
+. points - The points and point orientations, interleaved as pairs [p0, o0, p1, o1, ...], zeroed on exit
+- interp - PetscFEInterpolationRange for each point in the closure
+
+  Note:
+  If not using internal storage (points is not NULL on input), this call is unnecessary
+
+  Fortran Notes:
+  Since it returns an array, this routine is only available in Fortran 90, and you must
+  include petsc.h90 in your code.
+
+  The numPoints argument is not present in the Fortran 90 binding since it is internal to the array.
+
+  Level: advanced
+
+.seealso: DMPlexGetTransitiveClosureInterpolation(), DMPlexCreate(), DMPlexSetCone(), DMPlexSetChart(), DMPlexGetCone()
+@*/
+PetscErrorCode DMPlexRestoreTransitiveClosureInterpolation(DM dm, PetscInt p, PetscBool useCone, PetscInt *numPoints, PetscInt *points[], PetscInt *interp[])
+{
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
+  if (numPoints) PetscValidIntPointer(numPoints,4);
+  if (points) PetscValidPointer(points,5);
+  ierr = DMRestoreWorkArray(dm, 0, PETSC_INT, points);CHKERRQ(ierr);
+  if (interp) PetscValidPointer(interp,6);
+  ierr = DMRestoreWorkArray(dm, 0, PETSC_INT, interp);CHKERRQ(ierr);
+  if (numPoints) *numPoints = 0;
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
 #define __FUNCT__ "DMPlexGetMaxSizes"
 /*@
   DMPlexGetMaxSizes - Return the maximum number of in-edges (cone) and out-edges (support) for any point in the Sieve DAG
   ierr = PetscSectionGetStorageSize(mesh->coneSection, &size);CHKERRQ(ierr);
   ierr = PetscMalloc1(size, &mesh->cones);CHKERRQ(ierr);
   ierr = PetscCalloc1(size, &mesh->coneOrientations);CHKERRQ(ierr);
+  mesh->coneInterpolations = NULL;
   if (mesh->maxSupportSize) {
     ierr = PetscSectionSetUp(mesh->supportSection);CHKERRQ(ierr);
     ierr = PetscSectionGetStorageSize(mesh->supportSection, &size);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
 
+#undef __FUNCT__
+#define __FUNCT__ "DMPlexGetConeInterpolations"
+PetscErrorCode DMPlexGetConeInterpolations(DM dm, PetscInt *coneInterpolations[])
+{
+  DM_Plex *mesh = (DM_Plex*) dm->data;
+
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
+  if (coneInterpolations) *coneInterpolations = mesh->coneInterpolations;
+  PetscFunctionReturn(0);
+}
+
 /******************************** FEM Support **********************************/
 
 #undef __FUNCT__
 
 #undef __FUNCT__
 #define __FUNCT__ "DMPlexVecGetClosure_Fields_Static"
-PETSC_STATIC_INLINE PetscErrorCode DMPlexVecGetClosure_Fields_Static(PetscSection section, PetscInt numPoints, const PetscInt points[], PetscInt numFields, const PetscScalar vArray[], PetscInt *size, PetscScalar array[])
+PETSC_STATIC_INLINE PetscErrorCode DMPlexVecGetClosure_Fields_Static(PetscSection section, PetscInt numPoints, const PetscInt points[], PetscInt numFields, const PetscScalar vArray[], PetscInt *size, PetscScalar array[], PetscInt offsets[])
 {
   PetscInt       offset = 0, f;
   PetscErrorCode ierr;
     PetscInt fcomp, p;
 
     ierr = PetscSectionGetFieldComponents(section, f, &fcomp);CHKERRQ(ierr);
+    if (offsets) offsets[f] = offset;
     for (p = 0; p < numPoints*2; p += 2) {
       const PetscInt point = points[p];
       const PetscInt o     = points[p+1];
 }
 
 #undef __FUNCT__
-#define __FUNCT__ "DMPlexVecGetClosure"
+#define __FUNCT__ "DMPlexVecGetClosureInterpolated"
 /*@C
-  DMPlexVecGetClosure - Get an array of the values on the closure of 'point'
+  DMPlexVecGetClosureInterpolated - Get an array of values on the closure of 'point', and interpolate them to the
+                                    finite element bases of the point
 
   Not collective
 
   Input Parameters:
 + dm - The DM
+. fe - The finite elements that represent the closure on the point of the fields in the section
 . section - The section describing the layout in v, or NULL to use the default section
 . v - The local vector
-- point - The sieve point in the DM
-
-  Output Parameters:
-+ csize - The number of values in the closure, or NULL
+. point - The sieve point in the DM
+. csize - The number of values in the closure, or NULL
 - values - The array of values, which is a borrowed array and should not be freed
 
   Fortran Notes:
 
   The csize argument is not present in the Fortran 90 binding since it is internal to the array.
 
-  Level: intermediate
+  Level: advanced
 
-.seealso DMPlexVecRestoreClosure(), DMPlexVecSetClosure(), DMPlexMatSetClosure()
+.seealso DMPlexVecGetClosure(), DMPlexVecSetClosureInterpolated(), PetscFEInterpolateClosure(), PetscFEInterpolateClosureTranspose(), PetscSpaceInterpolationRange
 @*/
-PetscErrorCode DMPlexVecGetClosure(DM dm, PetscSection section, Vec v, PetscInt point, PetscInt *csize, PetscScalar *values[])
+PetscErrorCode DMPlexVecGetClosureInterpolated(DM dm, PetscFE fe[], PetscSection section, Vec v, PetscInt point, PetscInt *csize, PetscScalar *values[])
 {
   PetscSection    clSection;
   IS              clPoints;
   PetscScalar    *array, *vArray;
   PetscInt       *points = NULL;
+  PetscInt       *interp = NULL;
+  PetscInt      **interpptr;
   const PetscInt *clp;
-  PetscInt        depth, numFields, numPoints, size;
+  PetscInt        depth, f, numFields, numPoints, size;
   PetscErrorCode  ierr;
+  PetscBool       useinterp = (values != NULL && fe != NULL);
+  PetscInt       *offsets = NULL;
 
   PetscFunctionBeginHot;
   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
   PetscValidHeaderSpecific(v, VEC_CLASSID, 3);
   ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr);
   ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr);
-  if (depth == 1 && numFields < 2) {
+  if (depth == 1 && numFields < 2 && !useinterp) {
     ierr = DMPlexVecGetClosure_Depth1_Static(dm, section, v, point, csize, values);CHKERRQ(ierr);
     PetscFunctionReturn(0);
   }
+  if (useinterp && numFields > 1) {
+    ierr = PetscMalloc(numFields,&offsets);CHKERRQ(ierr);
+  }
   /* Get points */
-  ierr = PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);CHKERRQ(ierr);
+  if (!useinterp) {
+    ierr = PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);CHKERRQ(ierr);
+    interpptr = NULL;
+  }
+  else {
+    /* clSection and clPoints don't include interpolation values */
+    clPoints = NULL;
+    interpptr = &interp;
+  }
   if (!clPoints) {
     PetscInt pStart, pEnd, p, q;
 
     ierr = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr);
-    ierr = DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr);
+    ierr = DMPlexGetTransitiveClosureInterpolation(dm, point, PETSC_TRUE, &numPoints, &points, interpptr);CHKERRQ(ierr);
     /* Compress out points not in the section */
     for (p = 0, q = 0; p < numPoints*2; p += 2) {
       if ((points[p] >= pStart) && (points[p] < pEnd)) {
         points[q*2]   = points[p];
         points[q*2+1] = points[p+1];
+        if (useinterp) {
+          interp[q] = interp[p/2];
+        }
         ++q;
       }
     }
       asize += dof;
     }
     if (!values) {
-      if (!clPoints) {ierr = DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr);}
+      if (!clPoints) {ierr = DMPlexRestoreTransitiveClosureInterpolation(dm, point, PETSC_TRUE, &numPoints, &points, interpptr);CHKERRQ(ierr);}
       else           {ierr = ISRestoreIndices(clPoints, &clp);CHKERRQ(ierr);}
       if (csize) *csize = asize;
       PetscFunctionReturn(0);
   }
   ierr = VecGetArray(v, &vArray);CHKERRQ(ierr);
   /* Get values */
-  if (numFields > 1) {ierr = DMPlexVecGetClosure_Fields_Static(section, numPoints, points, numFields, vArray, &size, array);CHKERRQ(ierr);}
+  if (numFields > 1) {ierr = DMPlexVecGetClosure_Fields_Static(section, numPoints, points, numFields, vArray, &size, array, offsets);CHKERRQ(ierr);}
   else               {ierr = DMPlexVecGetClosure_Static(section, numPoints, points, vArray, &size, array);CHKERRQ(ierr);}
+  /* interpolate the values */
+  if (useinterp) {
+    if (numFields > 1) {
+      for (f = 0; f < numFields; f++) {
+        PetscScalar *varr = &array[offsets[f]];
+
+        ierr = PetscFEInterpolateClosure(fe[f], numPoints, interp, varr);CHKERRQ(ierr);
+      }
+      ierr = PetscFree(offsets);CHKERRQ(ierr);
+    }
+    else {
+      ierr = PetscFEInterpolateClosure(fe[0], numPoints, interp, array);CHKERRQ(ierr);
+    }
+  }
   /* Cleanup points */
-  if (!clPoints) {ierr = DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr);}
+  if (!clPoints) {ierr = DMPlexRestoreTransitiveClosureInterpolation(dm, point, PETSC_TRUE, &numPoints, &points, interpptr);CHKERRQ(ierr);}
   else           {ierr = ISRestoreIndices(clPoints, &clp);CHKERRQ(ierr);}
   /* Cleanup array */
   ierr = VecRestoreArray(v, &vArray);CHKERRQ(ierr);
 }
 
 #undef __FUNCT__
+#define __FUNCT__ "DMPlexVecGetClosure"
+/*@C
+  DMPlexVecGetClosure - Get an array of the values on the closure of 'point'
+
+  Not collective
+
+  Input Parameters:
++ dm - The DM
+. section - The section describing the layout in v, or NULL to use the default section
+. v - The local vector
+- point - The sieve point in the DM
+
+  Output Parameters:
++ csize - The number of values in the closure, or NULL
+- values - The array of values, which is a borrowed array and should not be freed
+
+  Fortran Notes:
+  Since it returns an array, this routine is only available in Fortran 90, and you must
+  include petsc.h90 in your code.
+
+  The csize argument is not present in the Fortran 90 binding since it is internal to the array.
+
+  Level: intermediate
+
+.seealso DMPlexVecRestoreClosure(), DMPlexVecSetClosure(), DMPlexMatSetClosure()
+@*/
+PetscErrorCode DMPlexVecGetClosure(DM dm, PetscSection section, Vec v, PetscInt point, PetscInt *csize, PetscScalar *values[])
+{
+  return DMPlexVecGetClosureInterpolated(dm,NULL,section,v,point,csize,values);
+}
+
+#undef __FUNCT__
 #define __FUNCT__ "DMPlexVecRestoreClosure"
 /*@C
   DMPlexVecRestoreClosure - Restore the array of the values on the closure of 'point'
 }
 
 #undef __FUNCT__
-#define __FUNCT__ "DMPlexVecSetClosure"
+#define __FUNCT__ "DMPlexVecSetClosureInterpolated"
 /*@C
-  DMPlexVecSetClosure - Set an array of the values on the closure of 'point'
+  DMPlexVecSetClosureInterpolated - Set an array of the values on the closure of 'point', after interpolating them
+                                    from the finite element dual bases on the point
 
   Not collective
 
   Input Parameters:
 + dm - The DM
+. fe - The finite elements that represent the closure on the point of the fields in the section
 . section - The section describing the layout in v, or NULL to use the default section
 . v - The local vector
 . point - The sieve point in the DM
   Fortran Notes:
   This routine is only available in Fortran 90, and you must include petsc.h90 in your code.
 
-  Level: intermediate
+  Level: advanced,
 
-.seealso DMPlexVecGetClosure(), DMPlexMatSetClosure()
+.seealso DMPlexVecGetClosureInterpolated(), PetscFEInterpolateClosure(), PetscFEInterpolateClosureTranspose(), PetscSpaceInterpolationRange
 @*/
-PetscErrorCode DMPlexVecSetClosure(DM dm, PetscSection section, Vec v, PetscInt point, const PetscScalar values[], InsertMode mode)
+PetscErrorCode DMPlexVecSetClosureInterpolated(DM dm, PetscFE fe[], PetscSection section, Vec v, PetscInt point, const PetscScalar values[], InsertMode mode)
 {
   PetscSection    clSection;
   IS              clPoints;
   PetscScalar    *array;
   PetscInt       *points = NULL;
+  PetscInt       *interp = NULL;
+  PetscInt      **interpptr;
   const PetscInt *clp;
-  PetscInt        depth, numFields, numPoints, p;
+  PetscInt        depth, f, numFields, numPoints, p, numDofs;
+  PetscBool       useinterp;
+  PetscInt       *offsets, count, offset;
   PetscErrorCode  ierr;
 
   PetscFunctionBeginHot;
   if (!section) {ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);}
   PetscValidHeaderSpecific(section, PETSC_SECTION_CLASSID, 2);
   PetscValidHeaderSpecific(v, VEC_CLASSID, 3);
+  useinterp = (fe != NULL);
   ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr);
   ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr);
-  if (depth == 1 && numFields < 2 && mode == ADD_VALUES) {
+  if (depth == 1 && numFields < 2 && mode == ADD_VALUES && !useinterp) {
     ierr = DMPlexVecSetClosure_Static(dm, section, v, point, values, mode);CHKERRQ(ierr);
     PetscFunctionReturn(0);
   }
   /* Get points */
-  ierr = PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);CHKERRQ(ierr);
+  if (!useinterp) {
+    ierr = PetscSectionGetClosureIndex(section, (PetscObject) dm, &clSection, &clPoints);CHKERRQ(ierr);
+    interpptr = NULL;
+  }
+  else {
+    /* clSection and clPoints don't include interpolation */
+    clPoints = NULL;
+    interpptr = &interp;
+  }
+  if (useinterp) {
+    ierr = PetscCalloc1(numFields,&offsets);CHKERRQ(ierr);
+  }
+  numDofs = 0;
   if (!clPoints) {
-    PetscInt pStart, pEnd, q;
+    PetscInt pStart, pEnd, q, dof;
 
     ierr = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr);
-    ierr = DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr);
+    ierr = DMPlexGetTransitiveClosureInterpolation(dm, point, PETSC_TRUE, &numPoints, &points, interpptr);CHKERRQ(ierr);
     /* Compress out points not in the section */
     for (p = 0, q = 0; p < numPoints*2; p += 2) {
       if ((points[p] >= pStart) && (points[p] < pEnd)) {
         points[q*2]   = points[p];
         points[q*2+1] = points[p+1];
+        if (useinterp) {
+          interp[q] = interp[p/2];
+        }
         ++q;
       }
+      else {
+        ierr = PetscSectionGetDof(section,points[p],&dof);
+        numDofs += dof;
+        for (f = 0; f < numFields; f++) {
+          ierr = PetscSectionGetFieldDof(section,points[p],f,&dof);
+          offsets[f] += dof;
+        }
+      }
+    }
+    offset = 0;
+    for (f = 0; f < numFields; f++) {
+      count = offsets[f];
+      offsets[f] = offset;
+      offset += count;
     }
     numPoints = q;
   } else {
     numPoints = dof/2;
     points    = (PetscInt *) &clp[off];
   }
+  /* perform interpolation in place before placing in array */
+  if (useinterp) {
+    PetscScalar *wvalues;
+
+    ierr = DMGetWorkArray(dm, numDofs, PETSC_SCALAR, &wvalues);CHKERRQ(ierr);
+    ierr = PetscMemcpy(wvalues, values, numDofs * sizeof(PetscScalar));CHKERRQ(ierr);
+    for (f = 0; f < numFields; f++) {
+      PetscScalar *varr = &wvalues[offsets[f]];
+
+      ierr = PetscFEInterpolateClosureTranspose(fe[f], numPoints, interp, varr);CHKERRQ(ierr);
+    }
+    values = wvalues;
+  }
   /* Get array */
   ierr = VecGetArray(v, &array);CHKERRQ(ierr);
   /* Get values */
       SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insert mode %D", mode);
     }
   }
+  if (interp) {
+    ierr = DMRestoreWorkArray(dm, numDofs, PETSC_SCALAR, &values);CHKERRQ(ierr);
+    ierr = PetscFree(offsets);CHKERRQ(ierr);
+  }
   /* Cleanup points */
-  if (!clPoints) {ierr = DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr);}
+  if (!clPoints) {ierr = DMPlexRestoreTransitiveClosureInterpolation(dm, point, PETSC_TRUE, &numPoints, &points, interpptr);CHKERRQ(ierr);}
   else           {ierr = ISRestoreIndices(clPoints, &clp);CHKERRQ(ierr);}
   /* Cleanup array */
   ierr = VecRestoreArray(v, &array);CHKERRQ(ierr);
 }
 
 #undef __FUNCT__
+#define __FUNCT__ "DMPlexVecSetClosure"
+/*@C
+  DMPlexVecSetClosure - Set an array of the values on the closure of 'point'
+
+  Not collective
+
+  Input Parameters:
++ dm - The DM
+. section - The section describing the layout in v, or NULL to use the default section
+. v - The local vector
+. point - The sieve point in the DM
+. values - The array of values
+- mode - The insert mode, where INSERT_ALL_VALUES and ADD_ALL_VALUES also overwrite boundary conditions
+
+  Fortran Notes:
+  This routine is only available in Fortran 90, and you must include petsc.h90 in your code.
+
+  Level: intermediate
+
+.seealso DMPlexVecGetClosure(), DMPlexMatSetClosure()
+@*/
+PetscErrorCode DMPlexVecSetClosure(DM dm, PetscSection section, Vec v, PetscInt point, const PetscScalar values[], InsertMode mode)
+{
+  return DMPlexVecSetClosureInterpolated(dm, NULL, section, v, point, values, mode);
+}
+
+#undef __FUNCT__
 #define __FUNCT__ "DMPlexPrintMatSetValues"
 PetscErrorCode DMPlexPrintMatSetValues(PetscViewer viewer, Mat A, PetscInt point, PetscInt numIndices, const PetscInt indices[], const PetscScalar values[])
 {
       for (f = 0; f < numFaces; ++f) {
         PetscInt *fclosure = NULL, fclosureSize, cl, fnumCorners = 0, v;
 
-        ierr = DMPlexGetTransitiveClosure_Internal(dm, cone[f], ornt[f], PETSC_TRUE, &fclosureSize, &fclosure);CHKERRQ(ierr);
+        ierr = DMPlexGetTransitiveClosure_Internal(dm, cone[f], ornt[f], PETSC_TRUE, &fclosureSize, &fclosure, NULL);CHKERRQ(ierr);
         for (cl = 0; cl < fclosureSize*2; cl += 2) {
           const PetscInt p = fclosure[cl];
           if ((p >= vStart) && (p < vEnd)) fclosure[fnumCorners++] = p;

src/dm/impls/plex/plexfem.c

         v += numComp;
       }
     }
-    ierr = DMPlexVecSetClosure(dm, section, localX, c, values, mode);CHKERRQ(ierr);
+    ierr = DMPlexVecSetClosureInterpolated (dm, fe, section, localX, c, values, mode);CHKERRQ(ierr);
   }
   ierr = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr);
   ierr = PetscFree2(v0,J);CHKERRQ(ierr);
 
     ierr = DMPlexComputeCellGeometry(dm, c, v0, J, invJ, &detJ);CHKERRQ(ierr);
     if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
-    ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr);
+    ierr = DMPlexVecGetClosureInterpolated(dm, fe, NULL, localX, c, NULL, &x);CHKERRQ(ierr);
 
     for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) {
       void * const     ctx           = ctxs ? ctxs[field] : NULL;
 
     ierr = DMPlexComputeCellGeometry(dm, c, v0, J, invJ, &detJ);CHKERRQ(ierr);
     if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
-    ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr);
+    ierr = DMPlexVecGetClosureInterpolated(dm, fe, NULL, localX, c, NULL, &x);CHKERRQ(ierr);
 
     for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) {
       void * const     ctx           = ctxs ? ctxs[field] : NULL;
     if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, "  elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);}
     localDiff += elemDiff;
   }
-    ierr  = PetscFree7(funcVal,coords,realSpaceDer,v0,J,invJ,interpolantVec);CHKERRQ(ierr);
+  ierr  = PetscFree7(funcVal,coords,realSpaceDer,v0,J,invJ,interpolantVec);CHKERRQ(ierr);
   ierr  = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr);
   ierr  = MPI_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr);
   *diff = PetscSqrtReal(*diff);
 
     ierr = DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr);
     if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c);
-    ierr = DMPlexVecGetClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr);
+    ierr = DMPlexVecGetClosureInterpolated(dm, fe, section, X, c, NULL, &x);CHKERRQ(ierr);
     for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i];
     ierr = DMPlexVecRestoreClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr);
     if (dmAux) {
-      ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr);
+      ierr = DMPlexVecGetClosureInterpolated(dmAux, fe, sectionAux, A, c, NULL, &x);CHKERRQ(ierr);
       for (i = 0; i < cellDofAux; ++i) a[c*cellDofAux+i] = x[i];
       ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr);
     }
   }
   for (c = cStart; c < cEnd; ++c) {
     if (mesh->printFEM > 1) {ierr = DMPrintCellVector(c, name, cellDof, &elemVec[c*cellDof]);CHKERRQ(ierr);}
-    ierr = DMPlexVecSetClosure(dm, section, F, c, &elemVec[c*cellDof], ADD_VALUES);CHKERRQ(ierr);
+    ierr = DMPlexVecSetClosureInterpolated(dm, fe, section, F, c, &elemVec[c*cellDof], ADD_VALUES);CHKERRQ(ierr);
   }
   ierr = PetscFree6(u,v0,J,invJ,detJ,elemVec);CHKERRQ(ierr);
   if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);}
       ierr = DMPlexComputeCellGeometry(dm, point, &v0[f*dim], &J[f*dim*dim], &invJ[f*dim*dim], &detJ[f]);CHKERRQ(ierr);
       ierr = DMPlexComputeCellGeometryFVM(dm, point, NULL, NULL, &n[f*dim]);
       if (detJ[f] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for face %d", detJ[f], point);
-      ierr = DMPlexVecGetClosure(dm, section, X, point, NULL, &x);CHKERRQ(ierr);
+      ierr = DMPlexVecGetClosureInterpolated(dm, fe, section, X, point, NULL, &x);CHKERRQ(ierr);
       for (i = 0; i < cellDof; ++i) u[f*cellDof+i] = x[i];
       ierr = DMPlexVecRestoreClosure(dm, section, X, point, NULL, &x);CHKERRQ(ierr);
       ++f;
       ierr = DMLabelGetValue(depth, point, &dep);CHKERRQ(ierr);
       if (dep != dim-1) continue;
       if (mesh->printFEM > 1) {ierr = DMPrintCellVector(point, "BdResidual", cellDof, &elemVec[f*cellDof]);CHKERRQ(ierr);}
-      ierr = DMPlexVecSetClosure(dm, NULL, F, point, &elemVec[f*cellDof], ADD_VALUES);CHKERRQ(ierr);
+      ierr = DMPlexVecSetClosureInterpolated(dm, fe, NULL, F, point, &elemVec[f*cellDof], ADD_VALUES);CHKERRQ(ierr);
       ++f;
     }
     ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr);
 
     ierr = DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr);
     if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c);
-    ierr = DMPlexVecGetClosure(dm, NULL, jctx->u, c, NULL, &x);CHKERRQ(ierr);
+    ierr = DMPlexVecGetClosureInterpolated(dm, fe, NULL, jctx->u, c, NULL, &x);CHKERRQ(ierr);
     for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i];
     ierr = DMPlexVecRestoreClosure(dm, NULL, jctx->u, c, NULL, &x);CHKERRQ(ierr);
-    ierr = DMPlexVecGetClosure(dm, NULL, X, c, NULL, &x);CHKERRQ(ierr);
+    ierr = DMPlexVecGetClosureInterpolated(dm, fe, NULL, X, c, NULL, &x);CHKERRQ(ierr);
     for (i = 0; i < cellDof; ++i) a[c*cellDof+i] = x[i];
     ierr = DMPlexVecRestoreClosure(dm, NULL, X, c, NULL, &x);CHKERRQ(ierr);
   }
   }
   for (c = cStart; c < cEnd; ++c) {
     if (mesh->printFEM > 1) {ierr = DMPrintCellVector(c, "Jacobian Action", cellDof, &elemVec[c*cellDof]);CHKERRQ(ierr);}
-    ierr = DMPlexVecSetClosure(dm, NULL, F, c, &elemVec[c*cellDof], ADD_VALUES);CHKERRQ(ierr);
+    ierr = DMPlexVecSetClosureInterpolated(dm, fe, NULL, F, c, &elemVec[c*cellDof], ADD_VALUES);CHKERRQ(ierr);
   }
   ierr = PetscFree7(u,a,v0,J,invJ,detJ,elemVec);CHKERRQ(ierr);
   if (mesh->printFEM) {