Commits

Jed Brown committed 1655452 Merge with conflicts

Merge branch 'knepley/plex'

Improved support for discretizing PDEs. Cell geometry and FVM. Little to
no topology.

* knepley/plex:
SNES ex12: Fix leak
TS ex11: error if dim != DIM because memory corruption occurs otherwise
DMPlex: All quadrilateral geometry tests working
DMPlex: remove unused variable
DMPlex: Hexahedron geometry tests now pass
DMPlex: More geometry work
DMPlex: use PETSC_HAVE_TGAMMA to get the build working on windows
configure: add test for tgamma() [this is missing on windows]
DMPlex: Added random transform tests for tets - Fixed memory leak
DMPlex: Fixed 3D geometry
DMPlex: FVM geometry tests now work for triangles
DMPlex: one more win32/complex fix
DMPlex: Enhanced geometry tests - Added random transforms of the reference triangle in 2D
DMPlex: Fixed complex problems in geometry calculation
DMPlex: Fixed Fortran test output
DMDT: M_PI --> PETSC_PI
DMPlex: Reorganized geometry testing
TS ex11: Put in 3D inflow and outflow ids from the command line
DMPlex: Added DMPlexRestoreCone/Support() to F90 interface definition header - Also added to test ex1f90
TS ex11: Need to make these ids input from the command line
DMPlex: Add Fortran bindings
DMPlex: Reorder tet output since I order tets differently
Builder: Link Fortran examples with Fortran linker
DMPlex: Fixed leak in example
DMPlex: Tutorial now shows how to output VTK
PetscViewer: Added VTK type for Fortran
DMPlex: Remove parameter names from prototype
DMPlex: Geometry fixes for C++ complex
DMAKKT: Remove from build
DMPlex: Now DMPlexCreateSection() only has an F90 interface
SNES: add Fortran bindings to DMSNESSet{Function,Jacobian}
SNES: remove superfluous casting in DM local Fortran bindings
SNES: compile zdmlocalsnesf.c
DMDT: Fix Fortran wrapper
DMPlex: Added Fortran interface for DMPlexMatSetClosure()
PetscSection: Added docs, enabled Fortran wrappers
DM: Added Fortran bindings for DMSNESSetFunction/JacobianLocal()
SNES ex12: Looks like 3D Neumann conditions are working
DMPlex: Fix volume determination for 2D in 3D
DMPlex: Better error reporting
SNES ex12: Boundary integration seems to work in 2D
DMPlex+FEM: Fix for boundary integration
TS ex11: Better error reporting
DMPlex ex8: Added FVM geometry tests
DMPlex: Fixed FVM geometry for 2D in 3D
SNES ex62: Now using PetscDT quadrature
DMDT: Translated FIAT's Gauss-Jacobi quadrature
DMPlex: Fix geometry tests - Fix tests for 2D projection
DMPlex: Geometry fixes
DMPlex: Small fixes for geometry
DMPlex: Fixed damn prototype
PetscSection: Added VecSetValuesSectionF90() - Fixed bug in Vec Fortran header
DMPlex: Added Fortran defines for additional insert modes
DMPlex: Silence warning in DMPlexInvertCells_Internal()
DMPlex: Functions passed to evaluation routines now return void and pass results in arguments
DMPlex: Now we give an explicit embedding dimension to the FVM geometry methods
DMPlex: Fixed complication with complex - Started to add 3D geometry stuff
TS ex11: Replace edge geometry with call to DMPlexComputeCellGeometryFVM()
DMPlex: Added normal argument to DMPlexComputeCellGeometryFVM() - This is only calculated for faces - Added calculation for faces in 2D
DMPlex: Damn makefiles
DMPlex: Added declaration for DMPlexCreateCGNS()
TS ex11: Reorganized computation of geometry in preparation for 3D
Builder: Now individual tests can have requirements - The key is 'requires', and it takes a list of package names - Also fixed up showSingleRun()
TS ex11: Use new DMPlexComputeCellGeometryFVM()
DMPlex: Added DMPlexComputeCellGeometryFVM() - Added internal volume methods
DMPlex: fix const and int vs. PetscInt
DMPlex: Better error reporting
DMPlex+ExodusII: Fixed reading of quads - Was broken by 2e1b13c25062c3c40593ce7412c5cd227259ade7
DMPlex: Hex cell geometry was broken - I don't see how this test passed before
DMPlex: Regression cleanup
DMPlex: Regression cleanup
SNES ex12: Fixed sign for Neumann BC
DMPlex+FEM: The dimension does not change for boundary elements
DMPlex: Geometry now works for boundary elements
DMPlex+FEM: Fixed quadrature coordinate handling for boundary integrals
SNES ex12: Reorganized tests, still working on Neumann conditions - Fixed boundary face label
DMPlex: Turn off CGNS by default
DMPlex: Added 1D cell geometry
DMPlex: DMPlexCreateCGNS()
DMPlex: Updated test ex1 output
DMPlex: Fixed compiler problems with complex
SNES ex12: Now does boundary integration - Still not verifying exact solution for Neumann conditions
DMPlex: Fix DMPlexComputeCellGeometry() to handle lower dimensional cells
DMPlex: Fixed error in projection from 3D to 2D - Need to check for case where normal is already z
DMPlex: Adding boundary integration
DMPlex: Fixed cell inversion for TetGen to commute with mesh interpolation
SNES ex12: Cleaning up ex12 testing
DMPlex: Added integration over boundary to DMPlexComputeResidualFEM() - Added quadBd, f0/1BdFuncs to PetscFEM struct - Needs to be tested
FEM Generation: We now allow a *_bd.h header that holds boundary discreization
Builder: Fixed regression requirements
Builder: Fixed specification of multiple tests to run
DMPlex ex7: Orientations now working properly
DMPlex: Changed orientation convention and fixed many bugs
SNEX ex12: Ignore generated headers
SNEX ex12 and ex52: Remove CTetGen warnings from test output - Fixed this in CTetGen repository
DMPlex: Fix default |J| for DMPlexComputeCellGeometry()
DMPlex: Fix for cell geometry (unbelievable screwup)
DMPlex ex7: Check that |J| is nonzero for each cell in interpolated mesh
DMPlex: Allow DMPlexComputeCellGeometry() to work with interpolated meshes
DMPlex: Added prototype for DMPlexCopyCoordinates(), and docs for that and DMPlexInterpolate()
DMPlex: Added F90 interface for DMPlexComputeCellGeometry()

Conflicts:
include/petscdmplex.h

  • Participants
  • Parent commits b30c3ef, 7436b39

Comments (0)

Files changed (94)

 src/dm/impls/complex/examples/tests/ex3_gpu_inline\.h
 src/dm/impls/complex/examples/tests/ex3_gpu\.h
 src/dm/impls/complex/examples/tests/ex3\.h
+src/snes/examples/tutorials/ex12.h
+src/snes/examples/tutorials/ex12_gpu.h
+src/snes/examples/tutorials/ex12_gpu_inline.h

File bin/pythonscripts/PetscGenerateFEMQuadrature.py

 generator  = PETSc.FEM.QuadratureGenerator()
 generator.setup()
 elements   = []
+bdElements = []
 if not (len(sys.argv)-2) % 5 == 0:
   sys.exit('Incomplete set of arguments')
 for n in range((len(sys.argv)-2) / 5):
   components = int(sys.argv[n*5+3])
   numBlocks  = int(sys.argv[n*5+4])
   operator   = sys.argv[n*5+5]
-  if order == 0:
-    element  = DiscontinuousLagrange(default_simplex(dim), order)
+  if operator == 'boundary':
+    if order == 0:
+      element  = DiscontinuousLagrange(default_simplex(dim-1), order)
+    else:
+      element  = Lagrange(default_simplex(dim-1), order)
+    element.numComponents = components
+    bdElements.append(element)
   else:
-    element  = Lagrange(default_simplex(dim), order)
-  element.numComponents = components
-  elements.append(element)
+    if order == 0:
+      element  = DiscontinuousLagrange(default_simplex(dim), order)
+    else:
+      element  = Lagrange(default_simplex(dim), order)
+    element.numComponents = components
+    elements.append(element)
 filename = sys.argv[-1]
 generator.quadDegree = max([e.order for e in elements])
-generator.run(elements, numBlocks, operator, filename)
+generator.run(elements, bdElements, numBlocks, operator, filename)

File bin/pythonscripts/PetscGenerateFEMQuadratureTensorProduct.py

   elements.append(element)
 filename = sys.argv[-1]
 generator.quadDegree = max([e.order+1 for e in elements])
-generator.runTensorProduct(dim, elements, numBlocks, operator, filename)
+generator.runTensorProduct(dim, elements, None, numBlocks, operator, filename)

File config/BuildSystem/config/libraries.py

       self.logPrint('Warning: erf() not found')
     return
 
+  def checkMathTgamma(self):
+    '''Check for tgama() in libm, the math library'''
+    if not self.math is None and self.check(self.math, ['tgamma'], prototype = ['double tgamma(double);'], call = ['double x = 0,y; y = tgamma(x);\n']):
+      self.logPrint('tgamma() found')
+      self.addDefine('HAVE_TGAMMA', 1)
+    else:
+      self.logPrint('Warning: tgamma() not found')
+    return
+
   def checkCompression(self):
     '''Check for libz, the compression library'''
     self.compression = None
     map(lambda args: self.executeTest(self.check, list(args)), self.libraries)
     self.executeTest(self.checkMath)
     self.executeTest(self.checkMathErf)
+    self.executeTest(self.checkMathTgamma)
     self.executeTest(self.checkCompression)
     self.executeTest(self.checkRealtime)
     self.executeTest(self.checkDynamic)

File config/PETSc/FEM.py

   /* ierr = PetscLogEventEnd(IntegrateJacobianEvent,0,0,0,0);CHKERRQ(ierr); */
   PetscFunctionReturn(0);
 }
+
+#undef __FUNCT__
+#define __FUNCT__ "FEMIntegrateBdResidualBatch"
+/*C
+  FEMIntegrateBdResidualBatch - Produce the element residual vector for a batch of elements by quadrature integration
+
+  Not collective
+
+  Input Parameters:
++ Ne                   - The number of elements in the batch
+. numFields            - The number of physical fields
+. field                - The field being integrated
+. quad                 - PetscQuadrature objects for each field
+. coefficients         - The array of FEM basis coefficients for the elements
+. v0s                  - The coordinates of the initial vertex for each element (the constant part of the transform from the reference element)
+. jacobians            - The Jacobian for each element (the linear part of the transform from the reference element)
+. jacobianInverses     - The Jacobian inverse for each element (the linear part of the transform to the reference element)
+. jacobianDeterminants - The Jacobian determinant for each element
+. f0_func              - f_0 function from the first order FEM model
+- f1_func              - f_1 function from the first order FEM model
+
+  Output Parameter
+. elemVec              - the element residual vectors from each element
+
+   Calling sequence of f0_func and f1_func:
+$    void f0(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], const PetscReal n[], PetscScalar f0[])
+
+  Note:
+$ Loop over batch of elements (e):
+$   Loop over quadrature points (q):
+$     Make u_q and gradU_q (loops over fields,Nb,Ncomp) and x_q
+$     Call f_0 and f_1
+$   Loop over element vector entries (f,fc --> i):
+$     elemVec[i] += \psi^{fc}_f(q) f0_{fc}(u, \nabla u) + \nabla\psi^{fc}_f(q) \cdot f1_{fc,df}(u, \nabla u)
+*/
+PetscErrorCode FEMIntegrateBdResidualBatch(PetscInt Ne, PetscInt numFields, PetscInt field, PetscQuadrature quad[], const PetscScalar coefficients[],
+                                           const PetscReal v0s[], const PetscReal normals[], const PetscReal jacobians[], const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[],
+                                           void (*f0_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], const PetscReal n[], PetscScalar f0[]),
+                                           void (*f1_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], const PetscReal n[], PetscScalar f1[]), PetscScalar elemVec[])
+{
+  const PetscInt debug   = 0;
+  const PetscInt dim     = SPATIAL_DIM_0;
+  const PetscInt numComponents = NUM_BASIS_COMPONENTS_TOTAL;
+  PetscInt       cOffset = 0;
+  PetscInt       eOffset = 0, e;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  /* ierr = PetscLogEventBegin(IntegrateResidualEvent,0,0,0,0);CHKERRQ(ierr); */
+  for (e = 0; e < Ne; ++e) {
+    const PetscReal  detJ = jacobianDeterminants[e];
+    const PetscReal *v0   = &v0s[e*dim];
+    const PetscReal *n    = &normals[e*dim];
+    const PetscReal *J    = &jacobians[e*dim*dim];
+    const PetscReal *invJ = &jacobianInverses[e*dim*dim];
+    const PetscInt   Nq   = quad[field].numQuadPoints;
+    PetscScalar      f0[NUM_QUADRATURE_POINTS_0*dim];
+    PetscScalar      f1[NUM_QUADRATURE_POINTS_0*dim*dim];
+    PetscInt         q, f;
+
+    if (Nq > NUM_QUADRATURE_POINTS_0) SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_LIB, "Number of quadrature points %d should be <= %d", Nq, NUM_QUADRATURE_POINTS_0);
+    if (debug > 1) {
+      ierr = PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", detJ);CHKERRQ(ierr);
+      ierr = DMPrintCellMatrix(e, "invJ", dim, dim, invJ);CHKERRQ(ierr);
+    }
+    for (q = 0; q < Nq; ++q) {
+      if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);CHKERRQ(ierr);}
+      PetscScalar      u[NUM_BASIS_COMPONENTS_TOTAL];
+      PetscScalar      gradU[dim*(NUM_BASIS_COMPONENTS_TOTAL)];
+      PetscReal        x[SPATIAL_DIM_0];
+      PetscInt         fOffset     = 0;
+      PetscInt         dOffset     = cOffset;
+      const PetscInt   Ncomp       = quad[field].numComponents;
+      const PetscReal *quadPoints  = quad[field].quadPoints;
+      const PetscReal *quadWeights = quad[field].quadWeights;
+      PetscInt         d, d2, f, i;
+
+      for (d = 0; d < numComponents; ++d)       {u[d]     = 0.0;}
+      for (d = 0; d < dim*(numComponents); ++d) {gradU[d] = 0.0;}
+      for (d = 0; d < dim; ++d) {
+        x[d] = v0[d];
+        for (d2 = 0; d2 < dim-1; ++d2) {
+          x[d] += J[d*dim+d2]*(quadPoints[q*(dim-1)+d2] + 1.0);
+        }
+      }
+      for (f = 0; f < numFields; ++f) {
+        const PetscInt   Nb       = quad[f].numBasisFuncs;
+        const PetscInt   Ncomp    = quad[f].numComponents;
+        const PetscReal *basis    = quad[f].basis;
+        const PetscReal *basisDer = quad[f].basisDer;
+        PetscInt         b, comp;
+
+        for (b = 0; b < Nb; ++b) {
+          for (comp = 0; comp < Ncomp; ++comp) {
+            const PetscInt cidx = b*Ncomp+comp;
+            PetscScalar    realSpaceDer[dim];
+            PetscInt       d, g;
+
+            u[fOffset+comp] += coefficients[dOffset+cidx]*basis[q*Nb*Ncomp+cidx];
+            for (d = 0; d < dim; ++d) {
+              realSpaceDer[d] = 0.0;
+              for (g = 0; g < dim-1; ++g) {
+                realSpaceDer[d] += invJ[g*dim+d]*basisDer[(q*Nb*Ncomp+cidx)*dim+g];
+              }
+              gradU[(fOffset+comp)*dim+d] += coefficients[dOffset+cidx]*realSpaceDer[d];
+            }
+          }
+        }
+        if (debug > 1) {
+          PetscInt d;
+          for (comp = 0; comp < Ncomp; ++comp) {
+            ierr = PetscPrintf(PETSC_COMM_SELF, "    u[%d,%d]: %g\n", f, comp, u[fOffset+comp]);CHKERRQ(ierr);
+            for (d = 0; d < dim; ++d) {
+              ierr = PetscPrintf(PETSC_COMM_SELF, "    gradU[%d,%d]_%c: %g\n", f, comp, 'x'+d, gradU[(fOffset+comp)*dim+d]);CHKERRQ(ierr);
+            }
+          }
+        }
+        fOffset += Ncomp;
+        dOffset += Nb*Ncomp;
+      }
+
+      f0_func(u, gradU, x, n, &f0[q*Ncomp]);
+      for (i = 0; i < Ncomp; ++i) {
+        f0[q*Ncomp+i] *= detJ*quadWeights[q];
+      }
+      f1_func(u, gradU, x, n, &f1[q*Ncomp*dim]);
+      for (i = 0; i < Ncomp*dim; ++i) {
+        f1[q*Ncomp*dim+i] *= detJ*quadWeights[q];
+      }
+      if (debug > 1) {
+        PetscInt c,d;
+        for (c = 0; c < Ncomp; ++c) {
+          ierr = PetscPrintf(PETSC_COMM_SELF, "    f0[%d]: %g\n", c, f0[q*Ncomp+c]);CHKERRQ(ierr);
+          for (d = 0; d < dim; ++d) {
+            ierr = PetscPrintf(PETSC_COMM_SELF, "    f1[%d]_%c: %g\n", c, 'x'+d, f1[(q*Ncomp + c)*dim+d]);CHKERRQ(ierr);
+          }
+        }
+      }
+      if (q == Nq-1) {cOffset = dOffset;}
+    }
+    for (f = 0; f < numFields; ++f) {
+      const PetscInt   Nq       = quad[f].numQuadPoints;
+      const PetscInt   Nb       = quad[f].numBasisFuncs;
+      const PetscInt   Ncomp    = quad[f].numComponents;
+      const PetscReal *basis    = quad[f].basis;
+      const PetscReal *basisDer = quad[f].basisDer;
+      PetscInt         b, comp;
+
+      if (f == field) {
+      for (b = 0; b < Nb; ++b) {
+        for (comp = 0; comp < Ncomp; ++comp) {
+          const PetscInt cidx = b*Ncomp+comp;
+          PetscInt       q;
+
+          elemVec[eOffset+cidx] = 0.0;
+          for (q = 0; q < Nq; ++q) {
+            PetscScalar realSpaceDer[dim];
+            PetscInt    d, g;
+
+            elemVec[eOffset+cidx] += basis[q*Nb*Ncomp+cidx]*f0[q*Ncomp+comp];
+            for (d = 0; d < dim; ++d) {
+              realSpaceDer[d] = 0.0;
+              for (g = 0; g < dim-1; ++g) {
+                realSpaceDer[d] += invJ[g*dim+d]*basisDer[(q*Nb*Ncomp+cidx)*dim+g];
+              }
+              elemVec[eOffset+cidx] += realSpaceDer[d]*f1[(q*Ncomp+comp)*dim+d];
+            }
+          }
+        }
+      }
+      if (debug > 1) {
+        PetscInt b, comp;
+
+        for (b = 0; b < Nb; ++b) {
+          for (comp = 0; comp < Ncomp; ++comp) {
+            ierr = PetscPrintf(PETSC_COMM_SELF, "    elemVec[%d,%d]: %g\n", b, comp, elemVec[eOffset+b*Ncomp+comp]);CHKERRQ(ierr);
+          }
+        }
+      }
+      }
+      eOffset += Nb*Ncomp;
+    }
+  }
+  /* TODO ierr = PetscLogFlops();CHKERRQ(ierr); */
+  /* ierr = PetscLogEventEnd(IntegrateResidualEvent,0,0,0,0);CHKERRQ(ierr); */
+  PetscFunctionReturn(0);
+}
 '''
 
 class QuadratureGenerator(script.Script):
     header.purpose    = CodePurpose.SKELETON
     return header
 
-  def getElementSource(self, elements, numBlocks = 1, operator = None, sourceType = 'CPU', tensor = 0):
+  def getElementSource(self, elements, numBlocks = 1, operator = None, sourceType = 'CPU', tensor = 0, isBoundary = False):
     from GenericCompiler import CompilerException
 
     self.logPrint('Generating element module', debugSection = 'codegen')
           defns.extend(self.getPhysicsRoutines(operator))
           defns.extend(self.getComputationLayoutStructs(numBlocks))
         elif sourceType == 'CPU':
-          defns.extend(self.getQuadratureStructs(2*len(quadrature.pts)-1, quadrature, n, tensor))
-          defns.extend(self.getBasisStructs(name, element, quadrature, n, tensor))
+          if isBoundary:
+            defns.extend(self.getQuadratureStructs(2*len(quadrature.pts)-1, quadrature, str(n)+'_BD', tensor))
+            defns.extend(self.getBasisStructs(name, element, quadrature, str(n)+'_BD', tensor))
+          else:
+            defns.extend(self.getQuadratureStructs(2*len(quadrature.pts)-1, quadrature, n, tensor))
+            defns.extend(self.getBasisStructs(name, element, quadrature, n, tensor))
           #defns.extend(self.getIntegratorPoints(n, element))
           #defns.extend(self.getIntegratorSetup(n, element))
           #defns.extend(self.getIntegratorSetup(n, element, True))
           n += element.value_shape()[0]
         else:
           n += 1
-      defns.extend(self.getAggregateBasisStructs(elements))
+      if not isBoundary:
+        defns.extend(self.getAggregateBasisStructs(elements))
       #defns.extend(self.getQuadratureSetup())
       #defns.extend(self.getElementIntegrals())
     except CompilerException, e:
         print e
     return
 
-  def run(self, elements, numBlocks, operator, filename = ''):
+  def run(self, elements, bdElements, numBlocks, operator, filename = ''):
     import os
     if elements is None:
       from FIAT.reference_element import default_simplex
     self.outputElementSource(self.getElementSource(elements), filename, extra = femIntegrationCode)
     self.outputElementSource(self.getElementSource(elements, numBlocks, operator, sourceType = 'GPU'), os.path.splitext(filename)[0]+'_gpu'+os.path.splitext(filename)[1])
     self.outputElementSource(self.getElementSource(elements, numBlocks, operator, sourceType = 'GPU_inline'), os.path.splitext(filename)[0]+'_gpu_inline'+os.path.splitext(filename)[1])
+    if bdElements:
+      self.outputElementSource(self.getElementSource(bdElements, isBoundary = True), os.path.splitext(filename)[0]+'_bd'+os.path.splitext(filename)[1])
     return
 
   def runTensorProduct(self, dim, elements, numBlocks, operator, filename = ''):

File config/builder.py

 
 import logger, script
 
-regressionRequirements = {'src/vec/vec/examples/tests/ex31':  set(['Matlab']),
-                          'src/vec/vec/examples/tests/ex52':  set(['Cuda'])
-                          }
-
-regressionParameters = {'src/sys/comm/examples/tests/ex1':    [{'numProcs': 2},
-                                                               {'numProcs': 5}],
-                        'src/vec/vec/examples/tests/ex1_2':    {'numProcs': 2},
-                        'src/vec/vec/examples/tests/ex3':      {'numProcs': 2},
-                        'src/vec/vec/examples/tests/ex4':      {'numProcs': 2},
-                        'src/vec/vec/examples/tests/ex5':      {'numProcs': 2},
-                        'src/vec/vec/examples/tests/ex9':      {'numProcs': 2},
-                        'src/vec/vec/examples/tests/ex10':     {'numProcs': 2},
-                        'src/vec/vec/examples/tests/ex11':     {'numProcs': 2},
-                        'src/vec/vec/examples/tests/ex12':     {'numProcs': 2},
-                        'src/vec/vec/examples/tests/ex13':     {'numProcs': 2},
-                        'src/vec/vec/examples/tests/ex14':     {'numProcs': 2},
-                        'src/vec/vec/examples/tests/ex16':     {'numProcs': 2},
-                        'src/vec/vec/examples/tests/ex17':     {'numProcs': 2},
-                        'src/vec/vec/examples/tests/ex17f':    {'numProcs': 3},
-                        'src/vec/vec/examples/tests/ex21_2':   {'numProcs': 2},
-                        'src/vec/vec/examples/tests/ex22':     {'numProcs': 4},
-                        'src/vec/vec/examples/tests/ex23':     {'numProcs': 2},
-                        'src/vec/vec/examples/tests/ex24':     {'numProcs': 3},
-                        'src/vec/vec/examples/tests/ex25':     {'numProcs': 3},
-                        'src/vec/vec/examples/tests/ex26':     {'numProcs': 4},
-                        'src/vec/vec/examples/tests/ex28':     {'numProcs': 3},
-                        'src/vec/vec/examples/tests/ex29':     {'numProcs': 3, 'args': '-n 126'},
-                        'src/vec/vec/examples/tests/ex30f':    {'numProcs': 4},
-                        'src/vec/vec/examples/tests/ex33':     {'numProcs': 4},
-                        'src/vec/vec/examples/tests/ex36':     {'numProcs': 2, 'args': '-set_option_negidx -set_values_negidx -get_values_negidx'},
-                        'src/vec/vec/examples/tutorials/ex10':  [{'numProcs': 1, 'args': '-hdf5'},
-                                                                 {'numProcs': 2, 'args': '-binary'},
-                                                                 {'numProcs': 3, 'args': '-binary'},
-                                                                 {'numProcs': 5, 'args': '-binary'}],
-                        'src/dm/impls/patch/examples/tests/ex1': [{'numProcs': 1, 'args': '-patch_size 2 -grid_size 6'},
+regressionParameters = {'src/dm/impls/patch/examples/tests/ex1': [{'numProcs': 1, 'args': '-patch_size 2 -grid_size 6'},
                                                                   {'numProcs': 4, 'args': '-patch_size 2 -grid_size 4'},
                                                                   {'numProcs': 4, 'args': '-patch_size 2 -grid_size 4 -comm_size 1'},
                                                                   {'numProcs': 4, 'args': '-patch_size 2 -grid_size 4 -comm_size 2'},
                                                                   {'numProcs': 16, 'args': '-patch_size 2 -grid_size 8 -comm_size 2'}],
-                        'src/dm/impls/plex/examples/tests/ex1': [{'numProcs': 1, 'args': '-dim 3 -ctetgen_verbose 4 -dm_view ::ascii_info_detail -info -info_exclude null'},
+                        'src/dm/impls/plex/examples/tests/ex1': [# CTetGen tests 0-1
+                                                                 {'numProcs': 1, 'args': '-dim 3 -ctetgen_verbose 4 -dm_view ::ascii_info_detail -info -info_exclude null'},
                                                                  {'numProcs': 1, 'args': '-dim 3 -ctetgen_verbose 4 -refinement_limit 0.0625 -dm_view ::ascii_info_detail -info -info_exclude null'},
+                                                                 # Test 2D LaTex and ASCII output 2-9
                                                                  {'numProcs': 1, 'args': '-dim 2 -dm_view ::ascii_latex'},
                                                                  {'numProcs': 1, 'args': '-dim 2 -refinement_uniform 1 -interpolate 1 -dm_view ::ascii_info_detail'},
                                                                  {'numProcs': 2, 'args': '-dim 2 -refinement_uniform 1 -interpolate 1 -dm_view ::ascii_info_detail'},
                                                                  {'numProcs': 1, 'args': '-dim 2 -cell_simplex 0 -dm_view ::ascii_info_detail'},
                                                                  {'numProcs': 1, 'args': '-dim 2 -cell_simplex 0 -refinement_uniform 1 -dm_view ::ascii_info_detail'},
                                                                  {'numProcs': 2, 'args': '-dim 2 -cell_simplex 0 -refinement_uniform 1 -interpolate 1 -dm_view ::ascii_info_detail'},
-                                                                 {'numProcs': 2, 'args': '-dim 2 -cell_simplex 0 -refinement_uniform 1 -interpolate 1 -dm_view ::ascii_latex'}],
+                                                                 {'numProcs': 2, 'args': '-dim 2 -cell_simplex 0 -refinement_uniform 1 -interpolate 1 -dm_view ::ascii_latex'},
+                                                                 # CGNS tests 10-11 (need to find smaller test meshes)
+                                                                 {'numProcs': 1, 'args': '-filename %(meshes)s/tut21.cgns -interpolate 1 -dm_view', 'requires': ['CGNS']},
+                                                                 {'numProcs': 1, 'args': '-filename %(meshes)s/StaticMixer.cgns -interpolate 1 -dm_view', 'requires': ['CGNS']}],
                         'src/dm/impls/plex/examples/tests/ex3': [{'numProcs': 1, 'args': '',
                                                                   'setup': './bin/pythonscripts/PetscGenerateFEMQuadrature.py 2 1 2 1 identity src/dm/impls/plex/examples/tests/ex3.h'},
                                                                  {'numProcs': 1, 'args': '-order 1'},
                                                                  {'numProcs': 1, 'args': '-malloc_dump -pend 10000'},
                                                                  {'numProcs': 1, 'args': '-malloc_dump -pend 10000 -fill 0.05'},
                                                                  {'numProcs': 1, 'args': '-malloc_dump -pend 10000 -fill 0.25'}],
-                        'src/dm/impls/plex/examples/tests/ex7': [{'numProcs': 1, 'args': '-dim 2 -dm_view ::ascii_info_detail'},
+                        'src/dm/impls/plex/examples/tests/ex7': [# Two cell test meshes 0-7
+                                                                 {'numProcs': 1, 'args': '-dim 2 -dm_view ::ascii_info_detail'},
                                                                  {'numProcs': 2, 'args': '-dim 2 -dm_view ::ascii_info_detail'},
                                                                  {'numProcs': 1, 'args': '-dim 2 -cell_simplex 0 -dm_view ::ascii_info_detail'},
                                                                  {'numProcs': 2, 'args': '-dim 2 -cell_simplex 0 -dm_view ::ascii_info_detail'},
                                                                  {'numProcs': 2, 'args': '-dim 3 -dm_view ::ascii_info_detail'},
                                                                  {'numProcs': 1, 'args': '-dim 3 -cell_simplex 0 -dm_view ::ascii_info_detail'},
                                                                  {'numProcs': 2, 'args': '-dim 3 -cell_simplex 0 -dm_view ::ascii_info_detail'},
-                                                                 {'numProcs': 1, 'args': '-dim 2 -cell_simplex 0 -testnum 1 -dm_view ::ascii_info_detail'}],
-                        'src/dm/impls/plex/examples/tests/ex8': [{'numProcs': 1, 'args': '-dm_view ::ascii_info_detail'}],
+                                                                 # 2D Hybrid Mesh 8
+                                                                 {'numProcs': 1, 'args': '-dim 2 -cell_simplex 0 -testnum 1 -dm_view ::ascii_info_detail'},
+                                                                 # TetGen meshes 9-10
+                                                                 {'numProcs': 1, 'args': '-dim 2 -use_generator -dm_view ::ascii_info_detail'},
+                                                                 {'numProcs': 1, 'args': '-dim 3 -use_generator -dm_view ::ascii_info_detail'},
+                                                                 # Cubit meshes 11
+                                                                 {'numProcs': 1, 'args': '-dim 3 -filename %(meshes)s/blockcylinder-50.exo -dm_view ::ascii_info_detail'},
+                                                                 #{'numProcs': 1, 'args': '-dim 3 -filename /PETSc3/petsc/blockcylinder-50.exo -dm_view ::ascii_info_detail'},
+                                                                 #{'numProcs': 1, 'args': '-dim 3 -filename /PETSc3/petsc/blockcylinder-20.exo'},
+                                                                 ],
+                        'src/dm/impls/plex/examples/tests/ex8': [{'numProcs': 1, 'args': '-dm_view ::ascii_info_detail'},
+                                                                 {'numProcs': 1, 'args': '-transform'}],
                         'src/dm/impls/plex/examples/tests/ex1f90': [{'numProcs': 1, 'args': ''}],
                         'src/dm/impls/plex/examples/tests/ex2f90': [{'numProcs': 1, 'args': ''}],
                         'src/dm/impls/plex/examples/tutorials/ex1': [{'numProcs': 1, 'args': ''},
                                                                             {'numProcs': 1, 'args': '-i src/dm/impls/mesh/examples/tests/meshes/Hex.gen'},
                                                                             {'numProcs': 1, 'args': '-i src/dm/impls/mesh/examples/tests/meshes/Tet.gen -interpolate'},
                                                                             {'numProcs': 1, 'args': '-i src/dm/impls/mesh/examples/tests/meshes/Hex.gen -interpolate'}],
-                        'src/ksp/ksp/examples/tests/ex10':    [{'numProcs': 1, 'args': '-m 3 -matconvert_type seqaij -ksp_monitor_short'}],
-                        'src/ksp/ksp/examples/tutorials/ex4': [{'numProcs': 1, 'args': '-info -info_exclude null,vec'},
-                                                               {'numProcs': 1, 'args': '-da_grid_x 10 -da_grid_y 10 -solve -ksp_monitor'},
-                                                               {'numProcs': 2, 'args': '-da_grid_x 10 -da_grid_y 10'}],
-                                                               #{'numProcs': 1, 'args': '-info'}],
-                        'src/ksp/ksp/examples/tutorials/ex12': {'numProcs': 2, 'args': '-ksp_gmres_cgs_refinement_type refine_always'},
-                        'src/ksp/ksp/examples/tutorials/ex40': {'numProcs': 1, 'args': '-mat_no_inode -ksp_monitor_short'},
-                        'src/ksp/ksp/examples/tutorials/ex54':[{'numProcs': 4, 'args': '-ne 40 -alpha 1.e-3 -ksp_monitor_short -ksp_type cg -ksp_norm_type unpreconditioned'},
-                                                               {'numProcs': 4, 'args': '-ne 40 -alpha 1.e-3 -ksp_monitor_short -ksp_type cg -ksp_norm_type unpreconditioned -pc_gamg_type sa'}],
-                        'src/ksp/ksp/examples/tutorials/ex55':[{'numProcs': 4, 'args': '-ne 40 -alpha 1.e-3 -ksp_monitor_short -ksp_type cg -ksp_norm_type unpreconditioned'},
-                                                               {'numProcs': 4, 'args': '-ne 40 -alpha 1.e-3 -ksp_monitor_short -ksp_type cg -ksp_norm_type unpreconditioned -pc_gamg_type sa'}],
-                        'src/ksp/ksp/examples/tutorials/ex56':[{'numProcs': 8, 'args': '-ne 11 -alpha 1.e-3 -ksp_monitor_short -ksp_type cg -ksp_norm_type unpreconditioned -pc_gamg_type sa'}],
-                        'src/snes/examples/tutorials/ex12':   [# 2D serial P1 test
+                        'src/snes/examples/tutorials/ex12':   [# 2D serial P1 test 0-6
                                                                {'numProcs': 1, 'args': '-run_type test -refinement_limit 0.0    -bc_type dirichlet -interpolate 0 -show_initial -dm_plex_print_fem 1',
-                                                                'setup': './bin/pythonscripts/PetscGenerateFEMQuadrature.py 2 1 1 1 laplacian src/snes/examples/tutorials/ex12.h'},
-                                                               # 3D serial P1 test
+                                                                'setup': './bin/pythonscripts/PetscGenerateFEMQuadrature.py 2 1 1 1 laplacian 2 1 1 1 boundary src/snes/examples/tutorials/ex12.h'},
+                                                               {'numProcs': 1, 'args': '-run_type test -refinement_limit 0.0    -bc_type dirichlet -interpolate 1 -show_initial -dm_plex_print_fem 1'},
+                                                               {'numProcs': 1, 'args': '-run_type test -refinement_limit 0.0625 -bc_type dirichlet -interpolate 1 -show_initial -dm_plex_print_fem 1'},
+                                                               {'numProcs': 1, 'args': '-run_type test -refinement_limit 0.0    -bc_type neumann   -interpolate 1 -show_initial -dm_plex_print_fem 1 -dm_view ::ascii_info_detail'},
+                                                               {'numProcs': 1, 'args': '-run_type test -refinement_limit 0.0625 -bc_type neumann   -interpolate 1 -show_initial -dm_plex_print_fem 1'},
+                                                               {'numProcs': 1, 'args': '-run_type test -refinement_limit 0.0    -bc_type neumann   -interpolate 1 -show_initial -dm_plex_print_fem 1 -dm_view ::ascii_info_detail',
+                                                                'setup': './bin/pythonscripts/PetscGenerateFEMQuadrature.py 2 2 1 1 laplacian 2 2 1 1 boundary src/snes/examples/tutorials/ex12.h'},
+                                                               {'numProcs': 1, 'args': '-run_type test -refinement_limit 0.0625 -bc_type neumann   -interpolate 1 -show_initial -dm_plex_print_fem 1 -dm_view ::ascii_info_detail'},
+                                                               # 3D serial P1 test 7-
                                                                {'numProcs': 1, 'args': '-run_type test -dim 3 -refinement_limit 0.0    -bc_type dirichlet -interpolate 0 -show_initial -dm_plex_print_fem 1 -dm_view',
-                                                                'setup': './bin/pythonscripts/PetscGenerateFEMQuadrature.py 3 1 1 1 laplacian src/snes/examples/tutorials/ex12.h'},                                                                                    {'numProcs': 1, 'args': '-run_type test -dim 3 -refinement_limit 0.0125 -bc_type dirichlet -interpolate 0 -show_initial -dm_plex_print_fem 1 -dm_view'},
-                                                               {'numProcs': 1, 'args': '-run_type test -dim 3 -refinement_limit 0.0    -bc_type neumann   -interpolate 0 -snes_fd -show_initial -dm_plex_print_fem 1 -dm_view'},
-                                                               # Using ExodusII mesh  ::ascii_info_detail
+                                                                'setup': './bin/pythonscripts/PetscGenerateFEMQuadrature.py 3 1 1 1 laplacian 3 1 1 1 boundary src/snes/examples/tutorials/ex12.h'},                                                                   {'numProcs': 1, 'args': '-run_type test -dim 3 -refinement_limit 0.0    -bc_type dirichlet -interpolate 1 -show_initial -dm_plex_print_fem 1 -dm_view'},
+                                                               {'numProcs': 1, 'args': '-run_type test -dim 3 -refinement_limit 0.0125 -bc_type dirichlet -interpolate 1 -show_initial -dm_plex_print_fem 1 -dm_view'},
+                                                               {'numProcs': 1, 'args': '-run_type test -dim 3 -refinement_limit 0.0    -bc_type neumann   -interpolate 1 -snes_fd -show_initial -dm_plex_print_fem 1 -dm_view'},
+                                                               # Using ExodusII mesh 10-11 BROKEN
                                                                {'numProcs': 1, 'args': '-run_type test -f %(meshes)s/sevenside.exo -refinement_limit 0.0    -bc_type dirichlet -interpolate 1 -show_initial -dm_plex_print_fem 1 -dm_view',
                                                                 'setup': './bin/pythonscripts/PetscGenerateFEMQuadrature.py 2 1 1 1 laplacian src/snes/examples/tutorials/ex12.h'},
                                                                {'numProcs': 1, 'args': '-run_type test -dim 3 -f /Users/knepley/Downloads/kis_modell_tet2.exo -refinement_limit 0.0    -bc_type dirichlet -interpolate 1 -show_initial -dm_plex_print_fem 1 -dm_view',
                         'src/snes/examples/tutorials/ex52':   [# 2D Laplacian 0-3
                                                                {'numProcs': 1, 'args': '-dm_view -refinement_limit 0.0 -compute_function',
                                                                 'setup': './bin/pythonscripts/PetscGenerateFEMQuadrature.py 2 1 1 1 laplacian src/snes/examples/tutorials/ex52.h',
-                                                                'source': ['src/snes/examples/tutorials/ex52_integrateElement.cu']},
-                                                               {'numProcs': 1, 'args': '-dm_view -refinement_limit 0.0 -compute_function -batch'},
-                                                               {'numProcs': 1, 'args': '-dm_view -refinement_limit 0.0 -compute_function -batch -gpu'},
-                                                               {'numProcs': 1, 'args': '-dm_view -refinement_limit 0.0 -compute_function -batch -gpu -gpu_batches 2'},
+                                                                'source': ['src/snes/examples/tutorials/ex52_integrateElement.cu'], 'requires': ['Cuda']},
+                                                               {'numProcs': 1, 'args': '-dm_view -refinement_limit 0.0 -compute_function -batch', 'requires': ['Cuda']},
+                                                               {'numProcs': 1, 'args': '-dm_view -refinement_limit 0.0 -compute_function -batch -gpu', 'requires': ['Cuda']},
+                                                               {'numProcs': 1, 'args': '-dm_view -refinement_limit 0.0 -compute_function -batch -gpu -gpu_batches 2', 'requires': ['Cuda']},
                                                                # 2D Laplacian refined 4-8
-                                                               {'numProcs': 1, 'args': '-dm_view -refinement_limit 0.0625 -compute_function'},
-                                                               {'numProcs': 1, 'args': '-dm_view -refinement_limit 0.0625 -compute_function -batch'},
-                                                               {'numProcs': 1, 'args': '-dm_view -refinement_limit 0.0625 -compute_function -batch -gpu'},
-                                                               {'numProcs': 1, 'args': '-dm_view -refinement_limit 0.0625 -compute_function -batch -gpu -gpu_batches 2'},
-                                                               {'numProcs': 1, 'args': '-dm_view -refinement_limit 0.0625 -compute_function -batch -gpu -gpu_batches 4'},
+                                                               {'numProcs': 1, 'args': '-dm_view -refinement_limit 0.0625 -compute_function', 'requires': ['Cuda']},
+                                                               {'numProcs': 1, 'args': '-dm_view -refinement_limit 0.0625 -compute_function -batch', 'requires': ['Cuda']},
+                                                               {'numProcs': 1, 'args': '-dm_view -refinement_limit 0.0625 -compute_function -batch -gpu', 'requires': ['Cuda']},
+                                                               {'numProcs': 1, 'args': '-dm_view -refinement_limit 0.0625 -compute_function -batch -gpu -gpu_batches 2', 'requires': ['Cuda']},
+                                                               {'numProcs': 1, 'args': '-dm_view -refinement_limit 0.0625 -compute_function -batch -gpu -gpu_batches 4', 'requires': ['Cuda']},
                                                                # 2D Elasticity 9-12
                                                                {'numProcs': 1, 'args': '-dm_view -refinement_limit 0.0 -compute_function -op_type elasticity',
-                                                                'setup': './bin/pythonscripts/PetscGenerateFEMQuadrature.py 2 1 2 1 elasticity src/snes/examples/tutorials/ex52.h'},
-                                                               {'numProcs': 1, 'args': '-dm_view -refinement_limit 0.0 -compute_function -op_type elasticity -batch'},
-                                                               {'numProcs': 1, 'args': '-dm_view -refinement_limit 0.0 -compute_function -op_type elasticity -batch -gpu'},
-                                                               {'numProcs': 1, 'args': '-dm_view -refinement_limit 0.0 -compute_function -op_type elasticity -batch -gpu -gpu_batches 2'},
+                                                                'setup': './bin/pythonscripts/PetscGenerateFEMQuadrature.py 2 1 2 1 elasticity src/snes/examples/tutorials/ex52.h',
+                                                                'requires': ['Cuda']},
+                                                               {'numProcs': 1, 'args': '-dm_view -refinement_limit 0.0 -compute_function -op_type elasticity -batch', 'requires': ['Cuda']},
+                                                               {'numProcs': 1, 'args': '-dm_view -refinement_limit 0.0 -compute_function -op_type elasticity -batch -gpu', 'requires': ['Cuda']},
+                                                               {'numProcs': 1, 'args': '-dm_view -refinement_limit 0.0 -compute_function -op_type elasticity -batch -gpu -gpu_batches 2', 'requires': ['Cuda']},
                                                                # 2D Elasticity refined 13-17
-                                                               {'numProcs': 1, 'args': '-dm_view -refinement_limit 0.0625 -compute_function -op_type elasticity'},
-                                                               {'numProcs': 1, 'args': '-dm_view -refinement_limit 0.0625 -compute_function -op_type elasticity -batch'},
-                                                               {'numProcs': 1, 'args': '-dm_view -refinement_limit 0.0625 -compute_function -op_type elasticity -batch -gpu'},
-                                                               {'numProcs': 1, 'args': '-dm_view -refinement_limit 0.0625 -compute_function -op_type elasticity -batch -gpu -gpu_batches 2'},
-                                                               {'numProcs': 1, 'args': '-dm_view -refinement_limit 0.0625 -compute_function -op_type elasticity -batch -gpu -gpu_batches 4'},
+                                                               {'numProcs': 1, 'args': '-dm_view -refinement_limit 0.0625 -compute_function -op_type elasticity', 'requires': ['Cuda']},
+                                                               {'numProcs': 1, 'args': '-dm_view -refinement_limit 0.0625 -compute_function -op_type elasticity -batch', 'requires': ['Cuda']},
+                                                               {'numProcs': 1, 'args': '-dm_view -refinement_limit 0.0625 -compute_function -op_type elasticity -batch -gpu', 'requires': ['Cuda']},
+                                                               {'numProcs': 1, 'args': '-dm_view -refinement_limit 0.0625 -compute_function -op_type elasticity -batch -gpu -gpu_batches 2', 'requires': ['Cuda']},
+                                                               {'numProcs': 1, 'args': '-dm_view -refinement_limit 0.0625 -compute_function -op_type elasticity -batch -gpu -gpu_batches 4', 'requires': ['Cuda']},
                                                                # 3D Laplacian 18-20
                                                                {'numProcs': 1, 'args': '-dim 3 -dm_view -refinement_limit 0.0 -compute_function',
-                                                                'setup': './bin/pythonscripts/PetscGenerateFEMQuadrature.py 3 1 1 1 laplacian src/snes/examples/tutorials/ex52.h'},
-                                                               {'numProcs': 1, 'args': '-dim 3 -dm_view -refinement_limit 0.0 -compute_function -batch'},
-                                                               {'numProcs': 1, 'args': '-dim 3 -dm_view -refinement_limit 0.0 -compute_function -batch -gpu'},
+                                                                'setup': './bin/pythonscripts/PetscGenerateFEMQuadrature.py 3 1 1 1 laplacian src/snes/examples/tutorials/ex52.h',
+                                                                'requires': ['Cuda']},
+                                                               {'numProcs': 1, 'args': '-dim 3 -dm_view -refinement_limit 0.0 -compute_function -batch', 'requires': ['Cuda']},
+                                                               {'numProcs': 1, 'args': '-dim 3 -dm_view -refinement_limit 0.0 -compute_function -batch -gpu', 'requires': ['Cuda']},
                                                                # 3D Laplacian refined 21-24
-                                                               {'numProcs': 1, 'args': '-dim 3 -dm_view -refinement_limit 0.0125 -compute_function'},
-                                                               {'numProcs': 1, 'args': '-dim 3 -dm_view -refinement_limit 0.0125 -compute_function -batch'},
-                                                               {'numProcs': 1, 'args': '-dim 3 -dm_view -refinement_limit 0.0125 -compute_function -batch -gpu'},
-                                                               {'numProcs': 1, 'args': '-dim 3 -dm_view -refinement_limit 0.0125 -compute_function -batch -gpu -gpu_batches 2'},
+                                                               {'numProcs': 1, 'args': '-dim 3 -dm_view -refinement_limit 0.0125 -compute_function', 'requires': ['Cuda']},
+                                                               {'numProcs': 1, 'args': '-dim 3 -dm_view -refinement_limit 0.0125 -compute_function -batch', 'requires': ['Cuda']},
+                                                               {'numProcs': 1, 'args': '-dim 3 -dm_view -refinement_limit 0.0125 -compute_function -batch -gpu', 'requires': ['Cuda']},
+                                                               {'numProcs': 1, 'args': '-dim 3 -dm_view -refinement_limit 0.0125 -compute_function -batch -gpu -gpu_batches 2', 'requires': ['Cuda']},
                                                                # 3D Elasticity 25-27
                                                                {'numProcs': 1, 'args': '-dim 3 -dm_view -refinement_limit 0.0 -compute_function -op_type elasticity',
-                                                                'setup': './bin/pythonscripts/PetscGenerateFEMQuadrature.py 3 1 3 1 elasticity src/snes/examples/tutorials/ex52.h'},
-                                                               {'numProcs': 1, 'args': '-dim 3 -dm_view -refinement_limit 0.0 -compute_function -op_type elasticity -batch'},
-                                                               {'numProcs': 1, 'args': '-dim 3 -dm_view -refinement_limit 0.0 -compute_function -op_type elasticity -batch -gpu'},
+                                                                'setup': './bin/pythonscripts/PetscGenerateFEMQuadrature.py 3 1 3 1 elasticity src/snes/examples/tutorials/ex52.h',
+                                                                'requires': ['Cuda']},
+                                                               {'numProcs': 1, 'args': '-dim 3 -dm_view -refinement_limit 0.0 -compute_function -op_type elasticity -batch', 'requires': ['Cuda']},
+                                                               {'numProcs': 1, 'args': '-dim 3 -dm_view -refinement_limit 0.0 -compute_function -op_type elasticity -batch -gpu', 'requires': ['Cuda']},
                                                                # 3D Elasticity refined 28-31
-                                                               {'numProcs': 1, 'args': '-dim 3 -dm_view -refinement_limit 0.0125 -compute_function -op_type elasticity'},
-                                                               {'numProcs': 1, 'args': '-dim 3 -dm_view -refinement_limit 0.0125 -compute_function -op_type elasticity -batch'},
-                                                               {'numProcs': 1, 'args': '-dim 3 -dm_view -refinement_limit 0.0125 -compute_function -op_type elasticity -batch -gpu'},
-                                                               {'numProcs': 1, 'args': '-dim 3 -dm_view -refinement_limit 0.0125 -compute_function -op_type elasticity -batch -gpu -gpu_batches 2'},
+                                                               {'numProcs': 1, 'args': '-dim 3 -dm_view -refinement_limit 0.0125 -compute_function -op_type elasticity', 'requires': ['Cuda']},
+                                                               {'numProcs': 1, 'args': '-dim 3 -dm_view -refinement_limit 0.0125 -compute_function -op_type elasticity -batch', 'requires': ['Cuda']},
+                                                               {'numProcs': 1, 'args': '-dim 3 -dm_view -refinement_limit 0.0125 -compute_function -op_type elasticity -batch -gpu', 'requires': ['Cuda']},
+                                                               {'numProcs': 1, 'args': '-dim 3 -dm_view -refinement_limit 0.0125 -compute_function -op_type elasticity -batch -gpu -gpu_batches 2', 'requires': ['Cuda']},
                                                                # 'source': ['src/snes/examples/tutorials/ex52_integrateElement.cu']},
                                                                ],
                         'src/snes/examples/tutorials/ex62':   [# 2D serial P1 tests 0-3
                                                                {'numProcs': 1, 'args': '-run_type test -refinement_limit 0.0625 -bc_type dirichlet -interpolate 0 -show_initial -dm_plex_print_fem 1'},
                                                                {'numProcs': 1, 'args': '-run_type test -refinement_limit 0.0625 -bc_type dirichlet -interpolate 1 -show_initial -dm_plex_print_fem 1'},
                                                                # 2D serial P2 tests 4-5
-                                                               {'numProcs': 1, 'args': '-run_type test -refinement_limit 0.0    -bc_type dirichlet -interpolate 1 -show_initial -dm_plex_print_fem 1',
+                                                               {'numProcs': 1, 'args': '-run_type test -refinement_limit 0.0    -bc_type dirichlet -interpolate 1 -order 2,1 -show_initial -dm_plex_print_fem 1',
                                                                 'setup': './bin/pythonscripts/PetscGenerateFEMQuadrature.py 2 2 2 1 laplacian 2 1 1 1 gradient src/snes/examples/tutorials/ex62.h'},
 
-                                                               {'numProcs': 1, 'args': '-run_type test -refinement_limit 0.0625 -bc_type dirichlet -interpolate 1 -show_initial -dm_plex_print_fem 1'},
+                                                               {'numProcs': 1, 'args': '-run_type test -refinement_limit 0.0625 -bc_type dirichlet -interpolate 1 -order 2,1 -show_initial -dm_plex_print_fem 1'},
                                                                # Parallel tests 6-17
                                                                {'numProcs': 2, 'args': '-run_type test -refinement_limit 0.0    -bc_type dirichlet -interpolate 0 -dm_plex_print_fem 1',
                                                                 'setup': './bin/pythonscripts/PetscGenerateFEMQuadrature.py 2 1 2 1 laplacian 2 1 1 1 gradient src/snes/examples/tutorials/ex62.h'},
                                                                {'numProcs': 2, 'args': '-run_type full -refinement_limit 0.0625 -bc_type dirichlet -interpolate 1 -pc_type jacobi -ksp_rtol 1.0e-9 -snes_monitor_short -ksp_monitor_short -snes_converged_reason -snes_view'},
                                                                {'numProcs': 3, 'args': '-run_type full -refinement_limit 0.0625 -bc_type dirichlet -interpolate 1 -pc_type jacobi -ksp_rtol 1.0e-9 -snes_converged_reason -snes_view'},
                                                                {'numProcs': 5, 'args': '-run_type full -refinement_limit 0.0625 -bc_type dirichlet -interpolate 1 -pc_type jacobi -ksp_rtol 1.0e-9 -snes_converged_reason -snes_view'},
-                                                               {'numProcs': 1, 'args': '-run_type full -refinement_limit 0.0625 -bc_type dirichlet -interpolate 1 -pc_type jacobi -ksp_rtol 1.0e-9 -snes_monitor_short -ksp_monitor_short -snes_converged_reason -snes_view',
+                                                               {'numProcs': 1, 'args': '-run_type full -refinement_limit 0.0625 -bc_type dirichlet -interpolate 1 -order 2,1 -pc_type jacobi -ksp_rtol 1.0e-9 -snes_monitor_short -ksp_monitor_short -snes_converged_reason -snes_view',
                                                                 'setup': './bin/pythonscripts/PetscGenerateFEMQuadrature.py 2 2 2 1 laplacian 2 1 1 1 gradient src/snes/examples/tutorials/ex62.h'},
-                                                               {'numProcs': 2, 'args': '-run_type full -refinement_limit 0.0625 -bc_type dirichlet -interpolate 1 -pc_type jacobi -ksp_rtol 1.0e-9 -snes_monitor_short -ksp_monitor_short -snes_converged_reason -snes_view'},
-                                                               {'numProcs': 3, 'args': '-run_type full -refinement_limit 0.0625 -bc_type dirichlet -interpolate 1 -pc_type jacobi -ksp_rtol 1.0e-9 -snes_converged_reason -snes_view'},
-                                                               {'numProcs': 5, 'args': '-run_type full -refinement_limit 0.0625 -bc_type dirichlet -interpolate 1 -pc_type jacobi -ksp_rtol 1.0e-9 -snes_converged_reason -snes_view'},
+                                                               {'numProcs': 2, 'args': '-run_type full -refinement_limit 0.0625 -bc_type dirichlet -interpolate 1 -order 2,1 -pc_type jacobi -ksp_rtol 1.0e-9 -snes_monitor_short -ksp_monitor_short -snes_converged_reason -snes_view'},
+                                                               {'numProcs': 3, 'args': '-run_type full -refinement_limit 0.0625 -bc_type dirichlet -interpolate 1 -order 2,1 -pc_type jacobi -ksp_rtol 1.0e-9 -snes_converged_reason -snes_view'},
+                                                               {'numProcs': 5, 'args': '-run_type full -refinement_limit 0.0625 -bc_type dirichlet -interpolate 1 -order 2,1 -pc_type jacobi -ksp_rtol 1.0e-9 -snes_converged_reason -snes_view'},
                                                                # Stokes preconditioners 30-36
                                                                #   Jacobi
-                                                               {'numProcs': 1, 'args': '-run_type full -refinement_limit 0.00625 -bc_type dirichlet -interpolate 1 -ksp_gmres_restart 100 -pc_type jacobi -ksp_rtol 1.0e-9 -snes_monitor_short -ksp_monitor_short -snes_converged_reason -snes_view -show_solution 0',
+                                                               {'numProcs': 1, 'args': '-run_type full -refinement_limit 0.00625 -bc_type dirichlet -interpolate 1 -order 2,1 -ksp_gmres_restart 100 -pc_type jacobi -ksp_rtol 1.0e-9 -snes_monitor_short -ksp_monitor_short -snes_converged_reason -snes_view -show_solution 0',
                                                                 'setup': './bin/pythonscripts/PetscGenerateFEMQuadrature.py 2 2 2 1 laplacian 2 1 1 1 gradient src/snes/examples/tutorials/ex62.h'},
                                                                #  Block diagonal \begin{pmatrix} A & 0 \\ 0 & I \end{pmatrix}
-                                                               {'numProcs': 1, 'args': '-run_type full -refinement_limit 0.00625 -bc_type dirichlet -interpolate 1 -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -pc_type fieldsplit -pc_fieldsplit_type additive -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_pc_type jacobi -snes_monitor_short -ksp_monitor_short -snes_converged_reason -snes_view -show_solution 0'},
+                                                               {'numProcs': 1, 'args': '-run_type full -refinement_limit 0.00625 -bc_type dirichlet -interpolate 1 -order 2,1 -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -pc_type fieldsplit -pc_fieldsplit_type additive -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_pc_type jacobi -snes_monitor_short -ksp_monitor_short -snes_converged_reason -snes_view -show_solution 0'},
                                                                #  Block triangular \begin{pmatrix} A & B \\ 0 & I \end{pmatrix}
-                                                               {'numProcs': 1, 'args': '-run_type full -refinement_limit 0.00625 -bc_type dirichlet -interpolate 1 -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -pc_type fieldsplit -pc_fieldsplit_type multiplicative -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_pc_type jacobi -snes_monitor_short -ksp_monitor_short -snes_converged_reason -snes_view -show_solution 0'},
+                                                               {'numProcs': 1, 'args': '-run_type full -refinement_limit 0.00625 -bc_type dirichlet -interpolate 1 -order 2,1 -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -pc_type fieldsplit -pc_fieldsplit_type multiplicative -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_pc_type jacobi -snes_monitor_short -ksp_monitor_short -snes_converged_reason -snes_view -show_solution 0'},
                                                                #  Diagonal Schur complement \begin{pmatrix} A & 0 \\ 0 & S \end{pmatrix}
-                                                               {'numProcs': 1, 'args': '-run_type full -refinement_limit 0.00625 -bc_type dirichlet -interpolate 1 -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type diag -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_velocity_ksp_type gmres -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_pc_type jacobi -snes_monitor_short -ksp_monitor_short -snes_converged_reason -snes_view -show_solution 0'},
+                                                               {'numProcs': 1, 'args': '-run_type full -refinement_limit 0.00625 -bc_type dirichlet -interpolate 1 -order 2,1 -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type diag -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_velocity_ksp_type gmres -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_pc_type jacobi -snes_monitor_short -ksp_monitor_short -snes_converged_reason -snes_view -show_solution 0'},
                                                                #  Upper triangular Schur complement \begin{pmatrix} A & B \\ 0 & S \end{pmatrix}
-                                                               {'numProcs': 1, 'args': '-run_type full -refinement_limit 0.00625 -bc_type dirichlet -interpolate 1 -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type upper -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_velocity_ksp_type gmres -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_pc_type jacobi -snes_monitor_short -ksp_monitor_short -snes_converged_reason -snes_view -show_solution 0'},
+                                                               {'numProcs': 1, 'args': '-run_type full -refinement_limit 0.00625 -bc_type dirichlet -interpolate 1 -order 2,1 -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type upper -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_velocity_ksp_type gmres -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_pc_type jacobi -snes_monitor_short -ksp_monitor_short -snes_converged_reason -snes_view -show_solution 0'},
                                                                #  Lower triangular Schur complement \begin{pmatrix} A & B \\ 0 & S \end{pmatrix}
-                                                               {'numProcs': 1, 'args': '-run_type full -refinement_limit 0.00625 -bc_type dirichlet -interpolate 1 -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type lower -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_velocity_ksp_type gmres -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_pc_type jacobi -snes_monitor_short -ksp_monitor_short -snes_converged_reason -snes_view -show_solution 0'},
+                                                               {'numProcs': 1, 'args': '-run_type full -refinement_limit 0.00625 -bc_type dirichlet -interpolate 1 -order 2,1 -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type lower -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_velocity_ksp_type gmres -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_pc_type jacobi -snes_monitor_short -ksp_monitor_short -snes_converged_reason -snes_view -show_solution 0'},
                                                                #  Full Schur complement \begin{pmatrix} I & 0 \\ B^T A^{-1} & I \end{pmatrix} \begin{pmatrix} A & 0 \\ 0 & S \end{pmatrix} \begin{pmatrix} I & A^{-1} B \\ 0 & I \end{pmatrix}
-                                                               {'numProcs': 1, 'args': '-run_type full -refinement_limit 0.00625 -bc_type dirichlet -interpolate 1 -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_velocity_ksp_type gmres -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_pc_type jacobi -snes_monitor_short -ksp_monitor_short -snes_converged_reason -snes_view -show_solution 0'},
+                                                               {'numProcs': 1, 'args': '-run_type full -refinement_limit 0.00625 -bc_type dirichlet -interpolate 1 -order 2,1 -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_velocity_ksp_type gmres -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_pc_type jacobi -snes_monitor_short -ksp_monitor_short -snes_converged_reason -snes_view -show_solution 0'},
                                                                #  SIMPLE \begin{pmatrix} I & 0 \\ B^T A^{-1} & I \end{pmatrix} \begin{pmatrix} A & 0 \\ 0 & B^T diag(A)^{-1} B \end{pmatrix} \begin{pmatrix} I & diag(A)^{-1} B \\ 0 & I \end{pmatrix}
-                                                               #{'numProcs': 1, 'args': '-run_type full -refinement_limit 0.00625 -bc_type dirichlet -interpolate 1 -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_velocity_ksp_type gmres -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_pc_type jacobi -fieldsplit_pressure_inner_ksp_type preonly -fieldsplit_pressure_inner_pc_type jacobi -fieldsplit_pressure_upper_ksp_type preonly -fieldsplit_pressure_upper_pc_type jacobi -snes_monitor_short -ksp_monitor_short -snes_converged_reason -snes_view -show_solution 0'},
+                                                               #{'numProcs': 1, 'args': '-run_type full -refinement_limit 0.00625 -bc_type dirichlet -interpolate 1 -order 2,1 -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_velocity_ksp_type gmres -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_pc_type jacobi -fieldsplit_pressure_inner_ksp_type preonly -fieldsplit_pressure_inner_pc_type jacobi -fieldsplit_pressure_upper_ksp_type preonly -fieldsplit_pressure_upper_pc_type jacobi -snes_monitor_short -ksp_monitor_short -snes_converged_reason -snes_view -show_solution 0'},
                                                                #  SIMPLEC \begin{pmatrix} I & 0 \\ B^T A^{-1} & I \end{pmatrix} \begin{pmatrix} A & 0 \\ 0 & B^T rowsum(A)^{-1} B \end{pmatrix} \begin{pmatrix} I & rowsum(A)^{-1} B \\ 0 & I \end{pmatrix}
-                                                               #{'numProcs': 1, 'args': '-run_type full -refinement_limit 0.00625 -bc_type dirichlet -interpolate 1 -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_velocity_ksp_type gmres -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_pc_type jacobi -fieldsplit_pressure_inner_ksp_type preonly -fieldsplit_pressure_inner_pc_type jacobi -fieldsplit_pressure_inner_pc_jacobi_rowsum -fieldsplit_pressure_upper_ksp_type preonly -fieldsplit_pressure_upper_pc_type jacobi -fieldsplit_pressure_upper_pc_jacobi_rowsum -snes_monitor_short -ksp_monitor_short -snes_converged_reason -snes_view -show_solution 0'},
+                                                               #{'numProcs': 1, 'args': '-run_type full -refinement_limit 0.00625 -bc_type dirichlet -interpolate 1 -order 2,1 -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_velocity_ksp_type gmres -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_pc_type jacobi -fieldsplit_pressure_inner_ksp_type preonly -fieldsplit_pressure_inner_pc_type jacobi -fieldsplit_pressure_inner_pc_jacobi_rowsum -fieldsplit_pressure_upper_ksp_type preonly -fieldsplit_pressure_upper_pc_type jacobi -fieldsplit_pressure_upper_pc_jacobi_rowsum -snes_monitor_short -ksp_monitor_short -snes_converged_reason -snes_view -show_solution 0'},
                                                                # Stokes preconditioners with MF Jacobian action 37-42
                                                                #   Jacobi
-                                                               {'numProcs': 1, 'args': '-run_type full -refinement_limit 0.00625 -bc_type dirichlet -interpolate 1 -jacobian_mf -ksp_gmres_restart 100 -pc_type jacobi -ksp_rtol 1.0e-9 -snes_monitor_short -ksp_monitor_short -snes_converged_reason -snes_view -show_solution 0'},
+                                                               {'numProcs': 1, 'args': '-run_type full -refinement_limit 0.00625 -bc_type dirichlet -interpolate 1 -order 2,1 -jacobian_mf -ksp_gmres_restart 100 -pc_type jacobi -ksp_rtol 1.0e-9 -snes_monitor_short -ksp_monitor_short -snes_converged_reason -snes_view -show_solution 0'},
                                                                #  Block diagonal \begin{pmatrix} A & 0 \\ 0 & I \end{pmatrix}
-                                                               {'numProcs': 1, 'args': '-run_type full -refinement_limit 0.00625 -bc_type dirichlet -interpolate 1 -jacobian_mf -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -pc_type fieldsplit -pc_fieldsplit_type additive -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_pc_type jacobi -snes_monitor_short -ksp_monitor_short -snes_converged_reason -snes_view -show_solution 0'},
+                                                               {'numProcs': 1, 'args': '-run_type full -refinement_limit 0.00625 -bc_type dirichlet -interpolate 1 -order 2,1 -jacobian_mf -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -pc_type fieldsplit -pc_fieldsplit_type additive -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_pc_type jacobi -snes_monitor_short -ksp_monitor_short -snes_converged_reason -snes_view -show_solution 0'},
                                                                #  Block triangular \begin{pmatrix} A & B \\ 0 & I \end{pmatrix}
-                                                               {'numProcs': 1, 'args': '-run_type full -refinement_limit 0.00625 -bc_type dirichlet -interpolate 1 -jacobian_mf -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -pc_type fieldsplit -pc_fieldsplit_type multiplicative -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_pc_type jacobi -snes_monitor_short -ksp_monitor_short -snes_converged_reason -snes_view -show_solution 0'},
+                                                               {'numProcs': 1, 'args': '-run_type full -refinement_limit 0.00625 -bc_type dirichlet -interpolate 1 -order 2,1 -jacobian_mf -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -pc_type fieldsplit -pc_fieldsplit_type multiplicative -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_pc_type jacobi -snes_monitor_short -ksp_monitor_short -snes_converged_reason -snes_view -show_solution 0'},
                                                                #  Diagonal Schur complement \begin{pmatrix} A & 0 \\ 0 & S \end{pmatrix}
-                                                               {'numProcs': 1, 'args': '-run_type full -refinement_limit 0.00625 -bc_type dirichlet -interpolate 1 -jacobian_mf -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type diag -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_velocity_ksp_type gmres -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_pc_type jacobi -snes_monitor_short -ksp_monitor_short -snes_converged_reason -snes_view -show_solution 0'},
+                                                               {'numProcs': 1, 'args': '-run_type full -refinement_limit 0.00625 -bc_type dirichlet -interpolate 1 -order 2,1 -jacobian_mf -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type diag -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_velocity_ksp_type gmres -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_pc_type jacobi -snes_monitor_short -ksp_monitor_short -snes_converged_reason -snes_view -show_solution 0'},
                                                                #  Upper triangular Schur complement \begin{pmatrix} A & B \\ 0 & S \end{pmatrix}
-                                                               {'numProcs': 1, 'args': '-run_type full -refinement_limit 0.00625 -bc_type dirichlet -interpolate 1 -jacobian_mf -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type upper -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_velocity_ksp_type gmres -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_pc_type jacobi -snes_monitor_short -ksp_monitor_short -snes_converged_reason -snes_view -show_solution 0'},
+                                                               {'numProcs': 1, 'args': '-run_type full -refinement_limit 0.00625 -bc_type dirichlet -interpolate 1 -order 2,1 -jacobian_mf -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type upper -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_velocity_ksp_type gmres -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_pc_type jacobi -snes_monitor_short -ksp_monitor_short -snes_converged_reason -snes_view -show_solution 0'},
                                                                #  Lower triangular Schur complement \begin{pmatrix} A & B \\ 0 & S \end{pmatrix}
-                                                               #{'numProcs': 1, 'args': '-run_type full -refinement_limit 0.00625 -bc_type dirichlet -interpolate 1 -jacobian_mf -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type lower -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_velocity_ksp_type gmres -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_pc_type jacobi -snes_monitor_short -ksp_monitor_short -snes_converged_reason -snes_view -show_solution 0'},
+                                                               #{'numProcs': 1, 'args': '-run_type full -refinement_limit 0.00625 -bc_type dirichlet -interpolate 1 -order 2,1 -jacobian_mf -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type lower -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_velocity_ksp_type gmres -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_pc_type jacobi -snes_monitor_short -ksp_monitor_short -snes_converged_reason -snes_view -show_solution 0'},
                                                                #  Full Schur complement \begin{pmatrix} A & B \\ B^T & S \end{pmatrix}
-                                                               {'numProcs': 1, 'args': '-run_type full -refinement_limit 0.00625 -bc_type dirichlet -interpolate 1 -jacobian_mf -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_velocity_ksp_type gmres -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_pc_type jacobi -snes_monitor_short -ksp_monitor_short -snes_converged_reason -snes_view -show_solution 0'},
+                                                               {'numProcs': 1, 'args': '-run_type full -refinement_limit 0.00625 -bc_type dirichlet -interpolate 1 -order 2,1 -jacobian_mf -ksp_type fgmres -ksp_gmres_restart 100 -ksp_rtol 1.0e-9 -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_velocity_ksp_type gmres -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_pc_type jacobi -snes_monitor_short -ksp_monitor_short -snes_converged_reason -snes_view -show_solution 0'},
                                                                # 3D serial P1 tests 43-45
                                                                {'numProcs': 1, 'args': '-run_type test -dim 3 -refinement_limit 0.0    -bc_type dirichlet -interpolate 0 -show_initial -dm_plex_print_fem 1',
                                                                 'setup': './bin/pythonscripts/PetscGenerateFEMQuadrature.py 3 1 3 1 laplacian 3 1 1 1 gradient src/snes/examples/tutorials/ex62.h'},
                                                                #{'numProcs': 3, 'args': '-run_type full -refinement_limit 0.0625 -bc_type dirichlet -pc_type jacobi -ksp_rtol 1.0e-9 -snes_converged_reason -snes_view'},
                                                                #{'numProcs': 5, 'args': '-run_type full -refinement_limit 0.0625 -bc_type dirichlet -pc_type jacobi -ksp_rtol 1.0e-9 -snes_converged_reason -snes_view'}
 ],
-                        'src/ts/examples/tutorials/ex11':      [{'numProcs': 1, 'args': '-ufv_vtk_interval 0 -f %(meshes)s/sevenside.exo'},
+                        'src/ts/examples/tutorials/ex11':      [# 2D 0-6
+                                                                {'numProcs': 1, 'args': '-ufv_vtk_interval 0 -f %(meshes)s/sevenside.exo'},
                                                                 {'numProcs': 1, 'args': '-ufv_vtk_interval 0 -f %(meshes)s/sevenside-quad-15.exo'},
                                                                 {'numProcs': 2, 'args': '-ufv_vtk_interval 0 -f %(meshes)s/sevenside.exo'},
                                                                 {'numProcs': 2, 'args': '-ufv_vtk_interval 0 -f %(meshes)s/sevenside-quad-15.exo'},
                                                                 {'numProcs': 8, 'args': '-ufv_vtk_interval 0 -f %(meshes)s/sevenside-quad.exo'},
                                                                 {'numProcs': 1, 'args': '-ufv_vtk_interval 0 -f %(meshes)s/sevenside.exo -ts_type rosw'},
                                                                 {'numProcs': 1, 'args': '-ufv_vtk_interval 0 -f %(meshes)s/squaremotor-30.exo -ufv_split_faces'},
-                                                                {'numProcs': 1, 'args': '-ufv_vtk_interval 0 -f %(meshes)s/blockcylinder-50.exo'}],
+                                                                # 3D 7
+                                                                {'numProcs': 1, 'args': '-ufv_vtk_interval 0 -f %(meshes)s/blockcylinder-50.exo -bc_inflow 100,101,200 -bc_outflow 201'}],
                         }
 
 def noCheckCommand(command, status, output, error):

File config/builder2.py

   os.mkdir(objDir)
   executable  = os.path.join(objDir, exampleName)
   paramKey    = os.path.join(os.path.relpath(exampleDir, maker.petscDir), os.path.basename(executable))
-  if paramKey in builder.regressionRequirements:
-    if not builder.regressionRequirements[paramKey].issubset(packageNames):
-      raise RuntimeError('This test requires packages: %s' % builder.regressionRequirements[paramKey])
   params = builder.regressionParameters.get(paramKey, {})
   if not params:
     params = builder.getRegressionParameters(maker, exampleDir).get(paramKey, {})
   # Process testnum
   if args.testnum is not None:
     if args.testnum[0] == '[':
-      testnum = args.testnum[1:-1].split(',')
+      validTestnum = args.testnum[1:-1].split(',')
     else:
-      testnum = [args.testnum]
-    numtests = len(testnum)
+      validTestnum = [args.testnum]
+    numtests = len(validTestnum)
   else:
     numtests = len(params)
   rebuildTest = True
   maker.logPrint('Running %d tests\n' % (numtests,), debugSection='screen', forceScroll=True)
   for testnum, param in enumerate(params):
     testnum = str(testnum)
+    if 'requires' in param:
+      if not set(param['requires']).issubset(packageNames):
+        maker.logPrint('Test %s requires packages %s\n' % (testnum, param['requires']), debugSection='screen', forceScroll=True)
+        continue
     if 'num' in param: testnum = param['num']
-    if not args.testnum is None and not testnum in str(args.testnum): continue
+    if not args.testnum is None and not testnum in validTestnum: continue
     if 'setup' in param:
       print(param['setup'])
       os.system(sys.executable+' '+param['setup'])
         ex = [ex]+param['source']
       else:
         ex = ex+param['source']
+    # TODO: Fix this hack
+    if ex[-1] == 'F':
+      linkLanguage = 'FC'
+    else:
+      linkLanguage = maker.configInfo.languages.clanguage
     if rebuildTest:
       objects = maker.buildFile(ex, objDir)
       if not len(objects):
         print('TEST BUILD FAILED (check example.log for details)')
         return 1
-      maker.link(executable, objects, maker.configInfo.languages.clanguage)
+      maker.link(executable, objects, linkLanguage)
     if not 'args' in param: param['args'] = ''
     param['args'] += extraArgs
     if maker.runTest(exampleDir, executable, testnum, replace, **param):
   objDir        = maker.getObjDir(exampleName)
   executable    = os.path.join(objDir, exampleName)
   paramKey      = os.path.join(os.path.relpath(exampleDir, maker.petscDir), os.path.basename(executable))
-  if paramKey in builder.regressionRequirements:
-    if not builder.regressionRequirements[paramKey].issubset(packageNames):
-      raise RuntimeError('This test requires packages: %s' % builder.regressionRequirements[paramKey])
   params = builder.regressionParameters.get(paramKey, {})
   if not params:
     params = builder.getRegressionParameters(maker, exampleDir).get(paramKey, {})
     maker.logPrint('Makefile params '+str(params))
   if not isinstance(params, list):
     params = [params]
+  # Process testnum
+  if args.testnum is not None:
+    if args.testnum[0] == '[':
+      validTestnum = args.testnum[1:-1].split(',')
+    else:
+      validTestnum = [args.testnum]
+    numtests = len(validTestnum)
+  else:
+    numtests = len(params)
+  maker.logPrint('Running %d tests\n' % (numtests,), debugSection='screen', forceScroll=True)
   for testnum, param in enumerate(params):
+    testnum = str(testnum)
+    if 'requires' in param:
+      if not set(param['requires']).issubset(packageNames):
+        maker.logPrint('Test %s requires packages %s\n' % (testnum, param['requires']), debugSection='screen', forceScroll=True)
+        continue
     if 'num' in param: testnum = param['num']
-    if not args.testnum is None and testnum != args.testnum: continue
+    if not args.testnum is None and not testnum in validTestnum: continue
     if not 'args' in param: param['args'] = ''
     param['args'] += extraArgs
     print(str(testnum)+':  '+maker.getTestCommand(executable, **param))

File include/finclude/ftn-custom/petscdmplex.h90

       End Interface
 
       Interface
+        Subroutine DMPlexRestoreCone(m,p,cone,ierr)
+          USE_DMPLEX_HIDE
+          PetscInt     p
+          PetscInt, pointer :: cone(:)
+          PetscErrorCode ierr
+          DM_HIDE m
+        End Subroutine
+      End Interface
+
+      Interface
         Subroutine DMPlexGetConeOrientation(m,p,coneOrient,ierr)
           USE_DMPLEX_HIDE
           PetscInt     p
       End Interface
 
       Interface
+        Subroutine DMPlexRestoreSupport(m,p,support,ierr)
+          USE_DMPLEX_HIDE
+          PetscInt     p
+          PetscInt, pointer :: support(:)
+          PetscErrorCode ierr
+          DM_HIDE m
+        End Subroutine
+      End Interface
+
+      Interface
         Subroutine DMPlexGetTransitiveClosure(m,p,useCone,points,ierr)
           USE_DMPLEX_HIDE
           PetscInt     p
       End Interface
 
       Interface
+        Subroutine DMPlexMatSetClosure(m,section,globalSection,A,point, &
+     &  values,mode,ierr)
+          USE_DMPLEX_HIDE
+          PETSCSECTION_HIDE section
+          PETSCSECTION_HIDE globalSection
+          MAT_HIDE     A
+          PetscInt     point
+          InsertMode   mode
+          PetscScalar, pointer :: values(:)
+          PetscErrorCode ierr
+          DM_HIDE m
+        End Subroutine
+      End Interface
+
+      Interface
         Subroutine DMPlexCreateSection(m,dim,numFields,numComp,         &
      &  numDof,numBC,bcField,bcPoints,section,ierr)
           USE_DMPLEX_HIDE
           DM_HIDE m
         End Subroutine
       End Interface
+
+      Interface
+        Subroutine DMPlexComputeCellGeometry(m,cell,v0,J,invJ,detJ,ierr)
+          USE_DMPLEX_HIDE
+          PetscInt   cell
+          PetscReal, pointer :: v0(:)
+          PetscReal, pointer :: J(:)
+          PetscReal, pointer :: invJ(:)
+          PetscReal  detJ
+          PetscErrorCode ierr
+          DM_HIDE m
+        End Subroutine
+      End Interface

File include/finclude/petscvec.h

       PetscEnum INSERT_VALUES
       PetscEnum ADD_VALUES
       PetscEnum MAX_VALUES
+      PetscEnum INSERT_ALL_VALUES
+      PetscEnum ADD_ALL_VALUES
+      PetscEnum INSERT_BC_VALUES
+      PetscEnum ADD_BC_VALUES
 
       parameter (NOT_SET_VALUES=0,INSERT_VALUES=1,ADD_VALUES=2)
       parameter (MAX_VALUES=3)
+      parameter (INSERT_ALL_VALUES=4,ADD_ALL_VALUES=5)
+      parameter (INSERT_BC_VALUES=6,ADD_BC_VALUES=7)
 !
 !  Types of vector scatters
 !

File include/petsc-private/dmpleximpl.h

                                          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 (*integrateBdResidualFEM)(PetscInt, PetscInt, PetscInt, PetscQuadrature[], const PetscScalar[],
+                                           const PetscReal[], const PetscReal[], const PetscReal[], const PetscReal[], const PetscReal[],
+                                           void (*)(const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]),
+                                           void (*)(const PetscScalar[], const PetscScalar[], const PetscReal[], 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[]),

File include/petscdm.h

   void (**g1Funcs)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]); /* The g_1 functions for each field pair */
   void (**g2Funcs)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]); /* The g_2 functions for each field pair */
   void (**g3Funcs)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]); /* The g_3 functions for each field pair */
-  PetscScalar (**bcFuncs)(const PetscReal x[]); /* The boundary condition function for each field component */
+  void (**bcFuncs)(const PetscReal[], PetscScalar *); /* The boundary condition function for each field component */
+  PetscQuadrature *quadBd;
+  void (**f0BdFuncs)(const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]); /* The f_0 functions for each field */
+  void (**f1BdFuncs)(const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]); /* The f_1 functions for each field */
 } PetscFEM;
 
 typedef enum {PETSC_UNIT_LENGTH, PETSC_UNIT_MASS, PETSC_UNIT_TIME, PETSC_UNIT_CURRENT, PETSC_UNIT_TEMPERATURE, PETSC_UNIT_AMOUNT, PETSC_UNIT_LUMINOSITY, NUM_PETSC_UNITS} PetscUnit;

File include/petscdmplex.h

 /* Support for cell-vertex meshes */
 PETSC_EXTERN PetscErrorCode DMPlexGetNumFaceVertices(DM, PetscInt, PetscInt, PetscInt *);
 
+/* FVM Support */
+PETSC_EXTERN PetscErrorCode DMPlexComputeCellGeometryFVM(DM, PetscInt, PetscReal *, PetscReal [], PetscReal []);
+
 /* FEM Support */
 PETSC_EXTERN PetscErrorCode DMPlexComputeCellGeometry(DM, PetscInt, PetscReal *, PetscReal *, PetscReal *, PetscReal *);
 PETSC_EXTERN PetscErrorCode DMPlexVecGetClosure(DM, PetscSection, Vec, PetscInt, PetscInt *, PetscScalar *[]);
 PETSC_EXTERN PetscErrorCode DMPlexMatSetClosure(DM, PetscSection, PetscSection, Mat, PetscInt, const PetscScalar[], InsertMode);
 
 PETSC_EXTERN PetscErrorCode DMPlexCreateExodus(MPI_Comm, PetscInt, PetscBool, DM *);
+PETSC_EXTERN PetscErrorCode DMPlexCreateCGNS(MPI_Comm, PetscInt, PetscBool, DM *);
 
 PETSC_EXTERN PetscErrorCode DMPlexConstructGhostCells(DM, const char [], PetscInt *, DM *);
 PETSC_EXTERN PetscErrorCode DMPlexConstructCohesiveCells(DM, DMLabel, DM *);
   void *user;
 } JacActionCtx;
 
-PETSC_EXTERN PetscErrorCode DMPlexProjectFunction(DM, PetscInt, PetscScalar (**)(const PetscReal []), InsertMode, Vec);
-PETSC_EXTERN PetscErrorCode DMPlexProjectFunctionLocal(DM, PetscInt, PetscScalar (**)(const PetscReal []), InsertMode, Vec);
-PETSC_EXTERN PetscErrorCode DMPlexComputeL2Diff(DM, PetscQuadrature[], PetscScalar (**)(const PetscReal []), Vec, PetscReal *);
+PETSC_EXTERN PetscErrorCode DMPlexProjectFunction(DM, PetscInt, void (**)(const PetscReal [], PetscScalar *), InsertMode, Vec);
+PETSC_EXTERN PetscErrorCode DMPlexProjectFunctionLocal(DM, PetscInt, void (**)(const PetscReal [], PetscScalar *), InsertMode, Vec);
+PETSC_EXTERN PetscErrorCode DMPlexComputeL2Diff(DM, PetscQuadrature[], void (**)(const PetscReal [], PetscScalar *), Vec, PetscReal *);
 PETSC_EXTERN PetscErrorCode DMPlexComputeResidualFEM(DM, Vec, Vec, void *);
 PETSC_EXTERN PetscErrorCode DMPlexComputeJacobianActionFEM(DM, Mat, Vec, Vec, void *);
 PETSC_EXTERN PetscErrorCode DMPlexComputeJacobianFEM(DM, Vec, Mat, Mat, MatStructure*,void *);
                                                                           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 (*)(PetscInt, PetscInt, PetscInt, PetscQuadrature[], const PetscScalar[],
+                                                                          const PetscReal[], const PetscReal[], const PetscReal[], const PetscReal[], const PetscReal[],
+                                                                          void (*)(const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]),
+                                                                          void (*)(const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]), PetscScalar[]),
                                                        PetscErrorCode (*)(PetscInt, PetscInt, PetscInt, PetscQuadrature[], const PetscScalar[], const PetscScalar[],
                                                                           const PetscReal[], const PetscReal[], const PetscReal[], const PetscReal[],
                                                                           void (**)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]),

File src/dm/dt/interface/dt.c

 }
 
 #undef __FUNCT__
+#define __FUNCT__ "PetscDTFactorial_Internal"
+/* Evaluates the nth jacobi polynomial with weight parameters a,b at a point x.
+   Recurrence relations implemented from the pseudocode given in Karniadakis and Sherwin, Appendix B */
+PETSC_STATIC_INLINE PetscErrorCode PetscDTFactorial_Internal(PetscInt n, PetscReal *factorial)
+{
+  PetscReal f = 1.0;
+  PetscInt  i;
+
+  PetscFunctionBegin;
+  for (i = 1; i < n+1; ++i) f *= i;
+  *factorial = f;
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "PetscDTComputeJacobi"
+/* Evaluates the nth jacobi polynomial with weight parameters a,b at a point x.
+   Recurrence relations implemented from the pseudocode given in Karniadakis and Sherwin, Appendix B */
+PETSC_STATIC_INLINE PetscErrorCode PetscDTComputeJacobi(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscReal *P)
+{
+  PetscReal apb, pn1, pn2;
+  PetscInt  k;
+
+  PetscFunctionBegin;
+  if (!n) {*P = 1.0; PetscFunctionReturn(0);}
+  if (n == 1) {*P = 0.5 * (a - b + (a + b + 2.0) * x); PetscFunctionReturn(0);}
+  apb = a + b;
+  pn2 = 1.0;
+  pn1 = 0.5 * (a - b + (apb + 2.0) * x);
+  *P  = 0.0;
+  for (k = 2; k < n+1; ++k) {
+    PetscReal a1 = 2.0 * k * (k + apb) * (2.0*k + apb - 2.0);
+    PetscReal a2 = (2.0 * k + apb - 1.0) * (a*a - b*b);
+    PetscReal a3 = (2.0 * k + apb - 2.0) * (2.0 * k + apb - 1.0) * (2.0 * k + apb);
+    PetscReal a4 = 2.0 * (k + a - 1.0) * (k + b - 1.0) * (2.0 * k + apb);
+
+    a2  = a2 / a1;
+    a3  = a3 / a1;
+    a4  = a4 / a1;
+    *P  = (a2 + a3 * x) * pn1 - a4 * pn2;
+    pn2 = pn1;
+    pn1 = *P;
+  }
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "PetscDTComputeJacobiDerivative"
+/* Evaluates the first derivative of P_{n}^{a,b} at a point x. */
+PETSC_STATIC_INLINE PetscErrorCode PetscDTComputeJacobiDerivative(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscReal *P)
+{
+  PetscReal      nP;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  if (!n) {*P = 0.0; PetscFunctionReturn(0);}
+  ierr = PetscDTComputeJacobi(a+1, b+1, n-1, x, &nP);CHKERRQ(ierr);
+  *P   = 0.5 * (a + b + n + 1) * nP;
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "PetscDTMapSquareToTriangle_Internal"
+/* Maps from [-1,1]^2 to the (-1,1) reference triangle */
+PETSC_STATIC_INLINE PetscErrorCode PetscDTMapSquareToTriangle_Internal(PetscReal x, PetscReal y, PetscReal *xi, PetscReal *eta)
+{
+  PetscFunctionBegin;
+  *xi  = 0.5 * (1.0 + x) * (1.0 - y) - 1.0;
+  *eta = y;
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "PetscDTMapCubeToTetrahedron_Internal"
+/* Maps from [-1,1]^2 to the (-1,1) reference triangle */
+PETSC_STATIC_INLINE PetscErrorCode PetscDTMapCubeToTetrahedron_Internal(PetscReal x, PetscReal y, PetscReal z, PetscReal *xi, PetscReal *eta, PetscReal *zeta)
+{
+  PetscFunctionBegin;
+  *xi   = 0.25 * (1.0 + x) * (1.0 - y) * (1.0 - z) - 1.0;
+  *eta  = 0.5  * (1.0 + y) * (1.0 - z) - 1.0;
+  *zeta = z;
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "PetscDTGaussJacobiQuadrature1D_Internal"
+static PetscErrorCode PetscDTGaussJacobiQuadrature1D_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscReal *x, PetscReal *w)
+{
+  PetscInt       maxIter = 100;
+  PetscReal      eps     = 1.0e-8;
+  PetscReal      a1, a2, a3, a4, a5, a6;
+  PetscInt       k;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+
+  a1      = pow(2, a+b+1);
+#if defined(PETSC_HAVE_TGAMMA)
+  a2      = tgamma(a + npoints + 1);
+  a3      = tgamma(b + npoints + 1);
+  a4      = tgamma(a + b + npoints + 1);
+#else
+  SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"tgamma() - math routine is unavailable.");
+#endif
+
+  ierr = PetscDTFactorial_Internal(npoints, &a5);CHKERRQ(ierr);
+  a6   = a1 * a2 * a3 / a4 / a5;
+  /* Computes the m roots of P_{m}^{a,b} on [-1,1] by Newton's method with Chebyshev points as initial guesses.
+   Algorithm implemented from the pseudocode given by Karniadakis and Sherwin and Python in FIAT */
+  for (k = 0; k < npoints; ++k) {
+    PetscReal r = -cos((2.0*k + 1.0) * PETSC_PI / (2.0 * npoints)), dP;
+    PetscInt  j;
+
+    if (k > 0) r = 0.5 * (r + x[k-1]);
+    for (j = 0; j < maxIter; ++j) {
+      PetscReal s = 0.0, delta, f, fp;
+      PetscInt  i;
+
+      for (i = 0; i < k; ++i) s = s + 1.0 / (r - x[i]);
+      ierr = PetscDTComputeJacobi(a, b, npoints, r, &f);CHKERRQ(ierr);
+      ierr = PetscDTComputeJacobiDerivative(a, b, npoints, r, &fp);CHKERRQ(ierr);
+      delta = f / (fp - f * s);
+      r     = r - delta;
+      if (fabs(delta) < eps) break;
+    }
+    x[k] = r;
+    ierr = PetscDTComputeJacobiDerivative(a, b, npoints, x[k], &dP);CHKERRQ(ierr);
+    w[k] = a6 / (1.0 - PetscSqr(x[k])) / PetscSqr(dP);
+  }
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "PetscDTGaussJacobiQuadrature"
+/*@C
+  PetscDTGaussJacobiQuadrature - create Gauss-Jacobi quadrature for a simplex
+
+  Not Collective
+
+  Input Arguments:
++ dim - The simplex dimension
+. npoints - number of points
+. a - left end of interval (often-1)
+- b - right end of interval (often +1)
+
+  Output Arguments:
++ points - quadrature points
+- weights - quadrature weights
+
+  Level: intermediate
+
+  References:
+  Karniadakis and Sherwin.
+  FIAT
+
+.seealso: PetscDTGaussQuadrature()
+@*/
+PetscErrorCode PetscDTGaussJacobiQuadrature(PetscInt dim, PetscInt npoints, PetscReal a, PetscReal b, PetscReal *points[], PetscReal *weights[])
+{
+  PetscReal     *px, *wx, *py, *wy, *pz, *wz, *x, *w;
+  PetscInt       i, j, k;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  if ((a != -1.0) || (b != 1.0)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must use default internal right now");
+  switch (dim) {
+  case 1:
+    ierr = PetscMalloc(npoints * sizeof(PetscReal), &x);CHKERRQ(ierr);
+    ierr = PetscMalloc(npoints * sizeof(PetscReal), &w);CHKERRQ(ierr);
+    ierr = PetscDTGaussJacobiQuadrature1D_Internal(npoints, 0.0, 0.0, x, w);CHKERRQ(ierr);
+    break;
+  case 2:
+    ierr = PetscMalloc(npoints*npoints*2 * sizeof(PetscReal), &x);CHKERRQ(ierr);
+    ierr = PetscMalloc(npoints*npoints   * sizeof(PetscReal), &w);CHKERRQ(ierr);
+    ierr = PetscMalloc4(npoints,PetscReal,&px,npoints,PetscReal,&wx,npoints,PetscReal,&py,npoints,PetscReal,&wy);CHKERRQ(ierr);
+    ierr = PetscDTGaussJacobiQuadrature1D_Internal(npoints, 0.0, 0.0, px, wx);CHKERRQ(ierr);
+    ierr = PetscDTGaussJacobiQuadrature1D_Internal(npoints, 1.0, 0.0, py, wy);CHKERRQ(ierr);
+    for (i = 0; i < npoints; ++i) {
+      for (j = 0; j < npoints; ++j) {
+        ierr = PetscDTMapSquareToTriangle_Internal(px[i], py[j], &x[(i*npoints+j)*2+0], &x[(i*npoints+j)*2+1]);CHKERRQ(ierr);
+        w[i*npoints+j] = 0.5 * wx[i] * wy[j];
+      }
+    }
+    ierr = PetscFree4(px,wx,py,wy);CHKERRQ(ierr);
+    break;
+  case 3:
+    ierr = PetscMalloc(npoints*npoints*3 * sizeof(PetscReal), &x);CHKERRQ(ierr);
+    ierr = PetscMalloc(npoints*npoints   * sizeof(PetscReal), &w);CHKERRQ(ierr);
+    ierr = PetscMalloc6(npoints,PetscReal,&px,npoints,PetscReal,&wx,npoints,PetscReal,&py,npoints,PetscReal,&wy,npoints,PetscReal,&pz,npoints,PetscReal,&wz);CHKERRQ(ierr);
+    ierr = PetscDTGaussJacobiQuadrature1D_Internal(npoints, 0.0, 0.0, px, wx);CHKERRQ(ierr);
+    ierr = PetscDTGaussJacobiQuadrature1D_Internal(npoints, 1.0, 0.0, py, wy);CHKERRQ(ierr);
+    ierr = PetscDTGaussJacobiQuadrature1D_Internal(npoints, 2.0, 0.0, pz, wz);CHKERRQ(ierr);
+    for (i = 0; i < npoints; ++i) {
+      for (j = 0; j < npoints; ++j) {
+        for (k = 0; k < npoints; ++k) {
+          ierr = PetscDTMapCubeToTetrahedron_Internal(px[i], py[j], pz[k], &x[((i*npoints+j)*npoints+k)*3+0], &x[((i*npoints+j)*npoints+k)*3+1], &x[((i*npoints+j)*npoints+k)*3+2]);CHKERRQ(ierr);
+          w[(i*npoints+j)*npoints+k] = 0.125 * wx[i] * wy[j] * wz[k];
+        }
+      }
+    }
+    ierr = PetscFree6(px,wx,py,wy,pz,wz);CHKERRQ(ierr);
+    break;
+  default:
+    SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct quadrature rule for dimension %d", dim);
+  }
+  if (points)  *points  = x;
+  if (weights) *weights = w;
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
 #define __FUNCT__ "PetscDTPseudoInverseQR"
 /* Overwrites A. Can only handle full-rank problems with m>=n
  * A in column-major format

File src/dm/impls/plex/examples/tests/ex1.c

 
 #include <petscdmplex.h>
 
+#if defined(PETSC_HAVE_CGNS)
+#include <cgnslib.h>
+#endif
+
 typedef struct {
   DM            dm;                /* REQUIRED in order to use SNES evaluation functions */
   PetscInt      debug;             /* The debugging level */
   PetscBool refinementUniform;     /* Uniformly refine the mesh */
   PetscReal refinementLimit;       /* The largest allowable cell volume */
   PetscBool cellSimplex;           /* Use simplices or hexes */
+  char      filename[PETSC_MAX_PATH_LEN]; /* Import mesh from file */
 } AppCtx;
 
 #undef __FUNCT__
   options->refinementUniform = PETSC_FALSE;
   options->refinementLimit   = 0.0;
   options->cellSimplex       = PETSC_TRUE;
+  options->filename[0]       = '\0';
 
   ierr = PetscOptionsBegin(comm, "", "Meshing Problem Options", "DMPLEX");CHKERRQ(ierr);
   ierr = PetscOptionsInt("-debug", "The debugging level", "ex1.c", options->debug, &options->debug, NULL);CHKERRQ(ierr);
   ierr = PetscOptionsBool("-refinement_uniform", "Uniformly refine the mesh", "ex1.c", options->refinementUniform, &options->refinementUniform, NULL);CHKERRQ(ierr);
   ierr = PetscOptionsReal("-refinement_limit", "The largest allowable cell volume", "ex1.c", options->refinementLimit, &options->refinementLimit, NULL);CHKERRQ(ierr);
   ierr = PetscOptionsBool("-cell_simplex", "Use simplices if true, otherwise hexes", "ex1.c", options->cellSimplex, &options->cellSimplex, NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsString("-filename", "The mesh file", "ex7.c", options->filename, options->filename, PETSC_MAX_PATH_LEN, NULL);CHKERRQ(ierr);
   ierr = PetscOptionsEnd();
 
   ierr = PetscLogEventRegister("CreateMesh",          DM_CLASSID,   &options->createMeshEvent);CHKERRQ(ierr);
   PetscBool      refinementUniform = user->refinementUniform;
   PetscReal      refinementLimit   = user->refinementLimit;
   PetscBool      cellSimplex       = user->cellSimplex;
-  const char     *partitioner      = "chaco";
+  const char    *filename          = user->filename;
+  const char    *partitioner       = "chaco";
+  size_t         len;
+  PetscMPIInt    rank;
   PetscErrorCode ierr;
 
   PetscFunctionBegin;
   ierr = PetscLogEventBegin(user->createMeshEvent,0,0,0,0);CHKERRQ(ierr);
-  if (cellSimplex) {
+  ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
+  ierr = PetscStrlen(filename, &len);CHKERRQ(ierr);
+  if (len) {
+#if defined(PETSC_HAVE_CGNS)
+    int cgid = -1;
+
+    if (!rank) {
+      ierr = cg_open(filename, CG_MODE_READ, &cgid);CHKERRQ(ierr);
+      if (cgid <= 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_LIB, "cg_open(\"%s\",...) did not return a valid file ID", filename);
+    }
+    ierr = DMPlexCreateCGNS(comm, cgid, interpolate, dm);CHKERRQ(ierr);
+    if (!rank) {ierr = cg_close(cgid);CHKERRQ(ierr);}
+#else
+    SETERRQ(comm, PETSC_ERR_SUP, "Loading meshes requires CGNS support. Reconfigure using --with-cgns-dir");
+#endif
+  } else if (cellSimplex) {
     ierr = DMPlexCreateBoxMesh(comm, dim, interpolate, dm);CHKERRQ(ierr);
   } else {
     const PetscInt cells[3] = {2, 2, 2};

File src/dm/impls/plex/examples/tests/ex1f90.F

       DM dm
       PetscInt, target, dimension(4) :: EC
       PetscInt, pointer :: pEC(:)
+      PetscInt, pointer :: pES(:)
       PetscInt c,d, firstCell, numCells
+      PetscInt v, numVertices, numPoints
       PetscErrorCode ierr
 
       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
       call DMPlexCreate(PETSC_COMM_WORLD, dm, ierr)
       firstCell = 0
-      numCells = 5
-      call DMPlexSetChart(dm, 0, numCells, ierr)
+      numCells = 2
+      numVertices = 6
+      numPoints = numCells+numVertices
+      call DMPlexSetChart(dm, 0, numPoints, ierr)
       do c=firstCell,numCells-1
          call DMPlexSetConeSize(dm, c, 4, ierr)
       end do
       call DMSetUp(dm, ierr)
 
-      EC(1) = 1
-      EC(2) = 2
-      EC(3) = 3
-      EC(4) = 4
+      EC(1) = 2
+      EC(2) = 3
+      EC(3) = 4
+      EC(4) = 5
       pEC => EC
       c = 0
       write(*,*) 'cell',c,pEC
       call DMPlexGetCone(dm, c , pEC, ierr)
       CHKERRQ(ierr)
       write(*,*) 'cell',c,pEC
+      EC(1) = 4
+      EC(2) = 5
+      EC(3) = 6
+      EC(4) = 7
+      pEC => EC
+      c = 1
+      write(*,*) 'cell',c,pEC
+      call DMPlexSetCone(dm, c , pEC, ierr)
+      CHKERRQ(ierr)
+      call DMPlexGetCone(dm, c , pEC, ierr)
+      CHKERRQ(ierr)
+      write(*,*) 'cell',c,pEC
+      call DMPlexRestoreCone(dm, c , pEC, ierr)
+      CHKERRQ(ierr)
+
+      call DMPlexSymmetrize(dm, ierr)
+      call DMPlexStratify(dm, ierr)
+
+      v = 4
+      call DMPlexGetSupport(dm, v , pES, ierr)
+      CHKERRQ(ierr)
+      write(*,*) 'vertex',v,pES
+      call DMPlexRestoreSupport(dm, v , pES, ierr)
+      CHKERRQ(ierr)
 
       call DMDestroy(dm,ierr)
       call PetscFinalize(ierr)

File src/dm/impls/plex/examples/tests/ex3.c

   PetscInt order;                  /* Order of polynomials to test */
 } AppCtx;
 
-PetscScalar constant(const PetscReal coords[])
+void constant(const PetscReal coords[], PetscScalar *u)
 {
-  return 1.0;
+  *u = 1.0;
 }
 
-PetscScalar linear_x(const PetscReal coords[])
+void linear_x(const PetscReal coords[], PetscScalar *u)
 {
-  return coords[0];
+  *u = coords[0];
 }
-PetscScalar linear_y(const PetscReal coords[])
+void linear_y(const PetscReal coords[], PetscScalar *u)
 {
-  return coords[1];
+  *u = coords[1];
 }
-PetscScalar linear_z(const PetscReal coords[])
+void linear_z(const PetscReal coords[], PetscScalar *u)
 {
-  return coords[2];
+  *u = coords[2];
 }
 
-PetscScalar quadratic_xx(const PetscReal coords[])
+void quadratic_xx(const PetscReal coords[], PetscScalar *u)
 {
-  return coords[0]*coords[0];
+  *u = coords[0]*coords[0];
 }
-PetscScalar quadratic_xy(const PetscReal coords[])
+void quadratic_xy(const PetscReal coords[], PetscScalar *u)
 {
-  return coords[0]*coords[1];
+  *u = coords[0]*coords[1];
 }
-PetscScalar quadratic_yz(const PetscReal coords[])
+void quadratic_yz(const PetscReal coords[], PetscScalar *u)
 {
-  return coords[1]*coords[2];
+  *u = coords[1]*coords[2];
 }
-PetscScalar quadratic_zx(const PetscReal coords[])
+void quadratic_zx(const PetscReal coords[], PetscScalar *u)
 {
-  return coords[2]*coords[0];
+  *u = coords[2]*coords[0];
 }
 
 #undef __FUNCT__

File src/dm/impls/plex/examples/tests/ex7.c

 */
 
 #include <petscdmplex.h>
+#if defined(PETSC_HAVE_EXODUSII)
+#include <exodusII.h>
+#endif
+
 /* List of test meshes
 
 Triangle
  cell   5 _______    cell
  0    / | \      \       1
      /  |  \      \
-   18   |   19 13  22
-   / 9 17 10 \      \
+   17   |   18 13  22
+   / 8 19 10 \      \
   /     |     \      \
  2---14-|------4--20--6
   \     |     /      /
-   \ 8  | 7  /      /
+   \ 9  | 7  /      /
    16   |   15 11  21
      \  |  /      /
       \ | /      /
         3-------
 
+
 Quadrilateral
 -------------
 Test 0:
 
 typedef struct {
   DM        dm;
-  PetscInt  debug;       /* The debugging level */
-  PetscInt  testNum;     /* Indicates the mesh to create */
-  PetscInt  dim;         /* The topological mesh dimension */
-  PetscBool cellSimplex; /* Use simplices or hexes */
+  PetscInt  debug;                        /* The debugging level */
+  PetscInt  testNum;                      /* Indicates the mesh to create */
+  PetscInt  dim;                          /* The topological mesh dimension */
+  PetscBool cellSimplex;                  /* Use simplices or hexes */
+  PetscBool useGenerator;                 /* Construct mesh with a mesh generator */
+  char      filename[PETSC_MAX_PATH_LEN]; /* Import mesh from file */
 } AppCtx;
 
 #undef __FUNCT__
   PetscErrorCode ierr;
 
   PetscFunctionBegin;
-  options->debug       = 0;
-  options->testNum     = 0;
-  options->dim         = 2;
-  options->cellSimplex = PETSC_TRUE;
-
-  ierr = PetscOptionsBegin(comm, "", "Meshing Problem Options", "DMPLEX");CHKERRQ(ierr);
+  options->debug        = 0;
+  options->testNum      = 0;
+  options->dim          = 2;
+  options->cellSimplex  = PETSC_TRUE;
+  options->useGenerator = PETSC_FALSE;
+  options->filename[0]  = '\0';
+
+  ierr = PetscOptionsBegin(comm, "", "Meshing Interpolation Test Options", "DMPLEX");CHKERRQ(ierr);
   ierr = PetscOptionsInt("-debug", "The debugging level", "ex7.c", options->debug, &options->debug, NULL);CHKERRQ(ierr);
   ierr = PetscOptionsInt("-testnum", "The mesh to create", "ex7.c", options->testNum, &options->testNum, NULL);CHKERRQ(ierr);
   ierr = PetscOptionsInt("-dim", "The topological mesh dimension", "ex7.c", options->dim, &options->dim, NULL);CHKERRQ(ierr);
   ierr = PetscOptionsBool("-cell_simplex", "Use simplices if true, otherwise hexes", "ex7.c", options->cellSimplex, &options->cellSimplex, NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsBool("-use_generator", "Use a mesh generator to build the mesh", "ex7.c", options->useGenerator, &options->useGenerator, NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsString("-filename", "The mesh file", "ex7.c", options->filename, options->filename, PETSC_MAX_PATH_LEN, NULL);CHKERRQ(ierr);
   ierr = PetscOptionsEnd();
   PetscFunctionReturn(0);
 };
     {
       PetscInt    numPoints[2]        = {5, 2};
       PetscInt    coneSize[23]        = {4, 4, 0, 0, 0, 0, 0};
-      PetscInt    cones[8]            = {2, 3, 4, 5,  3, 6, 4, 5};
+      PetscInt    cones[8]            = {2, 4, 3, 5,  3, 4, 6, 5};
       PetscInt    coneOrientations[8] = {0, 0, 0, 0,  0, 0, 0, 0};
       PetscScalar vertexCoords[15]    = {0.0, 0.0, -0.5,  0.0, -0.5, 0.0,  1.0, 0.0, 0.0,  0.0, 0.5, 0.0,  0.0, 0.0, 0.5};
       PetscInt    markerPoints[8]     = {2, 1, 3, 1, 4, 1, 5, 1};
     {
       PetscInt    numPoints[2]         = {12, 2};
       PetscInt    coneSize[14]         = {8, 8, 0,0,0,0,0,0,0,0,0,0,0,0};
-      PetscInt    cones[16]            = {2,3,4,5,6,7,8,9,  3,10,11,4,7,12,13,8};
-      PetscInt    coneOrientations[16] = {0,0,0,0,0,0,0,0,  0, 0, 0,0,0, 0, 0,0};
+      PetscInt    cones[16]            = {2,5,4,3,6,7,8,9,  3,4,11,10,7,12,13,8};
+      PetscInt    coneOrientations[16] = {0,0,0,0,0,0,0,0,  0,0, 0, 0,0, 0, 0,0};
       PetscScalar vertexCoords[36]     = {-0.5,0.0,0.0, 0.0,0.0,0.0, 0.0,1.0,0.0, -0.5,1.0,0.0,
                                           -0.5,0.0,1.0, 0.0,0.0,1.0, 0.0,1.0,1.0, -0.5,1.0,1.0,
                                            0.5,0.0,0.0, 0.5,1.0,0.0, 0.5,0.0,1.0,  0.5,1.0,1.0};
 }
 
 #undef __FUNCT__
+#define __FUNCT__ "CheckMesh"
+PetscErrorCode CheckMesh(DM dm)
+{
+  PetscReal      detJ, J[9];
+  PetscInt       cStart, cEnd, c;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
+  for (c = cStart; c < cEnd; ++c) {
+    ierr = DMPlexComputeCellGeometry(dm, c, NULL, J, NULL, &detJ);CHKERRQ(ierr);
+    if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mesh cell %d is inverted, |J| = %g", c, detJ);
+  }
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "CompareCones"
+PetscErrorCode CompareCones(DM dm, DM idm)
+{
+  PetscInt       cStart, cEnd, c, vStart, vEnd, v;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
+  ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
+  for (c = cStart; c < cEnd; ++c) {
+    const PetscInt *cone;
+    PetscInt       *points = NULL, numPoints, p, numVertices = 0, coneSize;
+
+    ierr = DMPlexGetCone(dm, c, &cone);CHKERRQ(ierr);
+    ierr = DMPlexGetConeSize(dm, c, &coneSize);CHKERRQ(ierr);
+    ierr = DMPlexGetTransitiveClosure(idm, c, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr);
+    for (p = 0; p < numPoints*2; p += 2) {
+      const PetscInt point = points[p];
+      if ((point >= vStart) && (point < vEnd)) points[numVertices++] = point;
+    }
+    if (numVertices != coneSize) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In cell %d, cone size %d != %d vertices in closure", c, coneSize, numVertices);
+    for (v = 0; v < numVertices; ++v) {
+      if (cone[v] != points[v]) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In cell %d, cone point %d is %d != %d vertex in closure", c, v, cone[v], points[v]);
+    }
+    ierr = DMPlexRestoreTransitiveClosure(idm, c, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr);
+  }
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
 #define __FUNCT__ "CreateMesh"
 PetscErrorCode CreateMesh(MPI_Comm comm, PetscInt testNum, AppCtx *user, DM *dm)
 {
-  PetscInt       dim         = user->dim;
-  PetscBool      cellSimplex = user->cellSimplex;
-  const char    *partitioner = "chaco";
+  PetscInt       dim          = user->dim;
+  PetscBool      cellSimplex  = user->cellSimplex;
+  PetscBool      useGenerator = user->useGenerator;
+  const char    *filename     = user->filename;
+  const char    *partitioner  = "chaco";
+  size_t         len;
   PetscMPIInt    rank;
   PetscErrorCode ierr;
 
   PetscFunctionBegin;
   ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
-  ierr = DMCreate(comm, dm);CHKERRQ(ierr);
-  ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr);
-  ierr = DMPlexSetDimension(*dm, dim);CHKERRQ(ierr);
-  switch (dim) {
-  case 2:
-    if (cellSimplex) {
-      ierr = CreateSimplex_2D(comm, *dm);CHKERRQ(ierr);
-    } else {
-      ierr = CreateQuad_2D(comm, testNum, *dm);CHKERRQ(ierr);
+  ierr = PetscStrlen(filename, &len);CHKERRQ(ierr);
+  if (len) {
+#if defined(PETSC_HAVE_EXODUSII)
+    int   CPU_word_size = 0, IO_word_size = 0, exoid = -1;
+    float version;
+
+    if (!rank) {
+      exoid = ex_open(filename, EX_READ, &CPU_word_size, &IO_word_size, &version);
+      if (exoid <= 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_LIB, "ex_open(\"%s\",...) did not return a valid file ID", filename);
     }
-    break;
-  case 3:
+    ierr = DMPlexCreateExodus(comm, exoid, PETSC_FALSE, dm);CHKERRQ(ierr);
+    if (!rank) {ierr = ex_close(exoid);CHKERRQ(ierr);}
+    ierr = DMPlexGetDimension(*dm, &dim);CHKERRQ(ierr);
+#else
+    SETERRQ(comm, PETSC_ERR_SUP, "Loading meshes requires ExodusII support. Reconfigure using --download-exodusii");
+#endif
+  } else if (useGenerator) {
     if (cellSimplex) {
-      ierr = CreateSimplex_3D(comm, *dm);CHKERRQ(ierr);
+      ierr = DMPlexCreateBoxMesh(comm, dim, PETSC_FALSE, dm);CHKERRQ(ierr);
     } else {
-      ierr = CreateHex_3D(comm, *dm);CHKERRQ(ierr);
+      ierr = DMPlexCreateHexBoxMesh(comm, dim, PETSC_FALSE, dm);CHKERRQ(ierr);
+    }
+  } else {
+    ierr = DMCreate(comm, dm);CHKERRQ(ierr);
+    ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr);
+    ierr = DMPlexSetDimension(*dm, dim);CHKERRQ(ierr);
+    switch (dim) {
+    case 2:
+      if (cellSimplex) {
+        ierr = CreateSimplex_2D(comm, *dm);CHKERRQ(ierr);
+      } else {
+        ierr = CreateQuad_2D(comm, testNum, *dm);CHKERRQ(ierr);
+      }
+      break;
+    case 3:
+      if (cellSimplex) {
+        ierr = CreateSimplex_3D(comm, *dm);CHKERRQ(ierr);
+      } else {
+        ierr = CreateHex_3D(comm, *dm);CHKERRQ(ierr);
+      }
+      break;
+    default:
+      SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make meshes for dimension %d", dim);
     }
-    break;
-  default:
-    SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make meshes for dimension %d", dim);
   }
   {
     DM interpolatedMesh = NULL;
 
+    ierr = CheckMesh(*dm);CHKERRQ(ierr);
     ierr = DMPlexInterpolate(*dm, &interpolatedMesh);CHKERRQ(ierr);
     ierr = DMPlexCopyCoordinates(*dm, interpolatedMesh);CHKERRQ(ierr);
+    ierr = CompareCones(*dm, interpolatedMesh);CHKERRQ(ierr);
     ierr = DMDestroy(dm);CHKERRQ(ierr);
     *dm  = interpolatedMesh;
   }
   ierr = PetscInitialize(&argc, &argv, NULL, help);CHKERRQ(ierr);
   ierr = ProcessOptions(PETSC_COMM_WORLD, &user);CHKERRQ(ierr);
   ierr = CreateMesh(PETSC_COMM_WORLD, user.testNum, &user, &user.dm);CHKERRQ(ierr);
+  ierr = CheckMesh(user.dm);CHKERRQ(ierr);
   ierr = DMDestroy(&user.dm);CHKERRQ(ierr);
   ierr = PetscFinalize();
   return 0;

File src/dm/impls/plex/examples/tests/ex8.c

 #include <petscdmplex.h>
 
 #undef __FUNCT__
+#define __FUNCT__ "ChangeCoordinates"
+PetscErrorCode ChangeCoordinates(DM dm, PetscInt spaceDim, PetscScalar vertexCoords[])
+{
+  PetscSection   coordSection;
+  Vec            coordinates;
+  PetscScalar   *coords;
+  PetscInt       vStart, vEnd, v, d, coordSize;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
+  ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
+  ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr);
+  ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr);
+  ierr = PetscSectionSetChart(coordSection, vStart, vEnd);CHKERRQ(ierr);
+  for (v = vStart; v < vEnd; ++v) {
+    ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr);
+    ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr);
+  }
+  ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr);
+  ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr);
+  ierr = VecCreate(PetscObjectComm((PetscObject) dm), &coordinates);CHKERRQ(ierr);
+  ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr);
+  ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr);
+  ierr = VecSetFromOptions(coordinates);CHKERRQ(ierr);
+  ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr);
+  for (v = vStart; v < vEnd; ++v) {
+    PetscInt off;
+
+    ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr);
+    for (d = 0; d < spaceDim; ++d) {
+      coords[off+d] = vertexCoords[(v-vStart)*spaceDim+d];
+    }
+  }
+  ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr);
+  ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr);
+  ierr = VecDestroy(&coordinates);CHKERRQ(ierr);
+  ierr = DMSetFromOptions(dm);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "CheckFEMGeometry"
+PetscErrorCode CheckFEMGeometry(DM dm, PetscInt cell, PetscInt spaceDim, PetscReal v0Ex[], PetscReal JEx[], PetscReal invJEx[], PetscReal detJEx)
+{
+  PetscReal      v0[3], J[9], invJ[9], detJ;
+  PetscInt       d, i, j;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  ierr = DMPlexComputeCellGeometry(dm, cell, v0, J, invJ, &detJ);CHKERRQ(ierr);
+  for (d = 0; d < spaceDim; ++d) {
+    if (v0[d] != v0Ex[d]) {
+      switch (spaceDim) {
+      case 2: SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid v0 (%g, %g) != (%g, %g)", v0[0], v0[1], v0Ex[0], v0Ex[1]);break;
+      case 3: SETERRQ6(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid v0 (%g, %g, %g) != (%g, %g, %g)", v0[0], v0[1], v0[2], v0Ex[0], v0Ex[1], v0Ex[2]);break;
+      default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid space dimension %d", spaceDim);
+      }
+    }
+  }
+  for (i = 0; i < spaceDim; ++i) {
+    for (j = 0; j < spaceDim; ++j) {
+      if (fabs(J[i*spaceDim+j]    - JEx[i*spaceDim+j])    > 1.0e-9) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid J[%d,%d]: %g != %g", i, j, J[i*spaceDim+j], JEx[i*spaceDim+j]);
+      if (fabs(invJ[i*spaceDim+j] - invJEx[i*spaceDim+j]) > 1.0e-9) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid invJ[%d,%d]: %g != %g", i, j, invJ[i*spaceDim+j], invJEx[i*spaceDim+j]);
+    }
+  }
+  if (fabs(detJ - detJEx) > 1.0e-9) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid |J| = %g != %g", detJ, detJEx);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "CheckFVMGeometry"
+PetscErrorCode CheckFVMGeometry(DM dm, PetscInt cell, PetscInt spaceDim, PetscReal centroidEx[], PetscReal normalEx[], PetscReal volEx)
+{
+  PetscReal      centroid[3], normal[3], vol;
+  PetscInt       d;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  ierr = DMPlexComputeCellGeometryFVM(dm, cell, &vol, centroid, normal);CHKERRQ(ierr);
+  for (d = 0; d < spaceDim; ++d) {
+    if (fabs(centroid[d] - centroidEx[d]) > 1.0e-9) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid centroid[%d]: %g != %g", d, centroid[d], centroidEx[d]);
+    if (fabs(normal[d] - normalEx[d]) > 1.0e-9) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid normal[%d]: %g != %g", d, normal[d], normalEx[d]);
+  }
+  if (fabs(volEx - vol) > 1.0e-9) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid volume = %g != %g", vol, volEx);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
 #define __FUNCT__ "TestTriangle"
-PetscErrorCode TestTriangle(MPI_Comm comm)
+PetscErrorCode TestTriangle(MPI_Comm comm, PetscBool interpolate, PetscBool transform)
 {
   DM             dm;
-  PetscReal      v0[3], J[9], invJ[9], detJ;
-  PetscInt       dim, i, j;
+  PetscRandom    r, ang, ang2;
+  PetscInt       dim, t;
   PetscErrorCode ierr;
 
   PetscFunctionBegin;
   /* Create reference triangle */
   dim  = 2;
   ierr = DMCreate(comm, &dm);CHKERRQ(ierr);
+  ierr = PetscObjectSetName((PetscObject) dm, "triangle");CHKERRQ(ierr);
   ierr = DMSetType(dm, DMPLEX);CHKERRQ(ierr);
   ierr = DMPlexSetDimension(dm, dim);CHKERRQ(ierr);
   {
     PetscScalar vertexCoords[6]     = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0};
 
     ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
+    if (interpolate) {
+      DM idm;
+
+      ierr = DMPlexInterpolate(dm, &idm);CHKERRQ(ierr);
+      ierr = PetscObjectSetName((PetscObject) idm, "triangle");CHKERRQ(ierr);
+      ierr = DMPlexCopyCoordinates(dm, idm);CHKERRQ(ierr);
+      ierr = DMDestroy(&dm);CHKERRQ(ierr);
+      dm   = idm;
+    }
     ierr = DMSetFromOptions(dm);CHKERRQ(ierr);
   }
   /* Check reference geometry: determinant is scaled by reference volume (2.0) */
-  ierr = DMPlexComputeCellGeometry(dm, 0, v0, J, invJ, &detJ);CHKERRQ(ierr);
-  if ((v0[0] != -1.0) || (v0[1] != -1.0)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid v0 (%g, %g)", v0[0], v0[1]);
-  for (i = 0; i < dim; ++i) {
-    for (j = 0; j < dim; ++j) {
-      if (fabs(J[i*dim+j] - (i == j ? 1.0 : 0.0)) > 1.0e-9) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid J");
-      if (fabs(invJ[i*dim+j] - (i == j ? 1.0 : 0.0)) > 1.0e-9) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid invJ");
+  {
+    PetscReal v0Ex[2]       = {-1.0, -1.0};
+    PetscReal JEx[4]        = {1.0, 0.0, 0.0, 1.0};
+    PetscReal invJEx[4]     = {1.0, 0.0, 0.0, 1.0};
+    PetscReal detJEx        = 1.0;
+    PetscReal centroidEx[2] = {-0.333333333333, -0.333333333333};
+    PetscReal normalEx[2]   = {0.0, 0.0};
+    PetscReal volEx         = 2.0;
+
+    ierr = CheckFEMGeometry(dm, 0, dim, v0Ex, JEx, invJEx, detJEx);CHKERRQ(ierr);
+    if (interpolate) {ierr = CheckFVMGeometry(dm, 0, dim, centroidEx, normalEx, volEx);CHKERRQ(ierr);}
+  }
+  /* Check random triangles: rotate, scale, then translate */
+  if (transform) {
+    ierr = PetscRandomCreate(PETSC_COMM_SELF, &r);CHKERRQ(ierr);
+    ierr = PetscRandomSetFromOptions(r);CHKERRQ(ierr);
+    ierr = PetscRandomSetInterval(r, 0.0, 10.0);CHKERRQ(ierr);
+    ierr = PetscRandomCreate(PETSC_COMM_SELF, &ang);CHKERRQ(ierr);
+    ierr = PetscRandomSetFromOptions(ang);CHKERRQ(ierr);
+    ierr = PetscRandomSetInterval(ang, 0.0, 2*PETSC_PI);CHKERRQ(ierr);
+    for (t = 0; t < 100; ++t) {
+      PetscScalar vertexCoords[6] = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0}, trans[2];
+      PetscReal   v0Ex[2]         = {-1.0, -1.0};
+      PetscReal   JEx[4]          = {1.0, 0.0, 0.0, 1.0}, R[4], rot[2], rotM[4];
+      PetscReal   invJEx[4]       = {1.0, 0.0, 0.0, 1.0};
+      PetscReal   detJEx          = 1.0, scale, phi;
+      PetscReal   centroidEx[2]   = {-0.333333333333, -0.333333333333};
+      PetscReal   normalEx[2]     = {0.0, 0.0};
+      PetscReal   volEx           = 2.0;
+      PetscInt    d, e, f, p;
+
+      ierr = PetscRandomGetValueReal(r, &scale);CHKERRQ(ierr);
+      ierr = PetscRandomGetValueReal(ang, &phi);CHKERRQ(ierr);
+      R[0] = cos(phi); R[1] = -sin(phi);
+      R[2] = sin(phi); R[3] =  cos(phi);
+      for (p = 0; p < 3; ++p) {
+        for (d = 0; d < dim; ++d) {
+          for (e = 0, rot[d] = 0.0; e < dim; ++e) {
+            rot[d] += R[d*dim+e] * vertexCoords[p*dim+e];
+          }
+        }
+        for (d = 0; d < dim; ++d) vertexCoords[p*dim+d] = rot[d];
+      }
+      for (d = 0; d < dim; ++d) {
+        for (e = 0, rot[d] = 0.0; e < dim; ++e) {
+          rot[d] += R[d*dim+e] * centroidEx[e];
+        }
+      }
+      for (d = 0; d < dim; ++d) centroidEx[d] = rot[d];
+      for (d = 0; d < dim; ++d) {
+        for (e = 0; e < dim; ++e) {
+          for (f = 0, rotM[d*dim+e] = 0.0; f < dim; ++f) {
+            rotM[d*dim+e] += R[d*dim+f] * JEx[f*dim+e];
+          }
+        }
+      }
+      for (d = 0; d < dim; ++d) {
+        for (e = 0; e < dim; ++e) {
+          JEx[d*dim+e] = rot