Commits

Matt Knepley committed cb1e121

DMPlex: Moved out Preallocation and FEM support to other files

Hg-commit: b2561291e574ff9a0f8a073426e136d4a44a4dde

Comments (0)

Files changed (5)

include/petscdmplex.h

 PETSC_EXTERN PetscErrorCode DMPlexGetCoordinateSection(DM, PetscSection *);
 PETSC_EXTERN PetscErrorCode DMPlexSetCoordinateSection(DM, PetscSection);
 PETSC_EXTERN PetscErrorCode DMPlexSetPreallocationCenterDimension(DM,PetscInt);
+PETSC_EXTERN PetscErrorCode DMPlexPreallocateOperator(DM, PetscInt, PetscSection, PetscSection, PetscInt[], PetscInt[], PetscInt[], PetscInt[], Mat, PetscBool);
 PETSC_EXTERN PetscErrorCode DMPlexGetPointLocal(DM,PetscInt,PetscInt*,PetscInt*);
 PETSC_EXTERN PetscErrorCode DMPlexPointLocalRef(DM,PetscInt,PetscScalar*,void*);
 PETSC_EXTERN PetscErrorCode DMPlexPointLocalRead(DM,PetscInt,const PetscScalar*,const void*);

src/dm/impls/plex/makefile

 CPPFLAGS =
 CFLAGS   =
 FFLAGS   =
-SOURCEC  = plexcreate.c plex.c plexinterpolate.c plexgeometry.c plexlabel.c plexsubmesh.c plexexodusii.c plexvtk.c plexpoint.c plexvtu.c
+SOURCEC  = plexcreate.c plex.c plexinterpolate.c plexpreallocate.c plexgeometry.c plexlabel.c plexsubmesh.c plexexodusii.c plexvtk.c plexpoint.c plexvtu.c plexfem.c
 SOURCEF  =
 SOURCEH  =
 DIRS     = examples

src/dm/impls/plex/plex.c

 #include <petsc-private/dmpleximpl.h>   /*I      "petscdmplex.h"   I*/
 #include <../src/sys/utils/hash.h>
-#include <petsc-private/vecimpl.h>
 #include <petsc-private/isimpl.h>
 
 /* Logging support */
 }
 
 #undef __FUNCT__
-#define __FUNCT__ "DMPlexGetAdjacencySingleLevel_Private"
-PetscErrorCode DMPlexGetAdjacencySingleLevel_Private(DM dm, PetscInt p, PetscBool useClosure, const PetscInt *tmpClosure, PetscInt *adjSize, PetscInt adj[])
-{
-  const PetscInt *support = NULL;
-  PetscInt        numAdj   = 0, maxAdjSize = *adjSize, supportSize, s;
-  PetscErrorCode  ierr;
-
-  PetscFunctionBegin;
-  if (useClosure) {
-    ierr = DMPlexGetConeSize(dm, p, &supportSize);CHKERRQ(ierr);
-    ierr = DMPlexGetCone(dm, p, &support);CHKERRQ(ierr);
-    for (s = 0; s < supportSize; ++s) {
-      const PetscInt *cone = NULL;
-      PetscInt        coneSize, c, q;
-
-      ierr = DMPlexGetSupportSize(dm, support[s], &coneSize);CHKERRQ(ierr);
-      ierr = DMPlexGetSupport(dm, support[s], &cone);CHKERRQ(ierr);
-      for (c = 0; c < coneSize; ++c) {
-        for (q = 0; q < numAdj || (adj[numAdj++] = cone[c],0); ++q) {
-          if (cone[c] == adj[q]) break;
-        }
-        if (numAdj > maxAdjSize) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid mesh exceeded adjacency allocation (%D)", maxAdjSize);
-      }
-    }
-  } else {
-    ierr = DMPlexGetSupportSize(dm, p, &supportSize);CHKERRQ(ierr);
-    ierr = DMPlexGetSupport(dm, p, &support);CHKERRQ(ierr);
-    for (s = 0; s < supportSize; ++s) {
-      const PetscInt *cone = NULL;
-      PetscInt        coneSize, c, q;
-
-      ierr = DMPlexGetConeSize(dm, support[s], &coneSize);CHKERRQ(ierr);
-      ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr);
-      for (c = 0; c < coneSize; ++c) {
-        for (q = 0; q < numAdj || (adj[numAdj++] = cone[c],0); ++q) {
-          if (cone[c] == adj[q]) break;
-        }
-        if (numAdj > maxAdjSize) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid mesh exceeded adjacency allocation (%D)", maxAdjSize);
-      }
-    }
-  }
-  *adjSize = numAdj;
-  PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "DMPlexGetAdjacency_Private"
-PetscErrorCode DMPlexGetAdjacency_Private(DM dm, PetscInt p, PetscBool useClosure, const PetscInt *tmpClosure, PetscInt *adjSize, PetscInt adj[])
-{
-  const PetscInt *star  = tmpClosure;
-  PetscInt        numAdj = 0, maxAdjSize = *adjSize, starSize, s;
-  PetscErrorCode  ierr;
-
-  PetscFunctionBegin;
-  ierr = DMPlexGetTransitiveClosure(dm, p, useClosure, &starSize, (PetscInt**) &star);CHKERRQ(ierr);
-  for (s = 2; s < starSize*2; s += 2) {
-    const PetscInt *closure = NULL;
-    PetscInt        closureSize, c, q;
-
-    ierr = DMPlexGetTransitiveClosure(dm, star[s], (PetscBool)!useClosure, &closureSize, (PetscInt**) &closure);CHKERRQ(ierr);
-    for (c = 0; c < closureSize*2; c += 2) {
-      for (q = 0; q < numAdj || (adj[numAdj++] = closure[c],0); ++q) {
-        if (closure[c] == adj[q]) break;
-      }
-      if (numAdj > maxAdjSize) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid mesh exceeded adjacency allocation (%D)", maxAdjSize);
-    }
-    ierr = DMPlexRestoreTransitiveClosure(dm, star[s], (PetscBool)!useClosure, &closureSize, (PetscInt**) &closure);CHKERRQ(ierr);
-  }
-  *adjSize = numAdj;
-  PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "DMPlexSetPreallocationCenterDimension"
-PetscErrorCode DMPlexSetPreallocationCenterDimension(DM dm, PetscInt preallocCenterDim)
-{
-  DM_Plex *mesh = (DM_Plex*) dm->data;
-
-  PetscFunctionBegin;
-  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
-  mesh->preallocCenterDim = preallocCenterDim;
-  PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "DMPlexPreallocateOperator"
-PetscErrorCode DMPlexPreallocateOperator(DM dm, PetscInt bs, PetscSection section, PetscSection sectionGlobal, PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[], Mat A, PetscBool fillMatrix)
-{
-  DM_Plex           *mesh = (DM_Plex*) dm->data;
-  MPI_Comm           comm;
-  PetscSF            sf, sfDof, sfAdj;
-  PetscSection       leafSectionAdj, rootSectionAdj, sectionAdj;
-  PetscInt           nleaves, l, p;
-  const PetscInt    *leaves;
-  const PetscSFNode *remotes;
-  PetscInt           dim, pStart, pEnd, numDof, globalOffStart, globalOffEnd, numCols;
-  PetscInt          *tmpClosure, *tmpAdj, *adj, *rootAdj, *cols, *remoteOffsets;
-  PetscInt           depth, maxConeSize, maxSupportSize, maxClosureSize, maxAdjSize, adjSize;
-  PetscLayout        rLayout;
-  PetscInt           locRows, rStart, rEnd, r;
-  PetscMPIInt        size;
-  PetscBool          useClosure, debug = PETSC_FALSE;
-  PetscErrorCode     ierr;
-
-  PetscFunctionBegin;
-  ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
-  ierr = PetscOptionsGetBool(NULL, "-dm_view_preallocation", &debug, NULL);CHKERRQ(ierr);
-  ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr);
-  ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
-  ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr);
-  /* Create dof SF based on point SF */
-  if (debug) {
-    ierr = PetscPrintf(comm, "Input Section for Preallocation:\n");CHKERRQ(ierr);
-    ierr = PetscSectionView(section, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
-    ierr = PetscPrintf(comm, "Input Global Section for Preallocation:\n");CHKERRQ(ierr);
-    ierr = PetscSectionView(sectionGlobal, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
-    ierr = PetscPrintf(comm, "Input SF for Preallocation:\n");CHKERRQ(ierr);
-    ierr = PetscSFView(sf, NULL);CHKERRQ(ierr);
-  }
-  ierr = PetscSFCreateRemoteOffsets(sf, section, section, &remoteOffsets);CHKERRQ(ierr);
-  ierr = PetscSFCreateSectionSF(sf, section, remoteOffsets, section, &sfDof);CHKERRQ(ierr);
-  if (debug) {
-    ierr = PetscPrintf(comm, "Dof SF for Preallocation:\n");CHKERRQ(ierr);
-    ierr = PetscSFView(sfDof, NULL);CHKERRQ(ierr);
-  }
-  /* Create section for dof adjacency (dof ==> # adj dof) */
-  /*   FEM:   Two points p and q are adjacent if q \in closure(star(p)), preallocCenterDim = dim   */
-  /*   FVM:   Two points p and q are adjacent if q \in star(cone(p)),    preallocCenterDim = dim-1 */
-  /*   FVM++: Two points p and q are adjacent if q \in star(closure(p)), preallocCenterDim = 0     */
-  if (mesh->preallocCenterDim == dim) {
-    useClosure = PETSC_FALSE;
-  } else if (mesh->preallocCenterDim == 0) {
-    useClosure = PETSC_TRUE;
-  } else SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Do not support preallocation with center points of dimension %d", mesh->preallocCenterDim);
-
-  ierr = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr);
-  ierr = PetscSectionGetStorageSize(section, &numDof);CHKERRQ(ierr);
-  ierr = PetscSectionCreate(comm, &leafSectionAdj);CHKERRQ(ierr);
-  ierr = PetscSectionSetChart(leafSectionAdj, 0, numDof);CHKERRQ(ierr);
-  ierr = PetscSectionCreate(comm, &rootSectionAdj);CHKERRQ(ierr);
-  ierr = PetscSectionSetChart(rootSectionAdj, 0, numDof);CHKERRQ(ierr);
-  /*   Fill in the ghost dofs on the interface */
-  ierr = PetscSFGetGraph(sf, NULL, &nleaves, &leaves, &remotes);CHKERRQ(ierr);
-  ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr);
-  ierr = DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);CHKERRQ(ierr);
-
-  maxClosureSize = 2*PetscMax(PetscPowInt(mesh->maxConeSize,depth),PetscPowInt(mesh->maxSupportSize,depth)) + 2;
-  maxAdjSize     = PetscPowInt(mesh->maxConeSize,depth) * PetscPowInt(mesh->maxSupportSize,depth) + 1;
-
-  ierr = PetscMalloc2(maxClosureSize,PetscInt,&tmpClosure,maxAdjSize,PetscInt,&tmpAdj);CHKERRQ(ierr);
-
-  /*
-   ** The bootstrapping process involves six rounds with similar structure of visiting neighbors of each point.
-    1. Visit unowned points on interface, count adjacencies placing in leafSectionAdj
-       Reduce those counts to rootSectionAdj (now redundantly counting some interface points)
-    2. Visit owned points on interface, count adjacencies placing in rootSectionAdj
-       Create sfAdj connecting rootSectionAdj and leafSectionAdj
-    3. Visit unowned points on interface, write adjacencies to adj
-       Gather adj to rootAdj (note that there is redundancy in rootAdj when multiple procs find the same adjacencies)
-    4. Visit owned points on interface, write adjacencies to rootAdj
-       Remove redundancy in rootAdj
-   ** The last two traversals use transitive closure
-    5. Visit all owned points in the subdomain, count dofs for each point (sectionAdj)
-       Allocate memory addressed by sectionAdj (cols)
-    6. Visit all owned points in the subdomain, insert dof adjacencies into cols
-   ** Knowing all the column adjacencies, check ownership and sum into dnz and onz
-  */
-
-  for (l = 0; l < nleaves; ++l) {
-    PetscInt dof, off, d, q;
-    PetscInt p = leaves[l], numAdj = maxAdjSize;
-
-    ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
-    ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
-    ierr = DMPlexGetAdjacency_Private(dm, p, useClosure, tmpClosure, &numAdj, tmpAdj);CHKERRQ(ierr);
-    for (q = 0; q < numAdj; ++q) {
-      const PetscInt padj = tmpAdj[q];
-      PetscInt ndof, ncdof;
-
-      if ((padj < pStart) || (padj >= pEnd)) continue;
-      ierr = PetscSectionGetDof(section, padj, &ndof);CHKERRQ(ierr);
-      ierr = PetscSectionGetConstraintDof(section, padj, &ncdof);CHKERRQ(ierr);
-      for (d = off; d < off+dof; ++d) {
-        ierr = PetscSectionAddDof(leafSectionAdj, d, ndof-ncdof);CHKERRQ(ierr);
-      }
-    }
-  }
-  ierr = PetscSectionSetUp(leafSectionAdj);CHKERRQ(ierr);
-  if (debug) {
-    ierr = PetscPrintf(comm, "Adjacency Section for Preallocation on Leaves:\n");CHKERRQ(ierr);
-    ierr = PetscSectionView(leafSectionAdj, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
-  }
-  /* Get maximum remote adjacency sizes for owned dofs on interface (roots) */
-  if (size > 1) {
-    ierr = PetscSFReduceBegin(sfDof, MPIU_INT, leafSectionAdj->atlasDof, rootSectionAdj->atlasDof, MPI_SUM);CHKERRQ(ierr);
-    ierr = PetscSFReduceEnd(sfDof, MPIU_INT, leafSectionAdj->atlasDof, rootSectionAdj->atlasDof, MPI_SUM);CHKERRQ(ierr);
-  }
-  if (debug) {
-    ierr = PetscPrintf(comm, "Adjancency Section for Preallocation on Roots:\n");CHKERRQ(ierr);
-    ierr = PetscSectionView(rootSectionAdj, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
-  }
-  /* Add in local adjacency sizes for owned dofs on interface (roots) */
-  for (p = pStart; p < pEnd; ++p) {
-    PetscInt numAdj = maxAdjSize, adof, dof, off, d, q;
-
-    ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
-    ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
-    if (!dof) continue;
-    ierr = PetscSectionGetDof(rootSectionAdj, off, &adof);CHKERRQ(ierr);
-    if (adof <= 0) continue;
-    ierr = DMPlexGetAdjacency_Private(dm, p, useClosure, tmpClosure, &numAdj, tmpAdj);CHKERRQ(ierr);
-    for (q = 0; q < numAdj; ++q) {
-      const PetscInt padj = tmpAdj[q];
-      PetscInt ndof, ncdof;
-
-      if ((padj < pStart) || (padj >= pEnd)) continue;
-      ierr = PetscSectionGetDof(section, padj, &ndof);CHKERRQ(ierr);
-      ierr = PetscSectionGetConstraintDof(section, padj, &ncdof);CHKERRQ(ierr);
-      for (d = off; d < off+dof; ++d) {
-        ierr = PetscSectionAddDof(rootSectionAdj, d, ndof-ncdof);CHKERRQ(ierr);
-      }
-    }
-  }
-  ierr = PetscSectionSetUp(rootSectionAdj);CHKERRQ(ierr);
-  if (debug) {
-    ierr = PetscPrintf(comm, "Adjancency Section for Preallocation on Roots after local additions:\n");CHKERRQ(ierr);
-    ierr = PetscSectionView(rootSectionAdj, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
-  }
-  /* Create adj SF based on dof SF */
-  ierr = PetscSFCreateRemoteOffsets(sfDof, rootSectionAdj, leafSectionAdj, &remoteOffsets);CHKERRQ(ierr);
-  ierr = PetscSFCreateSectionSF(sfDof, rootSectionAdj, remoteOffsets, leafSectionAdj, &sfAdj);CHKERRQ(ierr);
-  if (debug) {
-    ierr = PetscPrintf(comm, "Adjacency SF for Preallocation:\n");CHKERRQ(ierr);
-    ierr = PetscSFView(sfAdj, NULL);CHKERRQ(ierr);
-  }
-  ierr = PetscSFDestroy(&sfDof);CHKERRQ(ierr);
-  /* Create leaf adjacency */
-  ierr = PetscSectionSetUp(leafSectionAdj);CHKERRQ(ierr);
-  ierr = PetscSectionGetStorageSize(leafSectionAdj, &adjSize);CHKERRQ(ierr);
-  ierr = PetscMalloc(adjSize * sizeof(PetscInt), &adj);CHKERRQ(ierr);
-  ierr = PetscMemzero(adj, adjSize * sizeof(PetscInt));CHKERRQ(ierr);
-  for (l = 0; l < nleaves; ++l) {
-    PetscInt dof, off, d, q;
-    PetscInt p = leaves[l], numAdj = maxAdjSize;
-
-    ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
-    ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
-    ierr = DMPlexGetAdjacency_Private(dm, p, useClosure, tmpClosure, &numAdj, tmpAdj);CHKERRQ(ierr);
-    for (d = off; d < off+dof; ++d) {
-      PetscInt aoff, i = 0;
-
-      ierr = PetscSectionGetOffset(leafSectionAdj, d, &aoff);CHKERRQ(ierr);
-      for (q = 0; q < numAdj; ++q) {
-        const PetscInt padj = tmpAdj[q];
-        PetscInt ndof, ncdof, ngoff, nd;
-
-        if ((padj < pStart) || (padj >= pEnd)) continue;
-        ierr = PetscSectionGetDof(section, padj, &ndof);CHKERRQ(ierr);
-        ierr = PetscSectionGetConstraintDof(section, padj, &ncdof);CHKERRQ(ierr);
-        ierr = PetscSectionGetOffset(sectionGlobal, padj, &ngoff);CHKERRQ(ierr);
-        for (nd = 0; nd < ndof-ncdof; ++nd) {
-          adj[aoff+i] = (ngoff < 0 ? -(ngoff+1) : ngoff) + nd;
-          ++i;
-        }
-      }
-    }
-  }
-  /* Debugging */
-  if (debug) {
-    IS tmp;
-    ierr = PetscPrintf(comm, "Leaf adjacency indices\n");CHKERRQ(ierr);
-    ierr = ISCreateGeneral(comm, adjSize, adj, PETSC_USE_POINTER, &tmp);CHKERRQ(ierr);
-    ierr = ISView(tmp, NULL);CHKERRQ(ierr);
-  }
-  /* Gather adjacenct indices to root */
-  ierr = PetscSectionGetStorageSize(rootSectionAdj, &adjSize);CHKERRQ(ierr);
-  ierr = PetscMalloc(adjSize * sizeof(PetscInt), &rootAdj);CHKERRQ(ierr);
-  for (r = 0; r < adjSize; ++r) rootAdj[r] = -1;
-  if (size > 1) {
-    ierr = PetscSFGatherBegin(sfAdj, MPIU_INT, adj, rootAdj);CHKERRQ(ierr);
-    ierr = PetscSFGatherEnd(sfAdj, MPIU_INT, adj, rootAdj);CHKERRQ(ierr);
-  }
-  ierr = PetscSFDestroy(&sfAdj);CHKERRQ(ierr);
-  ierr = PetscFree(adj);CHKERRQ(ierr);
-  /* Debugging */
-  if (debug) {
-    IS tmp;
-    ierr = PetscPrintf(comm, "Root adjacency indices after gather\n");CHKERRQ(ierr);
-    ierr = ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp);CHKERRQ(ierr);
-    ierr = ISView(tmp, NULL);CHKERRQ(ierr);
-  }
-  /* Add in local adjacency indices for owned dofs on interface (roots) */
-  for (p = pStart; p < pEnd; ++p) {
-    PetscInt numAdj = maxAdjSize, adof, dof, off, d, q;
-
-    ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
-    ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
-    if (!dof) continue;
-    ierr = PetscSectionGetDof(rootSectionAdj, off, &adof);CHKERRQ(ierr);
-    if (adof <= 0) continue;
-    ierr = DMPlexGetAdjacency_Private(dm, p, useClosure, tmpClosure, &numAdj, tmpAdj);CHKERRQ(ierr);
-    for (d = off; d < off+dof; ++d) {
-      PetscInt adof, aoff, i;
-
-      ierr = PetscSectionGetDof(rootSectionAdj, d, &adof);CHKERRQ(ierr);
-      ierr = PetscSectionGetOffset(rootSectionAdj, d, &aoff);CHKERRQ(ierr);
-      i    = adof-1;
-      for (q = 0; q < numAdj; ++q) {
-        const PetscInt padj = tmpAdj[q];
-        PetscInt ndof, ncdof, ngoff, nd;
-
-        if ((padj < pStart) || (padj >= pEnd)) continue;
-        ierr = PetscSectionGetDof(section, padj, &ndof);CHKERRQ(ierr);
-        ierr = PetscSectionGetConstraintDof(section, padj, &ncdof);CHKERRQ(ierr);
-        ierr = PetscSectionGetOffset(sectionGlobal, padj, &ngoff);CHKERRQ(ierr);
-        for (nd = 0; nd < ndof-ncdof; ++nd) {
-          rootAdj[aoff+i] = ngoff < 0 ? -(ngoff+1)+nd : ngoff+nd;
-          --i;
-        }
-      }
-    }
-  }
-  /* Debugging */
-  if (debug) {
-    IS tmp;
-    ierr = PetscPrintf(comm, "Root adjacency indices\n");CHKERRQ(ierr);
-    ierr = ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp);CHKERRQ(ierr);
-    ierr = ISView(tmp, NULL);CHKERRQ(ierr);
-  }
-  /* Compress indices */
-  ierr = PetscSectionSetUp(rootSectionAdj);CHKERRQ(ierr);
-  for (p = pStart; p < pEnd; ++p) {
-    PetscInt dof, cdof, off, d;
-    PetscInt adof, aoff;
-
-    ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
-    ierr = PetscSectionGetConstraintDof(section, p, &cdof);CHKERRQ(ierr);
-    ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
-    if (!dof) continue;
-    ierr = PetscSectionGetDof(rootSectionAdj, off, &adof);CHKERRQ(ierr);
-    if (adof <= 0) continue;
-    for (d = off; d < off+dof-cdof; ++d) {
-      ierr = PetscSectionGetDof(rootSectionAdj, d, &adof);CHKERRQ(ierr);
-      ierr = PetscSectionGetOffset(rootSectionAdj, d, &aoff);CHKERRQ(ierr);
-      ierr = PetscSortRemoveDupsInt(&adof, &rootAdj[aoff]);CHKERRQ(ierr);
-      ierr = PetscSectionSetDof(rootSectionAdj, d, adof);CHKERRQ(ierr);
-    }
-  }
-  /* Debugging */
-  if (debug) {
-    IS tmp;
-    ierr = PetscPrintf(comm, "Adjancency Section for Preallocation on Roots after compression:\n");CHKERRQ(ierr);
-    ierr = PetscSectionView(rootSectionAdj, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
-    ierr = PetscPrintf(comm, "Root adjacency indices after compression\n");CHKERRQ(ierr);
-    ierr = ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp);CHKERRQ(ierr);
-    ierr = ISView(tmp, NULL);CHKERRQ(ierr);
-  }
-  /* Build adjacency section: Maps global indices to sets of adjacent global indices */
-  ierr = PetscSectionGetOffsetRange(sectionGlobal, &globalOffStart, &globalOffEnd);CHKERRQ(ierr);
-  ierr = PetscSectionCreate(comm, &sectionAdj);CHKERRQ(ierr);
-  ierr = PetscSectionSetChart(sectionAdj, globalOffStart, globalOffEnd);CHKERRQ(ierr);
-  for (p = pStart; p < pEnd; ++p) {
-    PetscInt  numAdj = maxAdjSize, dof, cdof, off, goff, d, q;
-    PetscBool found  = PETSC_TRUE;
-
-    ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
-    ierr = PetscSectionGetConstraintDof(section, p, &cdof);CHKERRQ(ierr);
-    ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
-    ierr = PetscSectionGetOffset(sectionGlobal, p, &goff);CHKERRQ(ierr);
-    for (d = 0; d < dof-cdof; ++d) {
-      PetscInt ldof, rdof;
-
-      ierr = PetscSectionGetDof(leafSectionAdj, off+d, &ldof);CHKERRQ(ierr);
-      ierr = PetscSectionGetDof(rootSectionAdj, off+d, &rdof);CHKERRQ(ierr);
-      if (ldof > 0) {
-        /* We do not own this point */
-      } else if (rdof > 0) {
-        ierr = PetscSectionSetDof(sectionAdj, goff+d, rdof);CHKERRQ(ierr);
-      } else {
-        found = PETSC_FALSE;
-      }
-    }
-    if (found) continue;
-    ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
-    ierr = PetscSectionGetOffset(sectionGlobal, p, &goff);CHKERRQ(ierr);
-    ierr = DMPlexGetAdjacency_Private(dm, p, useClosure, tmpClosure, &numAdj, tmpAdj);CHKERRQ(ierr);
-    for (q = 0; q < numAdj; ++q) {
-      const PetscInt padj = tmpAdj[q];
-      PetscInt ndof, ncdof, noff;
-
-      if ((padj < pStart) || (padj >= pEnd)) continue;
-      ierr = PetscSectionGetDof(section, padj, &ndof);CHKERRQ(ierr);
-      ierr = PetscSectionGetConstraintDof(section, padj, &ncdof);CHKERRQ(ierr);
-      ierr = PetscSectionGetOffset(section, padj, &noff);CHKERRQ(ierr);
-      for (d = goff; d < goff+dof-cdof; ++d) {
-        ierr = PetscSectionAddDof(sectionAdj, d, ndof-ncdof);CHKERRQ(ierr);
-      }
-    }
-  }
-  ierr = PetscSectionSetUp(sectionAdj);CHKERRQ(ierr);
-  if (debug) {
-    ierr = PetscPrintf(comm, "Adjacency Section for Preallocation:\n");CHKERRQ(ierr);
-    ierr = PetscSectionView(sectionAdj, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
-  }
-  /* Get adjacent indices */
-  ierr = PetscSectionGetStorageSize(sectionAdj, &numCols);CHKERRQ(ierr);
-  ierr = PetscMalloc(numCols * sizeof(PetscInt), &cols);CHKERRQ(ierr);
-  for (p = pStart; p < pEnd; ++p) {
-    PetscInt  numAdj = maxAdjSize, dof, cdof, off, goff, d, q;
-    PetscBool found  = PETSC_TRUE;
-
-    ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
-    ierr = PetscSectionGetConstraintDof(section, p, &cdof);CHKERRQ(ierr);
-    ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
-    ierr = PetscSectionGetOffset(sectionGlobal, p, &goff);CHKERRQ(ierr);
-    for (d = 0; d < dof-cdof; ++d) {
-      PetscInt ldof, rdof;
-
-      ierr = PetscSectionGetDof(leafSectionAdj, off+d, &ldof);CHKERRQ(ierr);
-      ierr = PetscSectionGetDof(rootSectionAdj, off+d, &rdof);CHKERRQ(ierr);
-      if (ldof > 0) {
-        /* We do not own this point */
-      } else if (rdof > 0) {
-        PetscInt aoff, roff;
-
-        ierr = PetscSectionGetOffset(sectionAdj, goff+d, &aoff);CHKERRQ(ierr);
-        ierr = PetscSectionGetOffset(rootSectionAdj, off+d, &roff);CHKERRQ(ierr);
-        ierr = PetscMemcpy(&cols[aoff], &rootAdj[roff], rdof * sizeof(PetscInt));CHKERRQ(ierr);
-      } else {
-        found = PETSC_FALSE;
-      }
-    }
-    if (found) continue;
-    ierr = DMPlexGetAdjacency_Private(dm, p, useClosure, tmpClosure, &numAdj, tmpAdj);CHKERRQ(ierr);
-    for (d = goff; d < goff+dof-cdof; ++d) {
-      PetscInt adof, aoff, i = 0;
-
-      ierr = PetscSectionGetDof(sectionAdj, d, &adof);CHKERRQ(ierr);
-      ierr = PetscSectionGetOffset(sectionAdj, d, &aoff);CHKERRQ(ierr);
-      for (q = 0; q < numAdj; ++q) {
-        const PetscInt  padj = tmpAdj[q];
-        PetscInt        ndof, ncdof, ngoff, nd;
-        const PetscInt *ncind;
-
-        /* Adjacent points may not be in the section chart */
-        if ((padj < pStart) || (padj >= pEnd)) continue;
-        ierr = PetscSectionGetDof(section, padj, &ndof);CHKERRQ(ierr);
-        ierr = PetscSectionGetConstraintDof(section, padj, &ncdof);CHKERRQ(ierr);
-        ierr = PetscSectionGetConstraintIndices(section, padj, &ncind);CHKERRQ(ierr);
-        ierr = PetscSectionGetOffset(sectionGlobal, padj, &ngoff);CHKERRQ(ierr);
-        for (nd = 0; nd < ndof-ncdof; ++nd, ++i) {
-          cols[aoff+i] = ngoff < 0 ? -(ngoff+1)+nd : ngoff+nd;
-        }
-      }
-      if (i != adof) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of entries %D != %D for dof %D (point %D)", i, adof, d, p);
-    }
-  }
-  ierr = PetscSectionDestroy(&leafSectionAdj);CHKERRQ(ierr);
-  ierr = PetscSectionDestroy(&rootSectionAdj);CHKERRQ(ierr);
-  ierr = PetscFree(rootAdj);CHKERRQ(ierr);
-  ierr = PetscFree2(tmpClosure, tmpAdj);CHKERRQ(ierr);
-  /* Debugging */
-  if (debug) {
-    IS tmp;
-    ierr = PetscPrintf(comm, "Column indices\n");CHKERRQ(ierr);
-    ierr = ISCreateGeneral(comm, numCols, cols, PETSC_USE_POINTER, &tmp);CHKERRQ(ierr);
-    ierr = ISView(tmp, NULL);CHKERRQ(ierr);
-  }
-  /* Create allocation vectors from adjacency graph */
-  ierr = MatGetLocalSize(A, &locRows, NULL);CHKERRQ(ierr);
-  ierr = PetscLayoutCreate(PetscObjectComm((PetscObject)A), &rLayout);CHKERRQ(ierr);
-  ierr = PetscLayoutSetLocalSize(rLayout, locRows);CHKERRQ(ierr);
-  ierr = PetscLayoutSetBlockSize(rLayout, 1);CHKERRQ(ierr);
-  ierr = PetscLayoutSetUp(rLayout);CHKERRQ(ierr);
-  ierr = PetscLayoutGetRange(rLayout, &rStart, &rEnd);CHKERRQ(ierr);
-  ierr = PetscLayoutDestroy(&rLayout);CHKERRQ(ierr);
-  /* Only loop over blocks of rows */
-  if (rStart%bs || rEnd%bs) SETERRQ3(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Invalid layout [%d, %d) for matrix, must be divisible by block size %d", rStart, rEnd, bs);
-  for (r = rStart/bs; r < rEnd/bs; ++r) {
-    const PetscInt row = r*bs;
-    PetscInt       numCols, cStart, c;
-
-    ierr = PetscSectionGetDof(sectionAdj, row, &numCols);CHKERRQ(ierr);
-    ierr = PetscSectionGetOffset(sectionAdj, row, &cStart);CHKERRQ(ierr);
-    for (c = cStart; c < cStart+numCols; ++c) {
-      if ((cols[c] >= rStart*bs) && (cols[c] < rEnd*bs)) {
-        ++dnz[r-rStart];
-        if (cols[c] >= row) ++dnzu[r-rStart];
-      } else {
-        ++onz[r-rStart];
-        if (cols[c] >= row) ++onzu[r-rStart];
-      }
-    }
-  }
-  if (bs > 1) {
-    for (r = 0; r < locRows/bs; ++r) {
-      dnz[r]  /= bs;
-      onz[r]  /= bs;
-      dnzu[r] /= bs;
-      onzu[r] /= bs;
-    }
-  }
-  /* Set matrix pattern */
-  ierr = MatXAIJSetPreallocation(A, bs, dnz, onz, dnzu, onzu);CHKERRQ(ierr);
-  ierr = MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);
-  /* Fill matrix with zeros */
-  if (fillMatrix) {
-    PetscScalar *values;
-    PetscInt     maxRowLen = 0;
-
-    for (r = rStart; r < rEnd; ++r) {
-      PetscInt len;
-
-      ierr      = PetscSectionGetDof(sectionAdj, r, &len);CHKERRQ(ierr);
-      maxRowLen = PetscMax(maxRowLen, len);
-    }
-    ierr = PetscMalloc(maxRowLen * sizeof(PetscScalar), &values);CHKERRQ(ierr);
-    ierr = PetscMemzero(values, maxRowLen * sizeof(PetscScalar));CHKERRQ(ierr);
-    for (r = rStart; r < rEnd; ++r) {
-      PetscInt numCols, cStart;
-
-      ierr = PetscSectionGetDof(sectionAdj, r, &numCols);CHKERRQ(ierr);
-      ierr = PetscSectionGetOffset(sectionAdj, r, &cStart);CHKERRQ(ierr);
-      ierr = MatSetValues(A, 1, &r, numCols, &cols[cStart], values, INSERT_VALUES);CHKERRQ(ierr);
-    }
-    ierr = PetscFree(values);CHKERRQ(ierr);
-    ierr = MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
-    ierr = MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
-  }
-  ierr = PetscSectionDestroy(&sectionAdj);CHKERRQ(ierr);
-  ierr = PetscFree(cols);CHKERRQ(ierr);
-  PetscFunctionReturn(0);
-}
-
-#if 0
-#undef __FUNCT__
-#define __FUNCT__ "DMPlexPreallocateOperator_2"
-PetscErrorCode DMPlexPreallocateOperator_2(DM dm, PetscInt bs, PetscSection section, PetscSection sectionGlobal, PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[], Mat A, PetscBool fillMatrix)
-{
-  PetscInt       *tmpClosure,*tmpAdj,*visits;
-  PetscInt        c,cStart,cEnd,pStart,pEnd;
-  PetscErrorCode  ierr;
-
-  PetscFunctionBegin;
-  ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
-  ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr);
-  ierr = DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);CHKERRQ(ierr);
-
-  maxClosureSize = 2*PetscMax(PetscPowInt(mesh->maxConeSize,depth),PetscPowInt(mesh->maxSupportSize,depth));
-
-  ierr    = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr);
-  npoints = pEnd - pStart;
-
-  ierr = PetscMalloc3(maxClosureSize,PetscInt,&tmpClosure,npoints,PetscInt,&lvisits,npoints,PetscInt,&visits);CHKERRQ(ierr);
-  ierr = PetscMemzero(lvisits,(pEnd-pStart)*sizeof(PetscInt));CHKERRQ(ierr);
-  ierr = PetscMemzero(visits,(pEnd-pStart)*sizeof(PetscInt));CHKERRQ(ierr);
-  ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
-  for (c=cStart; c<cEnd; c++) {
-    PetscInt *support = tmpClosure;
-    ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_FALSE, &supportSize, (PetscInt**)&support);CHKERRQ(ierr);
-    for (p=0; p<supportSize; p++) lvisits[support[p]]++;
-  }
-  ierr = PetscSFReduceBegin(sf,MPIU_INT,lvisits,visits,MPI_SUM);CHKERRQ(ierr);
-  ierr = PetscSFReduceEnd  (sf,MPIU_INT,lvisits,visits,MPI_SUM);CHKERRQ(ierr);
-  ierr = PetscSFBcastBegin(sf,MPIU_INT,visits,lvisits);CHKERRQ(ierr);
-  ierr = PetscSFBcastEnd  (sf,MPIU_INT,visits,lvisits);CHKERRQ(ierr);
-
-  ierr = PetscSFGetRanks();CHKERRQ(ierr);
-
-
-  ierr = PetscMalloc2(maxClosureSize*maxClosureSize,PetscInt,&cellmat,npoints,PetscInt,&owner);CHKERRQ(ierr);
-  for (c=cStart; c<cEnd; c++) {
-    ierr = PetscMemzero(cellmat,maxClosureSize*maxClosureSize*sizeof(PetscInt));CHKERRQ(ierr);
-    /*
-     Depth-first walk of transitive closure.
-     At each leaf frame f of transitive closure that we see, add 1/visits[f] to each pair (p,q) not marked as done in cellmat.
-     This contribution is added to dnz if owning ranks of p and q match, to onz otherwise.
-     */
-  }
-
-  ierr = PetscSFReduceBegin(sf,MPIU_INT,ldnz,dnz,MPI_SUM);CHKERRQ(ierr);
-  ierr = PetscSFReduceEnd  (sf,MPIU_INT,lonz,onz,MPI_SUM);CHKERRQ(ierr);
-  PetscFunctionReturn(0);
-}
-#endif
-
-#undef __FUNCT__
 #define __FUNCT__ "DMCreateMatrix_Plex"
 PetscErrorCode DMCreateMatrix_Plex(DM dm, MatType mtype, Mat *J)
 {
 }
 
 #undef __FUNCT__
+#define __FUNCT__ "DMPlexGetAdjacencySingleLevel_Internal"
+static PetscErrorCode DMPlexGetAdjacencySingleLevel_Internal(DM dm, PetscInt p, PetscBool useClosure, const PetscInt *tmpClosure, PetscInt *adjSize, PetscInt adj[])
+{
+  const PetscInt *support = NULL;
+  PetscInt        numAdj   = 0, maxAdjSize = *adjSize, supportSize, s;
+  PetscErrorCode  ierr;
+
+  PetscFunctionBegin;
+  if (useClosure) {
+    ierr = DMPlexGetConeSize(dm, p, &supportSize);CHKERRQ(ierr);
+    ierr = DMPlexGetCone(dm, p, &support);CHKERRQ(ierr);
+    for (s = 0; s < supportSize; ++s) {
+      const PetscInt *cone = NULL;
+      PetscInt        coneSize, c, q;
+
+      ierr = DMPlexGetSupportSize(dm, support[s], &coneSize);CHKERRQ(ierr);
+      ierr = DMPlexGetSupport(dm, support[s], &cone);CHKERRQ(ierr);
+      for (c = 0; c < coneSize; ++c) {
+        for (q = 0; q < numAdj || (adj[numAdj++] = cone[c],0); ++q) {
+          if (cone[c] == adj[q]) break;
+        }
+        if (numAdj > maxAdjSize) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid mesh exceeded adjacency allocation (%D)", maxAdjSize);
+      }
+    }
+  } else {
+    ierr = DMPlexGetSupportSize(dm, p, &supportSize);CHKERRQ(ierr);
+    ierr = DMPlexGetSupport(dm, p, &support);CHKERRQ(ierr);
+    for (s = 0; s < supportSize; ++s) {
+      const PetscInt *cone = NULL;
+      PetscInt        coneSize, c, q;
+
+      ierr = DMPlexGetConeSize(dm, support[s], &coneSize);CHKERRQ(ierr);
+      ierr = DMPlexGetCone(dm, support[s], &cone);CHKERRQ(ierr);
+      for (c = 0; c < coneSize; ++c) {
+        for (q = 0; q < numAdj || (adj[numAdj++] = cone[c],0); ++q) {
+          if (cone[c] == adj[q]) break;
+        }
+        if (numAdj > maxAdjSize) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid mesh exceeded adjacency allocation (%D)", maxAdjSize);
+      }
+    }
+  }
+  *adjSize = numAdj;
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
 #define __FUNCT__ "DMPlexCreateNeighborCSR"
 PetscErrorCode DMPlexCreateNeighborCSR(DM dm, PetscInt cellHeight, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency)
 {
   for (cell = cStart; cell < cEnd; ++cell) {
     PetscInt numNeighbors = maxNeighbors, n;
 
-    ierr = DMPlexGetAdjacencySingleLevel_Private(dm, cell, PETSC_TRUE, tmpClosure, &numNeighbors, neighborCells);CHKERRQ(ierr);
+    ierr = DMPlexGetAdjacencySingleLevel_Internal(dm, cell, PETSC_TRUE, tmpClosure, &numNeighbors, neighborCells);CHKERRQ(ierr);
     /* Get meet with each cell, and check with recognizer (could optimize to check each pair only once) */
     for (n = 0; n < numNeighbors; ++n) {
       PetscInt        cellPair[2];
       PetscInt numNeighbors = maxNeighbors, n;
       PetscInt cellOffset   = 0;
 
-      ierr = DMPlexGetAdjacencySingleLevel_Private(dm, cell, PETSC_TRUE, tmpClosure, &numNeighbors, neighborCells);CHKERRQ(ierr);
+      ierr = DMPlexGetAdjacencySingleLevel_Internal(dm, cell, PETSC_TRUE, tmpClosure, &numNeighbors, neighborCells);CHKERRQ(ierr);
       /* Get meet with each cell, and check with recognizer (could optimize to check each pair only once) */
       for (n = 0; n < numNeighbors; ++n) {
         PetscInt        cellPair[2];
   PetscFunctionReturn(0);
 }
 
-PETSC_STATIC_INLINE PetscInt epsilon(PetscInt i, PetscInt j, PetscInt k)
-{
-  switch (i) {
-  case 0:
-    switch (j) {
-    case 0: return 0;
-    case 1:
-      switch (k) {
-      case 0: return 0;
-      case 1: return 0;
-      case 2: return 1;
-      }
-    case 2:
-      switch (k) {
-      case 0: return 0;
-      case 1: return -1;
-      case 2: return 0;
-      }
-    }
-  case 1:
-    switch (j) {
-    case 0:
-      switch (k) {
-      case 0: return 0;
-      case 1: return 0;
-      case 2: return -1;
-      }
-    case 1: return 0;
-    case 2:
-      switch (k) {
-      case 0: return 1;
-      case 1: return 0;
-      case 2: return 0;
-      }
-    }
-  case 2:
-    switch (j) {
-    case 0:
-      switch (k) {
-      case 0: return 0;
-      case 1: return 1;
-      case 2: return 0;
-      }
-    case 1:
-      switch (k) {
-      case 0: return -1;
-      case 1: return 0;
-      case 2: return 0;
-      }
-    case 2: return 0;
-    }
-  }
-  return 0;
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "DMPlexCreateRigidBody"
-/*@C
-  DMPlexCreateRigidBody - create rigid body modes from coordinates
-
-  Collective on DM
-
-  Input Arguments:
-+ dm - the DM
-. section - the local section associated with the rigid field, or NULL for the default section
-- globalSection - the global section associated with the rigid field, or NULL for the default section
-
-  Output Argument:
-. sp - the null space
-
-  Note: This is necessary to take account of Dirichlet conditions on the displacements
-
-  Level: advanced
-
-.seealso: MatNullSpaceCreate()
-@*/
-PetscErrorCode DMPlexCreateRigidBody(DM dm, PetscSection section, PetscSection globalSection, MatNullSpace *sp)
-{
-  MPI_Comm       comm;
-  Vec            coordinates, localMode, mode[6];
-  PetscSection   coordSection;
-  PetscScalar   *coords;
-  PetscInt       dim, vStart, vEnd, v, n, m, d, i, j;
-  PetscErrorCode ierr;
-
-  PetscFunctionBegin;
-  ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
-  ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
-  if (dim == 1) {
-    ierr = MatNullSpaceCreate(comm, PETSC_TRUE, 0, NULL, sp);CHKERRQ(ierr);
-    PetscFunctionReturn(0);
-  }
-  if (!section)       {ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);}
-  if (!globalSection) {ierr = DMGetDefaultGlobalSection(dm, &globalSection);CHKERRQ(ierr);}
-  ierr = PetscSectionGetConstrainedStorageSize(globalSection, &n);CHKERRQ(ierr);
-  ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
-  ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
-  ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
-  m    = (dim*(dim+1))/2;
-  ierr = VecCreate(comm, &mode[0]);CHKERRQ(ierr);
-  ierr = VecSetSizes(mode[0], n, PETSC_DETERMINE);CHKERRQ(ierr);
-  ierr = VecSetUp(mode[0]);CHKERRQ(ierr);
-  for (i = 1; i < m; ++i) {ierr = VecDuplicate(mode[0], &mode[i]);CHKERRQ(ierr);}
-  /* Assume P1 */
-  ierr = DMGetLocalVector(dm, &localMode);CHKERRQ(ierr);
-  for (d = 0; d < dim; ++d) {
-    PetscScalar values[3] = {0.0, 0.0, 0.0};
-
-    values[d] = 1.0;
-    ierr      = VecSet(localMode, 0.0);CHKERRQ(ierr);
-    for (v = vStart; v < vEnd; ++v) {
-      ierr = DMPlexVecSetClosure(dm, section, localMode, v, values, INSERT_VALUES);CHKERRQ(ierr);
-    }
-    ierr = DMLocalToGlobalBegin(dm, localMode, INSERT_VALUES, mode[d]);CHKERRQ(ierr);
-    ierr = DMLocalToGlobalEnd(dm, localMode, INSERT_VALUES, mode[d]);CHKERRQ(ierr);
-  }
-  ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr);
-  for (d = dim; d < dim*(dim+1)/2; ++d) {
-    PetscInt i, j, k = dim > 2 ? d - dim : d;
-
-    ierr = VecSet(localMode, 0.0);CHKERRQ(ierr);
-    for (v = vStart; v < vEnd; ++v) {
-      PetscScalar values[3] = {0.0, 0.0, 0.0};
-      PetscInt    off;
-
-      ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr);
-      for (i = 0; i < dim; ++i) {
-        for (j = 0; j < dim; ++j) {
-          values[j] += epsilon(i, j, k)*PetscRealPart(coords[off+i]);
-        }
-      }
-      ierr = DMPlexVecSetClosure(dm, section, localMode, v, values, INSERT_VALUES);CHKERRQ(ierr);
-    }
-    ierr = DMLocalToGlobalBegin(dm, localMode, INSERT_VALUES, mode[d]);CHKERRQ(ierr);
-    ierr = DMLocalToGlobalEnd(dm, localMode, INSERT_VALUES, mode[d]);CHKERRQ(ierr);
-  }
-  ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr);
-  ierr = DMRestoreLocalVector(dm, &localMode);CHKERRQ(ierr);
-  for (i = 0; i < dim; ++i) {ierr = VecNormalize(mode[i], NULL);CHKERRQ(ierr);}
-  /* Orthonormalize system */
-  for (i = dim; i < m; ++i) {
-    PetscScalar dots[6];
-
-    ierr = VecMDot(mode[i], i, mode, dots);CHKERRQ(ierr);
-    for (j = 0; j < i; ++j) dots[j] *= -1.0;
-    ierr = VecMAXPY(mode[i], i, dots, mode);CHKERRQ(ierr);
-    ierr = VecNormalize(mode[i], NULL);CHKERRQ(ierr);
-  }
-  ierr = MatNullSpaceCreate(comm, PETSC_FALSE, m, mode, sp);CHKERRQ(ierr);
-  for (i = 0; i< m; ++i) {ierr = VecDestroy(&mode[i]);CHKERRQ(ierr);}
-  PetscFunctionReturn(0);
-}
-
 #undef __FUNCT__
 #define __FUNCT__ "DMPlexGetHybridBounds"
 PetscErrorCode DMPlexGetHybridBounds(DM dm, PetscInt *cMax, PetscInt *fMax, PetscInt *eMax, PetscInt *vMax)
   PetscFunctionReturn(0);
 }
 
-#undef __FUNCT__
-#define __FUNCT__ "DMPlexGetScale"
-PetscErrorCode DMPlexGetScale(DM dm, PetscUnit unit, PetscReal *scale)
-{
-  DM_Plex *mesh = (DM_Plex*) dm->data;
-
-  PetscFunctionBegin;
-  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
-  PetscValidPointer(scale, 3);
-  *scale = mesh->scale[unit];
-  PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "DMPlexSetScale"
-PetscErrorCode DMPlexSetScale(DM dm, PetscUnit unit, PetscReal scale)
-{
-  DM_Plex *mesh = (DM_Plex*) dm->data;
-
-  PetscFunctionBegin;
-  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
-  mesh->scale[unit] = scale;
-  PetscFunctionReturn(0);
-}
-
-
-/*******************************************************************************
-This should be in a separate Discretization object, but I am not sure how to lay
-it out yet, so I am stuffing things here while I experiment.
-*******************************************************************************/
-#undef __FUNCT__
-#define __FUNCT__ "DMPlexSetFEMIntegration"
-PetscErrorCode DMPlexSetFEMIntegration(DM dm,
-                                          PetscErrorCode (*integrateResidualFEM)(PetscInt, PetscInt, PetscInt, PetscQuadrature[], const PetscScalar[],
-                                                                                 const PetscReal[], const PetscReal[], const PetscReal[], const PetscReal[],
-                                                                                 void (*)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]),
-                                                                                 void (*)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]), PetscScalar[]),
-                                          PetscErrorCode (*integrateJacobianActionFEM)(PetscInt, PetscInt, PetscInt, PetscQuadrature[], const PetscScalar[], const PetscScalar[],
-                                                                                       const PetscReal[], const PetscReal[], const PetscReal[], const PetscReal[],
-                                                                                       void (**)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]),
-                                                                                       void (**)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]),
-                                                                                       void (**)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]),
-                                                                                       void (**)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]), PetscScalar[]),
-                                          PetscErrorCode (*integrateJacobianFEM)(PetscInt, PetscInt, PetscInt, PetscInt, PetscQuadrature[], const PetscScalar[],
-                                                                                 const PetscReal[], const PetscReal[], const PetscReal[], const PetscReal[],
-                                                                                 void (*)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]),
-                                                                                 void (*)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]),
-                                                                                 void (*)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]),
-                                                                                 void (*)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]), PetscScalar[]))
-{
-  DM_Plex *mesh = (DM_Plex*) dm->data;
-
-  PetscFunctionBegin;
-  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
-  mesh->integrateResidualFEM       = integrateResidualFEM;
-  mesh->integrateJacobianActionFEM = integrateJacobianActionFEM;
-  mesh->integrateJacobianFEM       = integrateJacobianFEM;
-  PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "DMPlexProjectFunctionLocal"
-PetscErrorCode DMPlexProjectFunctionLocal(DM dm, PetscInt numComp, PetscScalar (**funcs)(const PetscReal []), InsertMode mode, Vec localX)
-{
-  Vec            coordinates;
-  PetscSection   section, cSection;
-  PetscInt       dim, vStart, vEnd, v, c, d;
-  PetscScalar   *values, *cArray;
-  PetscReal     *coords;
-  PetscErrorCode ierr;
-
-  PetscFunctionBegin;
-  ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
-  ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
-  ierr = DMPlexGetCoordinateSection(dm, &cSection);CHKERRQ(ierr);
-  ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
-  ierr = PetscMalloc(numComp * sizeof(PetscScalar), &values);CHKERRQ(ierr);
-  ierr = VecGetArray(coordinates, &cArray);CHKERRQ(ierr);
-  ierr = PetscSectionGetDof(cSection, vStart, &dim);CHKERRQ(ierr);
-  ierr = PetscMalloc(dim * sizeof(PetscReal),&coords);CHKERRQ(ierr);
-  for (v = vStart; v < vEnd; ++v) {
-    PetscInt dof, off;
-
-    ierr = PetscSectionGetDof(cSection, v, &dof);CHKERRQ(ierr);
-    ierr = PetscSectionGetOffset(cSection, v, &off);CHKERRQ(ierr);
-    if (dof > dim) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Cannot have more coordinates %d then dimensions %d", dof, dim);
-    for (d = 0; d < dof; ++d) coords[d] = PetscRealPart(cArray[off+d]);
-    for (c = 0; c < numComp; ++c) values[c] = (*funcs[c])(coords);
-    ierr = VecSetValuesSection(localX, section, v, values, mode);CHKERRQ(ierr);
-  }
-  ierr = VecRestoreArray(coordinates, &cArray);CHKERRQ(ierr);
-  /* Temporary, must be replaced by a projection on the finite element basis */
-  {
-    PetscInt eStart = 0, eEnd = 0, e, depth;
-
-    ierr = DMPlexGetLabelSize(dm, "depth", &depth);CHKERRQ(ierr);
-    --depth;
-    if (depth > 1) {ierr = DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);CHKERRQ(ierr);}
-    for (e = eStart; e < eEnd; ++e) {
-      const PetscInt *cone = NULL;
-      PetscInt        coneSize, d;
-      PetscScalar    *coordsA, *coordsB;
-
-      ierr = DMPlexGetConeSize(dm, e, &coneSize);CHKERRQ(ierr);
-      ierr = DMPlexGetCone(dm, e, &cone);CHKERRQ(ierr);
-      if (coneSize != 2) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_SIZ, "Cone size %d for point %d should be 2", coneSize, e);
-      ierr = VecGetValuesSection(coordinates, cSection, cone[0], &coordsA);CHKERRQ(ierr);
-      ierr = VecGetValuesSection(coordinates, cSection, cone[1], &coordsB);CHKERRQ(ierr);
-      for (d = 0; d < dim; ++d) {
-        coords[d] = 0.5*(PetscRealPart(coordsA[d]) + PetscRealPart(coordsB[d]));
-      }
-      for (c = 0; c < numComp; ++c) values[c] = (*funcs[c])(coords);
-      ierr = VecSetValuesSection(localX, section, e, values, mode);CHKERRQ(ierr);
-    }
-  }
-
-  ierr = PetscFree(coords);CHKERRQ(ierr);
-  ierr = PetscFree(values);CHKERRQ(ierr);
-#if 0
-  const PetscInt localDof = this->_mesh->sizeWithBC(s, *cells->begin());
-  PetscReal      detJ;
-
-  ierr = PetscMalloc(localDof * sizeof(PetscScalar), &values);CHKERRQ(ierr);
-  ierr = PetscMalloc2(dim,PetscReal,&v0,dim*dim,PetscReal,&J);CHKERRQ(ierr);
-  ALE::ISieveVisitor::PointRetriever<PETSC_MESH_TYPE::sieve_type> pV(PetscPowInt(this->_mesh->getSieve()->getMaxConeSize(),dim+1), true);
-
-  for (PetscInt c = cStart; c < cEnd; ++c) {
-    ALE::ISieveTraversal<PETSC_MESH_TYPE::sieve_type>::orientedClosure(*this->_mesh->getSieve(), c, pV);
-    const PETSC_MESH_TYPE::point_type *oPoints = pV.getPoints();
-    const int                          oSize   = pV.getSize();
-    int                                v       = 0;
-
-    ierr = DMPlexComputeCellGeometry(dm, c, v0, J, NULL, &detJ);CHKERRQ(ierr);
-    for (PetscInt cl = 0; cl < oSize; ++cl) {
-      const PetscInt fDim;
-
-      ierr = PetscSectionGetDof(oPoints[cl], &fDim);CHKERRQ(ierr);
-      if (pointDim) {
-        for (PetscInt d = 0; d < fDim; ++d, ++v) {
-          values[v] = (*this->_options.integrate)(v0, J, v, initFunc);
-        }
-      }
-    }
-    ierr = DMPlexVecSetClosure(dm, NULL, localX, c, values);CHKERRQ(ierr);
-    pV.clear();
-  }
-  ierr = PetscFree2(v0,J);CHKERRQ(ierr);
-  ierr = PetscFree(values);CHKERRQ(ierr);
-#endif
-  PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "DMPlexProjectFunction"
-/*@C
-  DMPlexProjectFunction - This projects the given function into the function space provided.
-
-  Input Parameters:
-+ dm      - The DM
-. numComp - The number of components (functions)
-. funcs   - The coordinate functions to evaluate
-- mode    - The insertion mode for values
-
-  Output Parameter:
-. X - vector
-
-  Level: developer
-
-  Note:
-  This currently just calls the function with the coordinates of each vertex and edge midpoint, and stores the result in a vector.
-  We will eventually fix it.
-
-,seealso: DMPlexComputeL2Diff()
-*/
-PetscErrorCode DMPlexProjectFunction(DM dm, PetscInt numComp, PetscScalar (**funcs)(const PetscReal []), InsertMode mode, Vec X)
-{
-  Vec            localX;
-  PetscErrorCode ierr;
-
-  PetscFunctionBegin;
-  ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr);
-  ierr = DMPlexProjectFunctionLocal(dm, numComp, funcs, mode, localX);CHKERRQ(ierr);
-  ierr = DMLocalToGlobalBegin(dm, localX, mode, X);CHKERRQ(ierr);
-  ierr = DMLocalToGlobalEnd(dm, localX, mode, X);CHKERRQ(ierr);
-  ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr);
-  PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "DMPlexComputeL2Diff"
-/*@C
-  DMPlexComputeL2Diff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h.
-
-  Input Parameters:
-+ dm    - The DM
-. quad  - The PetscQuadrature object for each field
-. funcs - The functions to evaluate for each field component
-- X     - The coefficient vector u_h
-
-  Output Parameter:
-. diff - The diff ||u - u_h||_2
-
-  Level: developer
-
-.seealso: DMPlexProjectFunction()
-*/
-PetscErrorCode DMPlexComputeL2Diff(DM dm, PetscQuadrature quad[], PetscScalar (**funcs)(const PetscReal []), Vec X, PetscReal *diff)
-{
-  const PetscInt debug = 0;
-  PetscSection   section;
-  Vec            localX;
-  PetscReal     *coords, *v0, *J, *invJ, detJ;
-  PetscReal      localDiff = 0.0;
-  PetscInt       dim, numFields, numComponents = 0, cStart, cEnd, c, field, fieldOffset, comp;
-  PetscErrorCode ierr;
-
-  PetscFunctionBegin;
-  ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
-  ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
-  ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr);
-  ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr);
-  ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr);
-  ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr);
-  for (field = 0; field < numFields; ++field) {
-    numComponents += quad[field].numComponents;
-  }
-  ierr = DMPlexProjectFunctionLocal(dm, numComponents, funcs, INSERT_BC_VALUES, localX);CHKERRQ(ierr);
-  ierr = PetscMalloc4(dim,PetscReal,&coords,dim,PetscReal,&v0,dim*dim,PetscReal,&J,dim*dim,PetscReal,&invJ);CHKERRQ(ierr);
-  ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
-  for (c = cStart; c < cEnd; ++c) {
-    PetscScalar *x;
-    PetscReal    elemDiff = 0.0;
-
-    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);
-
-    for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) {
-      const PetscInt   numQuadPoints = quad[field].numQuadPoints;
-      const PetscReal *quadPoints    = quad[field].quadPoints;
-      const PetscReal *quadWeights   = quad[field].quadWeights;
-      const PetscInt   numBasisFuncs = quad[field].numBasisFuncs;
-      const PetscInt   numBasisComps = quad[field].numComponents;
-      const PetscReal *basis         = quad[field].basis;
-      PetscInt         q, d, e, fc, f;
-
-      if (debug) {
-        char title[1024];
-        ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr);
-        ierr = DMPrintCellVector(c, title, numBasisFuncs*numBasisComps, &x[fieldOffset]);CHKERRQ(ierr);
-      }
-      for (q = 0; q < numQuadPoints; ++q) {
-        for (d = 0; d < dim; d++) {
-          coords[d] = v0[d];
-          for (e = 0; e < dim; e++) {
-            coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0);
-          }
-        }
-        for (fc = 0; fc < numBasisComps; ++fc) {
-          const PetscReal funcVal     = PetscRealPart((*funcs[comp+fc])(coords));
-          PetscReal       interpolant = 0.0;
-          for (f = 0; f < numBasisFuncs; ++f) {
-            const PetscInt fidx = f*numBasisComps+fc;
-            interpolant += PetscRealPart(x[fieldOffset+fidx])*basis[q*numBasisFuncs*numBasisComps+fidx];
-          }
-          if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, "    elem %d field %d diff %g\n", c, field, PetscSqr(interpolant - funcVal)*quadWeights[q]*detJ);CHKERRQ(ierr);}
-          elemDiff += PetscSqr(interpolant - funcVal)*quadWeights[q]*detJ;
-        }
-      }
-      comp        += numBasisComps;
-      fieldOffset += numBasisFuncs*numBasisComps;
-    }
-    ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr);
-    if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, "  elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);}
-    localDiff += elemDiff;
-  }
-  ierr  = PetscFree4(coords,v0,J,invJ);CHKERRQ(ierr);
-  ierr  = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr);
-  ierr  = MPI_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPI_SUM, PETSC_COMM_WORLD);CHKERRQ(ierr);
-  *diff = PetscSqrtReal(*diff);
-  PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "DMPlexComputeResidualFEM"
-/*@
-  DMPlexComputeResidualFEM - Form the local residual F from the local input X using pointwise functions specified by the user
-
-  Input Parameters:
-+ dm - The mesh
-. X  - Local input vector
-- user - The user context
-
-  Output Parameter:
-. F  - Local output vector
-
-  Note:
-  The second member of the user context must be an FEMContext.
-
-  We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
-  like a GPU, or vectorize on a multicore machine.
-
-.seealso: DMPlexComputeJacobianActionFEM()
-*/
-PetscErrorCode DMPlexComputeResidualFEM(DM dm, Vec X, Vec F, void *user)
-{
-  DM_Plex         *mesh = (DM_Plex*) dm->data;
-  PetscFEM        *fem  = (PetscFEM*) &((DM*) user)[1];
-  PetscQuadrature *quad = fem->quad;
-  PetscSection     section;
-  PetscReal       *v0, *J, *invJ, *detJ;
-  PetscScalar     *elemVec, *u;
-  PetscInt         dim, numFields, field, numBatchesTmp = 1, numCells, cStart, cEnd, c;
-  PetscInt         cellDof = 0, numComponents = 0;
-  PetscErrorCode   ierr;
-
-  PetscFunctionBegin;
-  /* ierr = PetscLogEventBegin(ResidualFEMEvent,0,0,0,0);CHKERRQ(ierr); */
-  ierr     = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
-  ierr     = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
-  ierr     = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr);
-  ierr     = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
-  numCells = cEnd - cStart;
-  for (field = 0; field < numFields; ++field) {
-    cellDof       += quad[field].numBasisFuncs*quad[field].numComponents;
-    numComponents += quad[field].numComponents;
-  }
-  ierr = DMPlexProjectFunctionLocal(dm, numComponents, fem->bcFuncs, INSERT_BC_VALUES, X);CHKERRQ(ierr);
-  ierr = VecSet(F, 0.0);CHKERRQ(ierr);
-  ierr = PetscMalloc6(numCells*cellDof,PetscScalar,&u,numCells*dim,PetscReal,&v0,numCells*dim*dim,PetscReal,&J,numCells*dim*dim,PetscReal,&invJ,numCells,PetscReal,&detJ,numCells*cellDof,PetscScalar,&elemVec);CHKERRQ(ierr);
-  for (c = cStart; c < cEnd; ++c) {
-    PetscScalar *x;
-    PetscInt     i;
-
-    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, X, c, NULL, &x);CHKERRQ(ierr);
-
-    for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i];
-    ierr = DMPlexVecRestoreClosure(dm, NULL, X, c, NULL, &x);CHKERRQ(ierr);
-  }
-  for (field = 0; field < numFields; ++field) {
-    const PetscInt numQuadPoints = quad[field].numQuadPoints;
-    const PetscInt numBasisFuncs = quad[field].numBasisFuncs;
-    void           (*f0)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[]) = fem->f0Funcs[field];
-    void           (*f1)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f1[]) = fem->f1Funcs[field];
-    /* Conforming batches */
-    PetscInt blockSize  = numBasisFuncs*numQuadPoints;
-    PetscInt numBlocks  = 1;
-    PetscInt batchSize  = numBlocks * blockSize;
-    PetscInt numBatches = numBatchesTmp;
-    PetscInt numChunks  = numCells / (numBatches*batchSize);
-    /* Remainder */
-    PetscInt numRemainder = numCells % (numBatches * batchSize);
-    PetscInt offset       = numCells - numRemainder;
-
-    ierr = (*mesh->integrateResidualFEM)(numChunks*numBatches*batchSize, numFields, field, quad, u, v0, J, invJ, detJ, f0, f1, elemVec);CHKERRQ(ierr);
-    ierr = (*mesh->integrateResidualFEM)(numRemainder, numFields, field, quad, &u[offset*cellDof], &v0[offset*dim], &J[offset*dim*dim], &invJ[offset*dim*dim], &detJ[offset],
-                                         f0, f1, &elemVec[offset*cellDof]);CHKERRQ(ierr);
-  }
-  for (c = cStart; c < cEnd; ++c) {
-    if (mesh->printFEM > 1) {ierr = DMPrintCellVector(c, "Residual", cellDof, &elemVec[c*cellDof]);CHKERRQ(ierr);}
-    ierr = DMPlexVecSetClosure(dm, NULL, F, c, &elemVec[c*cellDof], ADD_VALUES);CHKERRQ(ierr);
-  }
-  ierr = PetscFree6(u,v0,J,invJ,detJ,elemVec);CHKERRQ(ierr);
-  if (mesh->printFEM) {
-    PetscMPIInt rank, numProcs;
-    PetscInt    p;
-
-    ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr);
-    ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm), &numProcs);CHKERRQ(ierr);
-    ierr = PetscPrintf(PETSC_COMM_WORLD, "Residual:\n");CHKERRQ(ierr);
-    for (p = 0; p < numProcs; ++p) {
-      if (p == rank) {
-        Vec f;
-
-        ierr = VecDuplicate(F, &f);CHKERRQ(ierr);
-        ierr = VecCopy(F, f);CHKERRQ(ierr);
-        ierr = VecChop(f, 1.0e-10);CHKERRQ(ierr);
-        ierr = VecView(f, PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);
-        ierr = VecDestroy(&f);CHKERRQ(ierr);
-        ierr = PetscViewerFlush(PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);
-      }
-      ierr = PetscBarrier((PetscObject) dm);CHKERRQ(ierr);
-    }
-  }
-  /* ierr = PetscLogEventEnd(ResidualFEMEvent,0,0,0,0);CHKERRQ(ierr); */
-  PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "DMPlexComputeJacobianActionFEM"
-/*@C
-  DMPlexComputeJacobianActionFEM - Form the local action of Jacobian J(u) on the local input X using pointwise functions specified by the user
-
-  Input Parameters:
-+ dm - The mesh
-. J  - The Jacobian shell matrix
-. X  - Local input vector
-- user - The user context
-
-  Output Parameter:
-. F  - Local output vector
-
-  Note:
-  The second member of the user context must be an FEMContext.
-
-  We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
-  like a GPU, or vectorize on a multicore machine.
-
-.seealso: DMPlexComputeResidualFEM()
-*/
-PetscErrorCode DMPlexComputeJacobianActionFEM(DM dm, Mat Jac, Vec X, Vec F, void *user)
-{
-  DM_Plex         *mesh = (DM_Plex*) dm->data;
-  PetscFEM        *fem  = (PetscFEM*) &((DM*) user)[1];
-  PetscQuadrature *quad = fem->quad;
-  PetscSection     section;
-  JacActionCtx    *jctx;
-  PetscReal       *v0, *J, *invJ, *detJ;
-  PetscScalar     *elemVec, *u, *a;
-  PetscInt         dim, numFields, field, numBatchesTmp = 1, numCells, cStart, cEnd, c;
-  PetscInt         cellDof = 0;
-  PetscErrorCode   ierr;
-
-  PetscFunctionBegin;
-  /* ierr = PetscLogEventBegin(JacobianActionFEMEvent,0,0,0,0);CHKERRQ(ierr); */
-  ierr     = MatShellGetContext(Jac, &jctx);CHKERRQ(ierr);
-  ierr     = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
-  ierr     = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
-  ierr     = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr);
-  ierr     = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
-  numCells = cEnd - cStart;
-  for (field = 0; field < numFields; ++field) {
-    cellDof += quad[field].numBasisFuncs*quad[field].numComponents;
-  }
-  ierr = VecSet(F, 0.0);CHKERRQ(ierr);
-  ierr = PetscMalloc7(numCells*cellDof,PetscScalar,&u,numCells*cellDof,PetscScalar,&a,numCells*dim,PetscReal,&v0,numCells*dim*dim,PetscReal,&J,numCells*dim*dim,PetscReal,&invJ,numCells,PetscReal,&detJ,numCells*cellDof,PetscScalar,&elemVec);CHKERRQ(ierr);
-  for (c = cStart; c < cEnd; ++c) {
-    PetscScalar *x;
-    PetscInt     i;
-
-    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);
-    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);
-    for (i = 0; i < cellDof; ++i) a[c*cellDof+i] = x[i];
-    ierr = DMPlexVecRestoreClosure(dm, NULL, X, c, NULL, &x);CHKERRQ(ierr);
-  }
-  for (field = 0; field < numFields; ++field) {
-    const PetscInt numQuadPoints = quad[field].numQuadPoints;
-    const PetscInt numBasisFuncs = quad[field].numBasisFuncs;
-    /* Conforming batches */
-    PetscInt blockSize  = numBasisFuncs*numQuadPoints;
-    PetscInt numBlocks  = 1;
-    PetscInt batchSize  = numBlocks * blockSize;
-    PetscInt numBatches = numBatchesTmp;
-    PetscInt numChunks  = numCells / (numBatches*batchSize);
-    /* Remainder */
-    PetscInt numRemainder = numCells % (numBatches * batchSize);
-    PetscInt offset       = numCells - numRemainder;
-
-    ierr = (*mesh->integrateJacobianActionFEM)(numChunks*numBatches*batchSize, numFields, field, quad, u, a, v0, J, invJ, detJ, fem->g0Funcs, fem->g1Funcs, fem->g2Funcs, fem->g3Funcs, elemVec);CHKERRQ(ierr);
-    ierr = (*mesh->integrateJacobianActionFEM)(numRemainder, numFields, field, quad, &u[offset*cellDof], &a[offset*cellDof], &v0[offset*dim], &J[offset*dim*dim], &invJ[offset*dim*dim], &detJ[offset],
-                                               fem->g0Funcs, fem->g1Funcs, fem->g2Funcs, fem->g3Funcs, &elemVec[offset*cellDof]);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 = PetscFree7(u,a,v0,J,invJ,detJ,elemVec);CHKERRQ(ierr);
-  if (mesh->printFEM) {
-    PetscMPIInt rank, numProcs;
-    PetscInt    p;
-
-    ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr);
-    ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm), &numProcs);CHKERRQ(ierr);
-    ierr = PetscPrintf(PETSC_COMM_WORLD, "Jacobian Action:\n");CHKERRQ(ierr);
-    for (p = 0; p < numProcs; ++p) {
-      if (p == rank) {ierr = VecView(F, PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);}
-      ierr = PetscBarrier((PetscObject) dm);CHKERRQ(ierr);
-    }
-  }
-  /* ierr = PetscLogEventEnd(JacobianActionFEMEvent,0,0,0,0);CHKERRQ(ierr); */
-  PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "DMPlexComputeJacobianFEM"
-/*@
-  DMPlexComputeJacobianFEM - Form the local portion of the Jacobian matrix J at the local solution X using pointwise functions specified by the user.
-
-  Input Parameters:
-+ dm - The mesh
-. X  - Local input vector
-- user - The user context
-
-  Output Parameter:
-. Jac  - Jacobian matrix
-
-  Note:
-  The second member of the user context must be an FEMContext.
-
-  We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
-  like a GPU, or vectorize on a multicore machine.
-
-.seealso: FormFunctionLocal()
-*/
-PetscErrorCode DMPlexComputeJacobianFEM(DM dm, Vec X, Mat Jac, Mat JacP, MatStructure *str,void *user)
-{
-  DM_Plex         *mesh = (DM_Plex*) dm->data;
-  PetscFEM        *fem  = (PetscFEM*) &((DM*) user)[1];
-  PetscQuadrature *quad = fem->quad;
-  PetscSection     section;
-  PetscReal       *v0, *J, *invJ, *detJ;
-  PetscScalar     *elemMat, *u;
-  PetscInt         dim, numFields, field, fieldI, numBatchesTmp = 1, numCells, cStart, cEnd, c;
-  PetscInt         cellDof = 0, numComponents = 0;
-  PetscBool        isShell;
-  PetscErrorCode   ierr;
-
-  PetscFunctionBegin;
-  /* ierr = PetscLogEventBegin(JacobianFEMEvent,0,0,0,0);CHKERRQ(ierr); */
-  ierr     = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
-  ierr     = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
-  ierr     = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr);
-  ierr     = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
-  numCells = cEnd - cStart;
-  for (field = 0; field < numFields; ++field) {
-    cellDof       += quad[field].numBasisFuncs*quad[field].numComponents;
-    numComponents += quad[field].numComponents;
-  }
-  ierr = DMPlexProjectFunctionLocal(dm, numComponents, fem->bcFuncs, INSERT_BC_VALUES, X);CHKERRQ(ierr);
-  ierr = MatZeroEntries(JacP);CHKERRQ(ierr);
-  ierr = PetscMalloc6(numCells*cellDof,PetscScalar,&u,numCells*dim,PetscReal,&v0,numCells*dim*dim,PetscReal,&J,numCells*dim*dim,PetscReal,&invJ,numCells,PetscReal,&detJ,numCells*cellDof*cellDof,PetscScalar,&elemMat);CHKERRQ(ierr);
-  for (c = cStart; c < cEnd; ++c) {
-    PetscScalar *x;
-    PetscInt     i;
-
-    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, X, c, NULL, &x);CHKERRQ(ierr);
-
-    for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i];
-    ierr = DMPlexVecRestoreClosure(dm, NULL, X, c, NULL, &x);CHKERRQ(ierr);
-  }
-  ierr = PetscMemzero(elemMat, numCells*cellDof*cellDof * sizeof(PetscScalar));CHKERRQ(ierr);
-  for (fieldI = 0; fieldI < numFields; ++fieldI) {
-    const PetscInt numQuadPoints = quad[fieldI].numQuadPoints;
-    const PetscInt numBasisFuncs = quad[fieldI].numBasisFuncs;
-    PetscInt       fieldJ;
-
-    for (fieldJ = 0; fieldJ < numFields; ++fieldJ) {
-      void (*g0)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar g0[]) = fem->g0Funcs[fieldI*numFields+fieldJ];
-      void (*g1)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar g1[]) = fem->g1Funcs[fieldI*numFields+fieldJ];
-      void (*g2)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar g2[]) = fem->g2Funcs[fieldI*numFields+fieldJ];
-      void (*g3)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar g3[]) = fem->g3Funcs[fieldI*numFields+fieldJ];
-      /* Conforming batches */
-      PetscInt blockSize  = numBasisFuncs*numQuadPoints;
-      PetscInt numBlocks  = 1;
-      PetscInt batchSize  = numBlocks * blockSize;
-      PetscInt numBatches = numBatchesTmp;
-      PetscInt numChunks  = numCells / (numBatches*batchSize);
-      /* Remainder */
-      PetscInt numRemainder = numCells % (numBatches * batchSize);
-      PetscInt offset       = numCells - numRemainder;
-
-      ierr = (*mesh->integrateJacobianFEM)(numChunks*numBatches*batchSize, numFields, fieldI, fieldJ, quad, u, v0, J, invJ, detJ, g0, g1, g2, g3, elemMat);CHKERRQ(ierr);
-      ierr = (*mesh->integrateJacobianFEM)(numRemainder, numFields, fieldI, fieldJ, quad, &u[offset*cellDof], &v0[offset*dim], &J[offset*dim*dim], &invJ[offset*dim*dim], &detJ[offset],
-                                           g0, g1, g2, g3, &elemMat[offset*cellDof*cellDof]);CHKERRQ(ierr);
-    }
-  }
-  for (c = cStart; c < cEnd; ++c) {
-    if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(c, "Jacobian", cellDof, cellDof, &elemMat[c*cellDof*cellDof]);CHKERRQ(ierr);}
-    ierr = DMPlexMatSetClosure(dm, NULL, NULL, JacP, c, &elemMat[c*cellDof*cellDof], ADD_VALUES);CHKERRQ(ierr);
-  }
-  ierr = PetscFree6(u,v0,J,invJ,detJ,elemMat);CHKERRQ(ierr);
-
-  /* Assemble matrix, using the 2-step process:
-       MatAssemblyBegin(), MatAssemblyEnd(). */
-  ierr = MatAssemblyBegin(JacP, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
-  ierr = MatAssemblyEnd(JacP, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
-
-  if (mesh->printFEM) {
-    ierr = PetscPrintf(PETSC_COMM_WORLD, "Jacobian:\n");CHKERRQ(ierr);
-    ierr = MatChop(JacP, 1.0e-10);CHKERRQ(ierr);
-    ierr = MatView(JacP, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
-  }
-  /* ierr = PetscLogEventEnd(JacobianFEMEvent,0,0,0,0);CHKERRQ(ierr); */
-  ierr = PetscObjectTypeCompare((PetscObject)Jac, MATSHELL, &isShell);CHKERRQ(ierr);
-  if (isShell) {
-    JacActionCtx *jctx;
-
-    ierr = MatShellGetContext(Jac, &jctx);CHKERRQ(ierr);
-    ierr = VecCopy(X, jctx->u);CHKERRQ(ierr);
-  }
-  *str = SAME_NONZERO_PATTERN;
-  PetscFunctionReturn(0);
-}
-
 
 #undef __FUNCT__
 #define __FUNCT__ "PetscSectionCreateGlobalSectionLabel"

src/dm/impls/plex/plexfem.c

+#include <petsc-private/dmpleximpl.h>   /*I      "petscdmplex.h"   I*/
+
+#undef __FUNCT__
+#define __FUNCT__ "DMPlexGetScale"
+PetscErrorCode DMPlexGetScale(DM dm, PetscUnit unit, PetscReal *scale)
+{
+  DM_Plex *mesh = (DM_Plex*) dm->data;
+
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
+  PetscValidPointer(scale, 3);
+  *scale = mesh->scale[unit];
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "DMPlexSetScale"
+PetscErrorCode DMPlexSetScale(DM dm, PetscUnit unit, PetscReal scale)
+{
+  DM_Plex *mesh = (DM_Plex*) dm->data;
+
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
+  mesh->scale[unit] = scale;
+  PetscFunctionReturn(0);
+}
+
+PETSC_STATIC_INLINE PetscInt epsilon(PetscInt i, PetscInt j, PetscInt k)
+{
+  switch (i) {
+  case 0:
+    switch (j) {
+    case 0: return 0;
+    case 1:
+      switch (k) {
+      case 0: return 0;
+      case 1: return 0;
+      case 2: return 1;
+      }
+    case 2:
+      switch (k) {
+      case 0: return 0;
+      case 1: return -1;
+      case 2: return 0;
+      }
+    }
+  case 1:
+    switch (j) {
+    case 0:
+      switch (k) {
+      case 0: return 0;
+      case 1: return 0;
+      case 2: return -1;
+      }
+    case 1: return 0;
+    case 2:
+      switch (k) {
+      case 0: return 1;
+      case 1: return 0;
+      case 2: return 0;
+      }
+    }
+  case 2:
+    switch (j) {
+    case 0:
+      switch (k) {
+      case 0: return 0;
+      case 1: return 1;
+      case 2: return 0;
+      }
+    case 1:
+      switch (k) {
+      case 0: return -1;
+      case 1: return 0;
+      case 2: return 0;
+      }
+    case 2: return 0;
+    }
+  }
+  return 0;
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "DMPlexCreateRigidBody"
+/*@C
+  DMPlexCreateRigidBody - create rigid body modes from coordinates
+
+  Collective on DM
+
+  Input Arguments:
++ dm - the DM
+. section - the local section associated with the rigid field, or NULL for the default section
+- globalSection - the global section associated with the rigid field, or NULL for the default section
+
+  Output Argument:
+. sp - the null space
+
+  Note: This is necessary to take account of Dirichlet conditions on the displacements
+
+  Level: advanced
+
+.seealso: MatNullSpaceCreate()
+@*/
+PetscErrorCode DMPlexCreateRigidBody(DM dm, PetscSection section, PetscSection globalSection, MatNullSpace *sp)
+{
+  MPI_Comm       comm;
+  Vec            coordinates, localMode, mode[6];
+  PetscSection   coordSection;
+  PetscScalar   *coords;
+  PetscInt       dim, vStart, vEnd, v, n, m, d, i, j;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
+  ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
+  if (dim == 1) {
+    ierr = MatNullSpaceCreate(comm, PETSC_TRUE, 0, NULL, sp);CHKERRQ(ierr);
+    PetscFunctionReturn(0);
+  }
+  if (!section)       {ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);}
+  if (!globalSection) {ierr = DMGetDefaultGlobalSection(dm, &globalSection);CHKERRQ(ierr);}
+  ierr = PetscSectionGetConstrainedStorageSize(globalSection, &n);CHKERRQ(ierr);
+  ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
+  ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
+  ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
+  m    = (dim*(dim+1))/2;
+  ierr = VecCreate(comm, &mode[0]);CHKERRQ(ierr);
+  ierr = VecSetSizes(mode[0], n, PETSC_DETERMINE);CHKERRQ(ierr);
+  ierr = VecSetUp(mode[0]);CHKERRQ(ierr);
+  for (i = 1; i < m; ++i) {ierr = VecDuplicate(mode[0], &mode[i]);CHKERRQ(ierr);}
+  /* Assume P1 */
+  ierr = DMGetLocalVector(dm, &localMode);CHKERRQ(ierr);
+  for (d = 0; d < dim; ++d) {
+    PetscScalar values[3] = {0.0, 0.0, 0.0};
+
+    values[d] = 1.0;
+    ierr      = VecSet(localMode, 0.0);CHKERRQ(ierr);
+    for (v = vStart; v < vEnd; ++v) {
+      ierr = DMPlexVecSetClosure(dm, section, localMode, v, values, INSERT_VALUES);CHKERRQ(ierr);
+    }
+    ierr = DMLocalToGlobalBegin(dm, localMode, INSERT_VALUES, mode[d]);CHKERRQ(ierr);
+    ierr = DMLocalToGlobalEnd(dm, localMode, INSERT_VALUES, mode[d]);CHKERRQ(ierr);
+  }
+  ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr);
+  for (d = dim; d < dim*(dim+1)/2; ++d) {
+    PetscInt i, j, k = dim > 2 ? d - dim : d;
+
+    ierr = VecSet(localMode, 0.0);CHKERRQ(ierr);
+    for (v = vStart; v < vEnd; ++v) {
+      PetscScalar values[3] = {0.0, 0.0, 0.0};
+      PetscInt    off;
+
+      ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr);
+      for (i = 0; i < dim; ++i) {
+        for (j = 0; j < dim; ++j) {
+          values[j] += epsilon(i, j, k)*PetscRealPart(coords[off+i]);
+        }
+      }
+      ierr = DMPlexVecSetClosure(dm, section, localMode, v, values, INSERT_VALUES);CHKERRQ(ierr);
+    }
+    ierr = DMLocalToGlobalBegin(dm, localMode, INSERT_VALUES, mode[d]);CHKERRQ(ierr);
+    ierr = DMLocalToGlobalEnd(dm, localMode, INSERT_VALUES, mode[d]);CHKERRQ(ierr);
+  }
+  ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr);
+  ierr = DMRestoreLocalVector(dm, &localMode);CHKERRQ(ierr);
+  for (i = 0; i < dim; ++i) {ierr = VecNormalize(mode[i], NULL);CHKERRQ(ierr);}
+  /* Orthonormalize system */
+  for (i = dim; i < m; ++i) {
+    PetscScalar dots[6];
+
+    ierr = VecMDot(mode[i], i, mode, dots);CHKERRQ(ierr);
+    for (j = 0; j < i; ++j) dots[j] *= -1.0;
+    ierr = VecMAXPY(mode[i], i, dots, mode);CHKERRQ(ierr);
+    ierr = VecNormalize(mode[i], NULL);CHKERRQ(ierr);
+  }
+  ierr = MatNullSpaceCreate(comm, PETSC_FALSE, m, mode, sp);CHKERRQ(ierr);
+  for (i = 0; i< m; ++i) {ierr = VecDestroy(&mode[i]);CHKERRQ(ierr);}
+  PetscFunctionReturn(0);
+}
+/*******************************************************************************
+This should be in a separate Discretization object, but I am not sure how to lay
+it out yet, so I am stuffing things here while I experiment.
+*******************************************************************************/
+#undef __FUNCT__
+#define __FUNCT__ "DMPlexSetFEMIntegration"
+PetscErrorCode DMPlexSetFEMIntegration(DM dm,
+                                          PetscErrorCode (*integrateResidualFEM)(PetscInt, PetscInt, PetscInt, PetscQuadrature[], const PetscScalar[],
+                                                                                 const PetscReal[], const PetscReal[], const PetscReal[], const PetscReal[],
+                                                                                 void (*)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]),
+                                                                                 void (*)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]), PetscScalar[]),
+                                          PetscErrorCode (*integrateJacobianActionFEM)(PetscInt, PetscInt, PetscInt, PetscQuadrature[], const PetscScalar[], const PetscScalar[],
+                                                                                       const PetscReal[], const PetscReal[], const PetscReal[], const PetscReal[],
+                                                                                       void (**)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]),
+                                                                                       void (**)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]),
+                                                                                       void (**)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]),
+                                                                                       void (**)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]), PetscScalar[]),
+                                          PetscErrorCode (*integrateJacobianFEM)(PetscInt, PetscInt, PetscInt, PetscInt, PetscQuadrature[], const PetscScalar[],
+                                                                                 const PetscReal[], const PetscReal[], const PetscReal[], const PetscReal[],
+                                                                                 void (*)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]),
+                                                                                 void (*)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]),
+                                                                                 void (*)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]),
+                                                                                 void (*)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]), PetscScalar[]))
+{
+  DM_Plex *mesh = (DM_Plex*) dm->data;
+
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
+  mesh->integrateResidualFEM       = integrateResidualFEM;
+  mesh->integrateJacobianActionFEM = integrateJacobianActionFEM;
+  mesh->integrateJacobianFEM       = integrateJacobianFEM;
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "DMPlexProjectFunctionLocal"
+PetscErrorCode DMPlexProjectFunctionLocal(DM dm, PetscInt numComp, PetscScalar (**funcs)(const PetscReal []), InsertMode mode, Vec localX)
+{
+  Vec            coordinates;
+  PetscSection   section, cSection;
+  PetscInt       dim, vStart, vEnd, v, c, d;
+  PetscScalar   *values, *cArray;
+  PetscReal     *coords;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
+  ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
+  ierr = DMPlexGetCoordinateSection(dm, &cSection);CHKERRQ(ierr);
+  ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
+  ierr = PetscMalloc(numComp * sizeof(PetscScalar), &values);CHKERRQ(ierr);
+  ierr = VecGetArray(coordinates, &cArray);CHKERRQ(ierr);
+  ierr = PetscSectionGetDof(cSection, vStart, &dim);CHKERRQ(ierr);
+  ierr = PetscMalloc(dim * sizeof(PetscReal),&coords);CHKERRQ(ierr);
+  for (v = vStart; v < vEnd; ++v) {
+    PetscInt dof, off;
+
+    ierr = PetscSectionGetDof(cSection, v, &dof);CHKERRQ(ierr);
+    ierr = PetscSectionGetOffset(cSection, v, &off);CHKERRQ(ierr);
+    if (dof > dim) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Cannot have more coordinates %d then dimensions %d", dof, dim);
+    for (d = 0; d < dof; ++d) coords[d] = PetscRealPart(cArray[off+d]);
+    for (c = 0; c < numComp; ++c) values[c] = (*funcs[c])(coords);
+    ierr = VecSetValuesSection(localX, section, v, values, mode);CHKERRQ(ierr);
+  }
+  ierr = VecRestoreArray(coordinates, &cArray);CHKERRQ(ierr);
+  /* Temporary, must be replaced by a projection on the finite element basis */
+  {
+    PetscInt eStart = 0, eEnd = 0, e, depth;
+
+    ierr = DMPlexGetLabelSize(dm, "depth", &depth);CHKERRQ(ierr);
+    --depth;
+    if (depth > 1) {ierr = DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);CHKERRQ(ierr);}
+    for (e = eStart; e < eEnd; ++e) {
+      const PetscInt *cone = NULL;
+      PetscInt        coneSize, d;
+      PetscScalar    *coordsA, *coordsB;
+
+      ierr = DMPlexGetConeSize(dm, e, &coneSize);CHKERRQ(ierr);
+      ierr = DMPlexGetCone(dm, e, &cone);CHKERRQ(ierr);
+      if (coneSize != 2) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_SIZ, "Cone size %d for point %d should be 2", coneSize, e);
+      ierr = VecGetValuesSection(coordinates, cSection, cone[0], &coordsA);CHKERRQ(ierr);
+      ierr = VecGetValuesSection(coordinates, cSection, cone[1], &coordsB);CHKERRQ(ierr);
+      for (d = 0; d < dim; ++d) {
+        coords[d] = 0.5*(PetscRealPart(coordsA[d]) + PetscRealPart(coordsB[d]));
+      }
+      for (c = 0; c < numComp; ++c) values[c] = (*funcs[c])(coords);
+      ierr = VecSetValuesSection(localX, section, e, values, mode);CHKERRQ(ierr);
+    }
+  }
+
+  ierr = PetscFree(coords);CHKERRQ(ierr);
+  ierr = PetscFree(values);CHKERRQ(ierr);
+#if 0
+  const PetscInt localDof = this->_mesh->sizeWithBC(s, *cells->begin());
+  PetscReal      detJ;
+
+  ierr = PetscMalloc(localDof * sizeof(PetscScalar), &values);CHKERRQ(ierr);
+  ierr = PetscMalloc2(dim,PetscReal,&v0,dim*dim,PetscReal,&J);CHKERRQ(ierr);
+  ALE::ISieveVisitor::PointRetriever<PETSC_MESH_TYPE::sieve_type> pV(PetscPowInt(this->_mesh->getSieve()->getMaxConeSize(),dim+1), true);
+
+  for (PetscInt c = cStart; c < cEnd; ++c) {
+    ALE::ISieveTraversal<PETSC_MESH_TYPE::sieve_type>::orientedClosure(*this->_mesh->getSieve(), c, pV);
+    const PETSC_MESH_TYPE::point_type *oPoints = pV.getPoints();
+    const int                          oSize   = pV.getSize();
+    int                                v       = 0;
+
+    ierr = DMPlexComputeCellGeometry(dm, c, v0, J, NULL, &detJ);CHKERRQ(ierr);
+    for (PetscInt cl = 0; cl < oSize; ++cl) {
+      const PetscInt fDim;
+
+      ierr = PetscSectionGetDof(oPoints[cl], &fDim);CHKERRQ(ierr);
+      if (pointDim) {
+        for (PetscInt d = 0; d < fDim; ++d, ++v) {
+          values[v] = (*this->_options.integrate)(v0, J, v, initFunc);
+        }
+      }
+    }
+    ierr = DMPlexVecSetClosure(dm, NULL, localX, c, values);CHKERRQ(ierr);
+    pV.clear();
+  }
+  ierr = PetscFree2(v0,J);CHKERRQ(ierr);
+  ierr = PetscFree(values);CHKERRQ(ierr);
+#endif
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "DMPlexProjectFunction"
+/*@C
+  DMPlexProjectFunction - This projects the given function into the function space provided.
+
+  Input Parameters:
++ dm      - The DM
+. numComp - The number of components (functions)
+. funcs   - The coordinate functions to evaluate
+- mode    - The insertion mode for values
+
+  Output Parameter:
+. X - vector
+
+  Level: developer
+
+  Note:
+  This currently just calls the function with the coordinates of each vertex and edge midpoint, and stores the result in a vector.
+  We will eventually fix it.
+
+,seealso: DMPlexComputeL2Diff()
+*/
+PetscErrorCode DMPlexProjectFunction(DM dm, PetscInt numComp, PetscScalar (**funcs)(const PetscReal []), InsertMode mode, Vec X)
+{
+  Vec            localX;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr);
+  ierr = DMPlexProjectFunctionLocal(dm, numComp, funcs, mode, localX);CHKERRQ(ierr);
+  ierr = DMLocalToGlobalBegin(dm, localX, mode, X);CHKERRQ(ierr);
+  ierr = DMLocalToGlobalEnd(dm, localX, mode, X);CHKERRQ(ierr);
+  ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "DMPlexComputeL2Diff"
+/*@C
+  DMPlexComputeL2Diff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h.
+
+  Input Parameters:
++ dm    - The DM
+. quad  - The PetscQuadrature object for each field
+. funcs - The functions to evaluate for each field component
+- X     - The coefficient vector u_h
+
+  Output Parameter:
+. diff - The diff ||u - u_h||_2
+
+  Level: developer
+
+.seealso: DMPlexProjectFunction()
+*/
+PetscErrorCode DMPlexComputeL2Diff(DM dm, PetscQuadrature quad[], PetscScalar (**funcs)(const PetscReal []), Vec X, PetscReal *diff)
+{
+  const PetscInt debug = 0;
+  PetscSection   section;
+  Vec            localX;
+  PetscReal     *coords, *v0, *J, *invJ, detJ;
+  PetscReal      localDiff = 0.0;
+  PetscInt       dim, numFields, numComponents = 0, cStart, cEnd, c, field, fieldOffset, comp;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
+  ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
+  ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr);
+  ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr);
+  ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr);
+  ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr);
+  for (field = 0; field < numFields; ++field) {
+    numComponents += quad[field].numComponents;
+  }
+  ierr = DMPlexProjectFunctionLocal(dm, numComponents, funcs, INSERT_BC_VALUES, localX);CHKERRQ(ierr);
+  ierr = PetscMalloc4(dim,PetscReal,&coords,dim,PetscReal,&v0,dim*dim,PetscReal,&J,dim*dim,PetscReal,&invJ);CHKERRQ(ierr);
+  ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
+  for (c = cStart; c < cEnd; ++c) {
+    PetscScalar *x;
+    PetscReal    elemDiff = 0.0;
+
+    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);
+
+    for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) {
+      const PetscInt   numQuadPoints = quad[field].numQuadPoints;
+      const PetscReal *quadPoints    = quad[field].quadPoints;
+      const PetscReal *quadWeights   = quad[field].quadWeights;
+      const PetscInt   numBasisFuncs = quad[field].numBasisFuncs;
+      const PetscInt   numBasisComps = quad[field].numComponents;
+      const PetscReal *basis         = quad[field].basis;
+      PetscInt         q, d, e, fc, f;
+
+      if (debug) {
+        char title[1024];
+        ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr);
+        ierr = DMPrintCellVector(c, title, numBasisFuncs*numBasisComps, &x[fieldOffset]);CHKERRQ(ierr);
+      }
+      for (q = 0; q < numQuadPoints; ++q) {
+        for (d = 0; d < dim; d++) {
+          coords[d] = v0[d];
+          for (e = 0; e < dim; e++) {
+            coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0);
+          }
+        }
+        for (fc = 0; fc < numBasisComps; ++fc) {
+          const PetscReal funcVal     = PetscRealPart((*funcs[comp+fc])(coords));
+          PetscReal       interpolant = 0.0;
+          for (f = 0; f < numBasisFuncs; ++f) {
+            const PetscInt fidx = f*numBasisComps+fc;
+            interpolant += PetscRealPart(x[fieldOffset+fidx])*basis[q*numBasisFuncs*numBasisComps+fidx];
+          }
+          if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, "    elem %d field %d diff %g\n", c, field, PetscSqr(interpolant - funcVal)*quadWeights[q]*detJ);CHKERRQ(ierr);}
+          elemDiff += PetscSqr(interpolant - funcVal)*quadWeights[q]*detJ;
+        }
+      }
+      comp        += numBasisComps;
+      fieldOffset += numBasisFuncs*numBasisComps;
+    }
+    ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr);
+    if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, "  elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);}
+    localDiff += elemDiff;
+  }
+  ierr  = PetscFree4(coords,v0,J,invJ);CHKERRQ(ierr);
+  ierr  = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr);
+  ierr  = MPI_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPI_SUM, PETSC_COMM_WORLD);CHKERRQ(ierr);
+  *diff = PetscSqrtReal(*diff);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "DMPlexComputeResidualFEM"
+/*@
+  DMPlexComputeResidualFEM - Form the local residual F from the local input X using pointwise functions specified by the user
+
+  Input Parameters:
++ dm - The mesh
+. X  - Local input vector
+- user - The user context
+
+  Output Parameter:
+. F  - Local output vector
+
+  Note:
+  The second member of the user context must be an FEMContext.
+
+  We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
+  like a GPU, or vectorize on a multicore machine.
+
+.seealso: DMPlexComputeJacobianActionFEM()
+*/
+PetscErrorCode DMPlexComputeResidualFEM(DM dm, Vec X, Vec F, void *user)
+{
+  DM_Plex         *mesh = (DM_Plex*) dm->data;
+  PetscFEM        *fem  = (PetscFEM*) &((DM*) user)[1];
+  PetscQuadrature *quad = fem->quad;
+  PetscSection     section;
+  PetscReal       *v0, *J, *invJ, *detJ;
+  PetscScalar     *elemVec, *u;
+  PetscInt         dim, numFields, field, numBatchesTmp = 1, numCells, cStart, cEnd, c;
+  PetscInt         cellDof = 0, numComponents = 0;
+  PetscErrorCode   ierr;
+
+  PetscFunctionBegin;
+  /* ierr = PetscLogEventBegin(ResidualFEMEvent,0,0,0,0);CHKERRQ(ierr); */
+  ierr     = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
+  ierr     = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
+  ierr     = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr);
+  ierr     = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
+  numCells = cEnd - cStart;
+  for (field = 0; field < numFields; ++field) {
+    cellDof       += quad[field].numBasisFuncs*quad[field].numComponents;
+    numComponents += quad[field].numComponents;
+  }
+  ierr = DMPlexProjectFunctionLocal(dm, numComponents, fem->bcFuncs, INSERT_BC_VALUES, X);CHKERRQ(ierr);
+  ierr = VecSet(F, 0.0);CHKERRQ(ierr);
+  ierr = PetscMalloc6(numCells*cellDof,PetscScalar,&u,numCells*dim,PetscReal,&v0,numCells*dim*dim,PetscReal,&J,numCells*dim*dim,PetscReal,&invJ,numCells,PetscReal,&detJ,numCells*cellDof,PetscScalar,&elemVec);CHKERRQ(ierr);
+  for (c = cStart; c < cEnd; ++c) {
+    PetscScalar *x;
+    PetscInt     i;
+
+    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, X, c, NULL, &x);CHKERRQ(ierr);
+
+    for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i];
+    ierr = DMPlexVecRestoreClosure(dm, NULL, X, c, NULL, &x);CHKERRQ(ierr);
+  }
+  for (field = 0; field < numFields; ++field) {
+    const PetscInt numQuadPoints = quad[field].numQuadPoints;
+    const PetscInt numBasisFuncs = quad[field].numBasisFuncs;
+    void           (*f0)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[]) = fem->f0Funcs[field];
+    void           (*f1)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f1[]) = fem->f1Funcs[field];
+    /* Conforming batches */
+    PetscInt blockSize  = numBasisFuncs*numQuadPoints;
+    PetscInt numBlocks  = 1;
+    PetscInt batchSize  = numBlocks * blockSize;
+    PetscInt numBatches = numBatchesTmp;
+    PetscInt numChunks  = numCells / (numBatches*batchSize);
+    /* Remainder */
+    PetscInt numRemainder = numCells % (numBatches * batchSize);
+    PetscInt offset       = numCells - numRemainder;
+
+    ierr = (*mesh->integrateResidualFEM)(numChunks*numBatches*batchSize, numFields, field, quad, u, v0, J, invJ, detJ, f0, f1, elemVec);CHKERRQ(ierr);
+    ierr = (*mesh->integrateResidualFEM)(numRemainder, numFields, field, quad, &u[offset*cellDof], &v0[offset*dim], &J[offset*dim*dim], &invJ[offset*dim*dim], &detJ[offset],
+                                         f0, f1, &elemVec[offset*cellDof]);CHKERRQ(ierr);
+  }
+  for (c = cStart; c < cEnd; ++c) {
+    if (mesh->printFEM > 1) {ierr = DMPrintCellVector(c, "Residual", cellDof, &elemVec[c*cellDof]);CHKERRQ(ierr);}
+    ierr = DMPlexVecSetClosure(dm, NULL, F, c, &elemVec[c*cellDof], ADD_VALUES);CHKERRQ(ierr);
+  }
+  ierr = PetscFree6(u,v0,J,invJ,detJ,elemVec);CHKERRQ(ierr);
+  if (mesh->printFEM) {
+    PetscMPIInt rank, numProcs;
+    PetscInt    p;
+
+    ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr);
+    ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm), &numProcs);CHKERRQ(ierr);
+    ierr = PetscPrintf(PETSC_COMM_WORLD, "Residual:\n");CHKERRQ(ierr);
+    for (p = 0; p < numProcs; ++p) {
+      if (p == rank) {
+        Vec f;
+
+        ierr = VecDuplicate(F, &f);CHKERRQ(ierr);
+        ierr = VecCopy(F, f);CHKERRQ(ierr);
+        ierr = VecChop(f, 1.0e-10);CHKERRQ(ierr);
+        ierr = VecView(f, PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);
+        ierr = VecDestroy(&f);CHKERRQ(ierr);
+        ierr = PetscViewerFlush(PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);
+      }
+      ierr = PetscBarrier((PetscObject) dm);CHKERRQ(ierr);
+    }
+  }
+  /* ierr = PetscLogEventEnd(ResidualFEMEvent,0,0,0,0);CHKERRQ(ierr); */
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "DMPlexComputeJacobianActionFEM"
+/*@C
+  DMPlexComputeJacobianActionFEM - Form the local action of Jacobian J(u) on the local input X using pointwise functions specified by the user
+
+  Input Parameters:
++ dm - The mesh
+. J  - The Jacobian shell matrix
+. X  - Local input vector
+- user - The user context
+
+  Output Parameter:
+. F  - Local output vector
+
+  Note:
+  The second member of the user context must be an FEMContext.
+
+  We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
+  like a GPU, or vectorize on a multicore machine.
+
+.seealso: DMPlexComputeResidualFEM()
+*/
+PetscErrorCode DMPlexComputeJacobianActionFEM(DM dm, Mat Jac, Vec X, Vec F, void *user)
+{
+  DM_Plex         *mesh = (DM_Plex*) dm->data;
+  PetscFEM        *fem  = (PetscFEM*) &((DM*) user)[1];
+  PetscQuadrature *quad = fem->quad;
+  PetscSection     section;
+  JacActionCtx    *jctx;
+  PetscReal       *v0, *J, *invJ, *detJ;
+  PetscScalar     *elemVec, *u, *a;
+  PetscInt         dim, numFields, field, numBatchesTmp = 1, numCells, cStart, cEnd, c;
+  PetscInt         cellDof = 0;
+  PetscErrorCode   ierr;
+
+  PetscFunctionBegin;
+  /* ierr = PetscLogEventBegin(JacobianActionFEMEvent,0,0,0,0);CHKERRQ(ierr); */
+  ierr     = MatShellGetContext(Jac, &jctx);CHKERRQ(ierr);
+  ierr     = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
+  ierr     = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
+  ierr     = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr);
+  ierr     = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
+  numCells = cEnd - cStart;
+  for (field = 0; field < numFields; ++field) {
+    cellDof += quad[field].numBasisFuncs*quad[field].numComponents;
+  }
+  ierr = VecSet(F, 0.0);CHKERRQ(ierr);
+  ierr = PetscMalloc7(numCells*cellDof,PetscScalar,&u,numCells*cellDof,PetscScalar,&a,numCells*dim,PetscReal,&v0,numCells*dim*dim,PetscReal,&J,numCells*dim*dim,PetscReal,&invJ,numCells,PetscReal,&detJ,numCells*cellDof,PetscScalar,&elemVec);CHKERRQ(ierr);
+  for (c = cStart; c < cEnd; ++c) {
+    PetscScalar *x;
+    PetscInt     i;
+
+    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);
+    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);
+    for (i = 0; i < cellDof; ++i) a[c*cellDof+i] = x[i];
+    ierr = DMPlexVecRestoreClosure(dm, NULL, X, c, NULL, &x);CHKERRQ(ierr);
+  }
+  for (field = 0; field < numFields; ++field) {
+    const PetscInt numQuadPoints = quad[field].numQuadPoints;
+    const PetscInt numBasisFuncs = quad[field].numBasisFuncs;
+    /* Conforming batches */
+    PetscInt blockSize  = numBasisFuncs*numQuadPoints;
+    PetscInt numBlocks  = 1;
+    PetscInt batchSize  = numBlocks * blockSize;
+    PetscInt numBatches = numBatchesTmp;
+    PetscInt numChunks  = numCells / (numBatches*batchSize);
+    /* Remainder */
+    PetscInt numRemainder = numCells % (numBatches * batchSize);
+    PetscInt offset       = numCells - numRemainder;
+
+    ierr = (*mesh->integrateJacobianActionFEM)(numChunks*numBatches*batchSize, numFields, field, quad, u, a, v0, J, invJ, detJ, fem->g0Funcs, fem->g1Funcs, fem->g2Funcs, fem->g3Funcs, elemVec);CHKERRQ(ierr);
+    ierr = (*mesh->integrateJacobianActionFEM)(numRemainder, numFields, field, quad, &u[offset*cellDof], &a[offset*cellDof], &v0[offset*dim], &J[offset*dim*dim], &invJ[offset*dim*dim], &detJ[offset],
+                                               fem->g0Funcs, fem->g1Funcs, fem->g2Funcs, fem->g3Funcs, &elemVec[offset*cellDof]);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 = PetscFree7(u,a,v0,J,invJ,detJ,elemVec);CHKERRQ(ierr);
+  if (mesh->printFEM) {
+    PetscMPIInt rank, numProcs;
+    PetscInt    p;
+
+    ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr);
+    ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm), &numProcs);CHKERRQ(ierr);
+    ierr = PetscPrintf(PETSC_COMM_WORLD, "Jacobian Action:\n");CHKERRQ(ierr);
+    for (p = 0; p < numProcs; ++p) {
+      if (p == rank) {ierr = VecView(F, PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);}
+      ierr = PetscBarrier((PetscObject) dm);CHKERRQ(ierr);
+    }
+  }
+  /* ierr = PetscLogEventEnd(JacobianActionFEMEvent,0,0,0,0);CHKERRQ(ierr); */
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "DMPlexComputeJacobianFEM"
+/*@
+  DMPlexComputeJacobianFEM - Form the local portion of the Jacobian matrix J at the local solution X using pointwise functions specified by the user.
+
+  Input Parameters:
++ dm - The mesh
+. X  - Local input vector
+- user - The user context
+
+  Output Parameter:
+. Jac  - Jacobian matrix
+
+  Note:
+  The second member of the user context must be an FEMContext.
+
+  We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
+  like a GPU, or vectorize on a multicore machine.
+
+.seealso: FormFunctionLocal()