Commits

Matt Knepley committed cf6d85e Merge

Merge branch 'knepley/fix-dm-vec-blocksize' into knepley/pylith

* knepley/fix-dm-vec-blocksize: (84 commits)
DM: Handle empty vectors correctly for blocksize determination
SNES & TS: remove stale comments about DMDAComputeFunction
DMPlex: Fixed #ifdef for HDF5
DMPlex ex2: Added an example that checks a mesh
DMPlex: #ifdef support functions for HDF5
DMPlex: Damn this branch to Hell
Doc: Bib fix
Doc: Added tutorial
TS: Fix output sequence
DMPlex: Check for HDF5 attributes before writing them
DMPlex: Label with timestep number since we do not have times
Viewer+HDF5: Added PetscViewerHDF5HasAttribute()
TS: Set output sequence for timesteps
DMPlex: Handle output sequences - Added static to some functions
DM: Add idea of an output sequence - Added DMGet/SetOutputSequenceNumber() - Added doc for outputDM
REVERT TS: Timestep gets set into DMTS
DMPlex: Fix up HDF5 segregation
added MPI_MAX_PROCESSOR_NAME for MPIUni fixed streams/process.py script to work cleanly with NPMAX=1
DMPlex:Fix compiler warnings
DMPlex: Replace MAXPATHLEN with PETSC_MAX_PATH_LEN
...

Comments (0)

Files changed (108)

config/BuildSystem/config/packages/__init__.py

-all = ['make','fblaslapack','f2cblaslapack','BlasLapack','exodusii', 'scientificpython', 'fiat', 'MOAB', 'MPI', 'netcdf', 'PETSc','boost', 'cusp', 'thrust', 'hdf5', 'netcdf-cxx']
+all = ['make','fblaslapack','f2cblaslapack','BlasLapack','exodusii', 'scientificpython', 'fiat', 'MOAB', 'MPI', 'netcdf', 'PETSc','boost', 'cusp', 'thrust', 'hdf5', 'netcdf-cxx', 'cgns']

config/BuildSystem/config/packages/cgns.py

+#!/usr/bin/env python
+import config.base
+import config.package
+import os
+
+class Configure(config.package.Package):
+  def __init__(self, framework):
+    config.package.Package.__init__(self, framework)
+    self.functions  = ['cg_close']
+    self.includes   = ['cgnslib.h']
+    self.liblist    = [['libcgns.a'],
+                       ['libcgns.a', 'libhdf5.a']] # Sometimes they over-link

config/builder.py

                                                                {'numProcs': 1, 'args': '-run_type full -refinement_limit 0.015625 -interpolate 1 -petscspace_order 2 -pc_type gamg -ksp_rtol 1.0e-10 -ksp_monitor_short -ksp_converged_reason -snes_monitor_short -snes_converged_reason'},
                                                                {'numProcs': 1, 'args': '-run_type full -refinement_limit 0.015625 -variable_coefficient nonlinear -interpolate 1 -petscspace_order 2 -pc_type svd -ksp_rtol 1.0e-10 -snes_monitor_short -snes_converged_reason'},
                                                                {'numProcs': 1, 'args': '-run_type full -refinement_limit 0.03125 -variable_coefficient nonlinear -interpolate 1 -petscspace_order 1 -snes_type fas -snes_fas_levels 2 -pc_type svd -ksp_rtol 1.0e-10 -fas_coarse_pc_type svd -fas_coarse_ksp_rtol 1.0e-10 -fas_coarse_snes_monitor_short -snes_monitor_short -snes_fas_monitor -snes_linesearch_type basic -fas_coarse_snes_linesearch_type basic -snes_converged_reason -dm_refine_hierarchy 1 -snes_view -fas_levels_1_snes_type newtonls -fas_levels_1_pc_type svd -fas_levels_1_ksp_rtol 1.0e-10 -fas_levels_1_snes_monitor'},
-                                                               {'numProcs': 1, 'args': '-run_type full -refinement_limit 0.0625 -variable_coefficient nonlinear -interpolate 1 -petscspace_order 1 -snes_type fas -snes_fas_levels 3 -pc_type svd -ksp_rtol 1.0e-10 -fas_coarse_pc_type svd -fas_coarse_ksp_rtol 1.0e-10 -fas_coarse_snes_monitor -snes_monitor_short -snes_fas_monitor -snes_linesearch_type basic -fas_coarse_snes_linesearch_type basic -snes_converged_reason -dm_refine_hierarchy 2 -dm_plex_print_fem 0 -snes_view -fas_levels_1_snes_type newtonls -fas_levels_1_pc_type svd -fas_levels_1_ksp_rtol 1.0e-10 -fas_levels_1_snes_monitor -fas_levels_2_snes_type newtonls -fas_levels_2_pc_type svd -fas_levels_2_ksp_rtol 1.0e-10 -fas_levels_2_snes_monitor'}],
+                                                               {'numProcs': 1, 'args': '-run_type full -refinement_limit 0.0625 -variable_coefficient nonlinear -interpolate 1 -petscspace_order 1 -snes_type fas -snes_fas_levels 3 -pc_type svd -ksp_rtol 1.0e-10 -fas_coarse_pc_type svd -fas_coarse_ksp_rtol 1.0e-10 -fas_coarse_snes_monitor -snes_monitor_short -snes_fas_monitor -snes_linesearch_type basic -fas_coarse_snes_linesearch_type basic -snes_converged_reason -dm_refine_hierarchy 2 -dm_plex_print_fem 0 -snes_view -fas_levels_1_snes_type newtonls -fas_levels_1_pc_type svd -fas_levels_1_ksp_rtol 1.0e-10 -fas_levels_1_snes_monitor -fas_levels_2_snes_type newtonls -fas_levels_2_pc_type svd -fas_levels_2_ksp_rtol 1.0e-10 -fas_levels_2_snes_monitor'},
+                                                               # Restarting 43
+                                                               {'numProcs': 1, 'args': '-run_type test -refinement_limit 0.0    -bc_type dirichlet -interpolate 1 -petscspace_order 1 -dm_view hdf5:sol.h5 -vec_view hdf5:sol.h5::append'},
+                                                               {'numProcs': 1, 'args': '-run_type test -refinement_limit 0.0    -bc_type dirichlet -interpolate 1 -petscspace_order 1 -f sol.h5 -restart'}],
                         'src/snes/examples/tutorials/ex31':   [# Decoupled field Dirichlet tests 0-6
                                                                {'numProcs': 1, 'args': '-run_type test -refinement_limit 0.0     -forcing_type constant -bc_type dirichlet -interpolate 1 -show_initial -dm_plex_print_fem 1',
                                                                 'setup': './bin/pythonscripts/PetscGenerateFEMQuadrature.py 2 2 2 1 laplacian 2 1 1 1 gradient 2 1 1 1 identity src/snes/examples/tutorials/ex31.h'},

include/finclude/ftn-custom/petscdmplex.h90

 
       Interface
         Subroutine DMPlexCreateSection(m,dim,numFields,numComp,         &
-     &  numDof,numBC,bcField,bcPoints,section,ierr)
+     &  numDof,numBC,bcField,bcPoints,perm,section,ierr)
           USE_DMPLEX_HIDE
           PETSCSECTION_HIDE section
           PetscInt     dim, numFields, numBC
           PetscInt, pointer :: numDof(:)
           PetscInt, pointer :: bcField(:)
           IS_HIDE,  pointer :: bcPoints(:)
+          IS_HIDE perm
           PetscErrorCode ierr
           DM_HIDE m
         End Subroutine
           DM_HIDE m
         End Subroutine
       End Interface
+
+      Interface
+        Subroutine DMPlexComputeCellGeometryFVM(m,cell,vol,centroid,    &
+     &    normal,ierr)
+          USE_DMPLEX_HIDE
+          PetscInt   cell
+          PetscReal  vol
+          PetscReal, pointer :: centroid(:)
+          PetscReal, pointer :: normal(:)
+          PetscErrorCode ierr
+          DM_HIDE m
+        End Subroutine
+      End Interface

include/mpiuni/mpi.h

 #endif
 extern int MPIUNI_Memcpy(void*,const void*,int);
 
+#define MPI_MAX_PROCESSOR_NAME 1024
 
 #define MPI_REQUEST_NULL     ((MPI_Request)0)
 #define MPI_GROUP_NULL       ((MPI_Group)0)

include/petsc-private/dmimpl.h

   /* Fields are represented by objects */
   PetscInt                numFields;
   PetscFE                *fields;
+  /* Output structures */
+  DM                      dmBC;                 /* The DM with boundary conditions in the global DM */
+  PetscInt                outputSequenceNum;    /* The current sequence number for output */
 
   PetscObject             dmksp,dmsnes,dmts;
 };

include/petsc-private/dmpleximpl.h

 PETSC_EXTERN PetscErrorCode DMPlexVTKGetCellType(DM,PetscInt,PetscInt,PetscInt*);
 PETSC_EXTERN PetscErrorCode VecView_Plex_Local(Vec,PetscViewer);
 PETSC_EXTERN PetscErrorCode VecView_Plex(Vec,PetscViewer);
+PETSC_EXTERN PetscErrorCode VecLoad_Plex_Local(Vec,PetscViewer);
+PETSC_EXTERN PetscErrorCode VecLoad_Plex(Vec,PetscViewer);
 
-PETSC_EXTERN PetscErrorCode DMPlexGetAdjacency_Internal(DM,PetscInt,PetscBool,PetscBool,PetscInt*,PetscInt[]);
+PETSC_EXTERN PetscErrorCode DMPlexGetAdjacency_Internal(DM,PetscInt,PetscBool,PetscBool,PetscInt*,PetscInt*[]);
 PETSC_EXTERN PetscErrorCode DMPlexGetFaces_Internal(DM,PetscInt,PetscInt,PetscInt*,PetscInt*,const PetscInt*[]);
 PETSC_EXTERN PetscErrorCode DMPlexGetRawFaces_Internal(DM,PetscInt,PetscInt,const PetscInt[], PetscInt*,PetscInt*,const PetscInt*[]);
 PETSC_EXTERN PetscErrorCode DMPlexRestoreFaces_Internal(DM,PetscInt,PetscInt,PetscInt*,PetscInt*,const PetscInt*[]);
 PETSC_EXTERN PetscErrorCode DMPlexGetCellRefiner_Internal(DM,CellRefiner*);
 PETSC_EXTERN PetscErrorCode CellRefinerGetAffineTransforms_Internal(CellRefiner, PetscInt *, PetscReal *[], PetscReal *[], PetscReal *[]);
 PETSC_EXTERN PetscErrorCode CellRefinerRestoreAffineTransforms_Internal(CellRefiner, PetscInt *, PetscReal *[], PetscReal *[], PetscReal *[]);
+PETSC_EXTERN PetscErrorCode DMPlexCreateGmsh_ReadElement(FILE *, PetscInt *, int *, PetscInt *, int[], int[]);
 
 #undef __FUNCT__
 #define __FUNCT__ "DMPlex_Invert2D_Internal"

include/petsc-private/isimpl.h

   PetscErrorCode (*duplicate)(IS,IS*);
   PetscErrorCode (*destroy)(IS);
   PetscErrorCode (*view)(IS,PetscViewer);
+  PetscErrorCode (*load)(IS,PetscViewer);
   PetscErrorCode (*identity)(IS,PetscBool*);
   PetscErrorCode (*copy)(IS,IS);
   PetscErrorCode (*togeneral)(IS);
 
 struct _p_IS {
   PETSCHEADER(struct _ISOps);
+  PetscLayout  map;
   PetscBool    isperm;          /* if is a permutation */
   PetscInt     max,min;         /* range of possible values */
-  PetscInt     bs;              /* block size */
   void         *data;
   PetscBool    isidentity;
   PetscInt     *total, *nonlocal;   /* local representation of ALL indices across the comm as well as the nonlocal part. */
   IS           complement;          /* IS wrapping nonlocal indices. */
 };
 
+extern PetscErrorCode ISLoad_Default(IS, PetscViewer);
+
 struct _p_ISLocalToGlobalMapping{
   PETSCHEADER(int);
   PetscInt n;                  /* number of local indices */
 };
 
 /* ----------------------------------------------------------------------------*/
-typedef struct _n_PetscUniformSection *PetscUniformSection;
-struct _n_PetscUniformSection {
-  MPI_Comm comm;
-  PetscInt pStart, pEnd; /* The chart: all points are contained in [pStart, pEnd) */
-  PetscInt numDof;       /* Describes layout of storage, point --> (constant # of values, (p - pStart)*constant # of values) */
-};
-
 struct _p_PetscSection {
   PETSCHEADER(int);
-  struct _n_PetscUniformSection atlasLayout;  /* Layout for the atlas */
+  PetscInt                      pStart, pEnd; /* The chart: all points are contained in [pStart, pEnd) */
+  IS                            perm;         /* A permutation of [0, pEnd-pStart) */
   PetscInt                     *atlasDof;     /* Describes layout of storage, point --> # of values */
   PetscInt                     *atlasOff;     /* Describes layout of storage, point --> offset into storage */
   PetscInt                      maxDof;       /* Maximum dof on any point */

include/petscdm.h

 PETSC_EXTERN PetscErrorCode DMGetPointSF(DM, PetscSF *);
 PETSC_EXTERN PetscErrorCode DMSetPointSF(DM, PetscSF);
 
+PETSC_EXTERN PetscErrorCode DMGetOutputDM(DM, DM *);
+PETSC_EXTERN PetscErrorCode DMGetOutputSequenceNumber(DM, PetscInt *);
+PETSC_EXTERN PetscErrorCode DMSetOutputSequenceNumber(DM, PetscInt);
+
 PETSC_EXTERN PetscErrorCode DMGetNumFields(DM, PetscInt *);
 PETSC_EXTERN PetscErrorCode DMSetNumFields(DM, PetscInt);
 PETSC_EXTERN PetscErrorCode DMGetField(DM, PetscInt, PetscObject *);

include/petscdmplex.h

 PETSC_EXTERN PetscErrorCode DMPlexComputeCellGeometryFVM(DM, PetscInt, PetscReal *, PetscReal [], PetscReal []);
 
 /* FEM Support */
-PETSC_EXTERN PetscErrorCode DMPlexCreateSection(DM, PetscInt, PetscInt,const PetscInt [],const PetscInt [], PetscInt,const PetscInt [],const IS [], PetscSection *);
+PETSC_EXTERN PetscErrorCode DMPlexCreateSection(DM, PetscInt, PetscInt,const PetscInt [],const PetscInt [], PetscInt,const PetscInt [],const IS [], IS, PetscSection *);
 
 PETSC_EXTERN PetscErrorCode DMPlexComputeCellGeometry(DM, PetscInt, PetscReal *, PetscReal *, PetscReal *, PetscReal *);
 PETSC_EXTERN PetscErrorCode DMPlexVecGetClosure(DM, PetscSection, Vec, PetscInt, PetscInt *, PetscScalar *[]);
   void *user;
 } JacActionCtx;
 
+PETSC_EXTERN PetscErrorCode DMPlexInsertBoundaryValues(DM, Vec);
 PETSC_EXTERN PetscErrorCode DMPlexProjectFunction(DM, PetscFE[], void (**)(const PetscReal [], PetscScalar *, void *), void **, InsertMode, Vec);
 PETSC_EXTERN PetscErrorCode DMPlexProjectFunctionLocal(DM, PetscFE[], void (**)(const PetscReal [], PetscScalar *, void *), void **, InsertMode, Vec);
 PETSC_EXTERN PetscErrorCode DMPlexComputeL2Diff(DM, PetscFE[], void (**)(const PetscReal [], PetscScalar *, void *), void **, Vec, PetscReal *);

include/petscis.h

 PETSC_EXTERN PetscErrorCode ISInvertPermutation(IS,PetscInt,IS*);
 PETSC_EXTERN PetscErrorCode ISView(IS,PetscViewer);
 PETSC_STATIC_INLINE PetscErrorCode ISViewFromOptions(IS A,const char prefix[],const char name[]) {return PetscObjectViewFromOptions((PetscObject)A,prefix,name);}
+PETSC_EXTERN PetscErrorCode ISLoad(IS,PetscViewer);
 PETSC_EXTERN PetscErrorCode ISEqual(IS,IS,PetscBool  *);
 PETSC_EXTERN PetscErrorCode ISSort(IS);
 PETSC_EXTERN PetscErrorCode ISSorted(IS,PetscBool  *);
 PETSC_EXTERN PetscErrorCode PetscSectionSetFieldComponents(PetscSection, PetscInt, PetscInt);
 PETSC_EXTERN PetscErrorCode PetscSectionGetChart(PetscSection, PetscInt *, PetscInt *);
 PETSC_EXTERN PetscErrorCode PetscSectionSetChart(PetscSection, PetscInt, PetscInt);
+PETSC_EXTERN PetscErrorCode PetscSectionGetPermutation(PetscSection, IS *);
+PETSC_EXTERN PetscErrorCode PetscSectionSetPermutation(PetscSection, IS);
 PETSC_EXTERN PetscErrorCode PetscSectionGetDof(PetscSection, PetscInt, PetscInt*);
 PETSC_EXTERN PetscErrorCode PetscSectionSetDof(PetscSection, PetscInt, PetscInt);
 PETSC_EXTERN PetscErrorCode PetscSectionAddDof(PetscSection, PetscInt, PetscInt);
 PETSC_EXTERN PetscErrorCode PetscSectionGetFieldDof(PetscSection, PetscInt, PetscInt, PetscInt*);
 PETSC_EXTERN PetscErrorCode PetscSectionSetFieldDof(PetscSection, PetscInt, PetscInt, PetscInt);
 PETSC_EXTERN PetscErrorCode PetscSectionAddFieldDof(PetscSection, PetscInt, PetscInt, PetscInt);
+PETSC_EXTERN PetscErrorCode PetscSectionHasConstraints(PetscSection, PetscBool *);
 PETSC_EXTERN PetscErrorCode PetscSectionGetConstraintDof(PetscSection, PetscInt, PetscInt*);
 PETSC_EXTERN PetscErrorCode PetscSectionSetConstraintDof(PetscSection, PetscInt, PetscInt);
 PETSC_EXTERN PetscErrorCode PetscSectionAddConstraintDof(PetscSection, PetscInt, PetscInt);

include/petscsys.h

 
 E*/
 typedef enum {PETSC_INT = 0,PETSC_DOUBLE = 1,PETSC_COMPLEX = 2, PETSC_LONG = 3 ,PETSC_SHORT = 4,PETSC_FLOAT = 5,
-              PETSC_CHAR = 6,PETSC_BIT_LOGICAL = 7,PETSC_ENUM = 8,PETSC_BOOL=9, PETSC___FLOAT128 = 10,PETSC_OBJECT = 11, PETSC_FUNCTION = 12} PetscDataType;
+              PETSC_CHAR = 6,PETSC_BIT_LOGICAL = 7,PETSC_ENUM = 8,PETSC_BOOL=9, PETSC___FLOAT128 = 10,PETSC_OBJECT = 11, PETSC_FUNCTION = 12, PETSC_STRING = 12} PetscDataType;
 PETSC_EXTERN const char *const PetscDataTypes[];
 
 #if defined(PETSC_USE_COMPLEX)

include/petscviewer.h

   PETSC_VIEWER_VTK_VTU,
   PETSC_VIEWER_BINARY_MATLAB,
   PETSC_VIEWER_NATIVE,
+  PETSC_VIEWER_HDF5_VIZ,
   PETSC_VIEWER_NOFORMAT
   } PetscViewerFormat;
 PETSC_EXTERN const char *const PetscViewerFormats[];
 PETSC_EXTERN PetscViewer    PETSC_VIEWER_SOCKET_(MPI_Comm);
 PETSC_EXTERN PetscViewer    PETSC_VIEWER_BINARY_(MPI_Comm);
 PETSC_EXTERN PetscViewer    PETSC_VIEWER_MATLAB_(MPI_Comm);
+PETSC_EXTERN PetscViewer    PETSC_VIEWER_HDF5_(MPI_Comm);
 PETSC_EXTERN PetscViewer   PETSC_VIEWER_MATHEMATICA_WORLD_PRIVATE;
 
 #define PETSC_VIEWER_STDERR_SELF  PETSC_VIEWER_STDERR_(PETSC_COMM_SELF)

include/petscviewerhdf5.h

   *b =  (hsize_t)(a);
   PetscFunctionReturn(0);
 }
+PETSC_EXTERN PetscErrorCode PetscDataTypeToHDF5DataType(PetscDataType,hid_t*);
+PETSC_EXTERN PetscErrorCode PetscHDF5DataTypeToPetscDataType(hid_t,PetscDataType*);
 
 #endif  /* defined(PETSC_HAVE_HDF5) */
 
+PETSC_EXTERN PetscErrorCode PetscViewerHDF5ReadAttribute(PetscViewer,const char[],const char[],PetscDataType,void*);
+PETSC_EXTERN PetscErrorCode PetscViewerHDF5WriteAttribute(PetscViewer,const char[],const char[],PetscDataType,const void*);
+PETSC_EXTERN PetscErrorCode PetscViewerHDF5HasAttribute(PetscViewer, const char[], const char[], PetscBool *);
 PETSC_EXTERN PetscErrorCode PetscViewerHDF5WriteSDS(PetscViewer,float *,int,int *,int);
 
 PETSC_EXTERN PetscErrorCode PetscViewerHDF5Open(MPI_Comm,const char[],PetscFileMode,PetscViewer*);
 check: test
 test:
 	-@${OMAKE} PETSC_ARCH=${PETSC_ARCH}  PETSC_DIR=${PETSC_DIR} test_build 2>&1 | tee ./${PETSC_ARCH}/conf/test.log
-	-@if [ "${PETSC_WITH_BATCH}" == "" ]; then printf "Now run make streams NPMAX=<number of expect MPI processes you intend to use>\nto evaluate the computer systems you plan use\n";fi
+	-@if [ "${PETSC_WITH_BATCH}" = "" ]; then \
+          printf "=========================================\n"; \
+          printf "Now to evaluate the computer systems you plan use - do:\n"; \
+          printf "make PETSC_DIR=${PETSC_DIR} PETSC_ARCH=${PETSC_ARCH} streams NPMAX=<number of MPI processes you intend to use>\n"; \
+        fi
 testx:
 	-@${OMAKE} PETSC_ARCH=${PETSC_ARCH}  PETSC_DIR=${PETSC_DIR} testx_build 2>&1 | tee ./${PETSC_ARCH}/conf/testx.log
 test_build:

src/benchmarks/streams/makefile

 
 ALL:
 
-CFLAGS	      = -fopenmp
+CFLAGS	      = 
 FFLAGS	      =
 CPPFLAGS      =
 FPPFLAGS      =

src/benchmarks/streams/process.py

       hosts[size] = fields[1:size+1]
       triads[size] = float(fields[size+5].split()[1])
 
+  if len(hosts) < 2: return
+
   ff = open('scaling.log','a')
   if fileoutput: print 'np  speedup'
   if fileoutput: ff.write('np  speedup\n')

src/dm/dt/interface/dtfe.c

   default:
     SETERRQ1(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_WRONG, "Cannot create reference cell for dimension %d", dim);
   }
+  *refdm = NULL;
   ierr = DMPlexInterpolate(rdm, refdm);CHKERRQ(ierr);
   ierr = DMPlexCopyCoordinates(rdm, *refdm);CHKERRQ(ierr);
   ierr = DMDestroy(&rdm);CHKERRQ(ierr);

src/dm/impls/plex/examples/tests/ex10.c

   ierr = ProcessOptions(&user);CHKERRQ(ierr);
   ierr = DMPlexCreateDoublet(PETSC_COMM_WORLD, user.dim, user.cellSimplex, user.interpolate, user.refinementUniform, user.refinementLimit, &dm);CHKERRQ(ierr);
   ierr = DMSetFromOptions(dm);CHKERRQ(ierr);
-  ierr = DMPlexCreateSection(dm, user.dim, user.numFields, user.numComponents, user.numDof, 0, NULL, NULL, &s);CHKERRQ(ierr);
+  ierr = DMPlexCreateSection(dm, user.dim, user.numFields, user.numComponents, user.numDof, 0, NULL, NULL, NULL, &s);CHKERRQ(ierr);
   ierr = DMSetDefaultSection(dm, s);CHKERRQ(ierr);
   ierr = PetscSectionDestroy(&s);CHKERRQ(ierr);
   ierr = TestReordering(dm, &user);CHKERRQ(ierr);

src/dm/impls/plex/examples/tests/ex4.c

 */
 PetscErrorCode CreateSimplexHybrid_2D(MPI_Comm comm, PetscInt testNum, DM *dm)
 {
-  DM             idm, hdm;
+  DM             idm = NULL, hdm = NULL;
   DMLabel        faultLabel, hybridLabel;
   PetscInt       p;
   PetscMPIInt    rank;
 #define __FUNCT__ "CreateTensorProductHybrid_2D"
 PetscErrorCode CreateTensorProductHybrid_2D(MPI_Comm comm, PetscInt testNum, DM *dm)
 {
-  DM             idm, hdm;
+  DM             idm = NULL, hdm = NULL;
   DMLabel        faultLabel, hybridLabel;
   PetscInt       p;
   PetscMPIInt    rank;
 */
 PetscErrorCode CreateSimplex_3D(MPI_Comm comm, PetscInt testNum, DM *dm)
 {
-  DM             idm;
+  DM             idm = NULL;
   PetscInt       depth = 3;
   PetscMPIInt    rank;
   PetscErrorCode ierr;
 */
 PetscErrorCode CreateSimplexHybrid_3D(MPI_Comm comm, PetscInt testNum, DM *dm)
 {
-  DM             idm, hdm;
+  DM             idm = NULL, hdm = NULL;
   DMLabel        faultLabel, hybridLabel;
   PetscInt       p;
   PetscMPIInt    rank;
 #define __FUNCT__ "CreateTensorProduct_3D"
 PetscErrorCode CreateTensorProduct_3D(MPI_Comm comm, PetscInt testNum, DM *dm)
 {
-  DM             idm;
+  DM             idm = NULL;
   PetscMPIInt    rank;
   PetscErrorCode ierr;
 
 #define __FUNCT__ "CreateTensorProductHybrid_3D"
 PetscErrorCode CreateTensorProductHybrid_3D(MPI_Comm comm, PetscInt testNum, DM *dm)
 {
-  DM             idm, hdm;
+  DM             idm = NULL, hdm = NULL;
   DMLabel        faultLabel;
   PetscInt       p;
   PetscMPIInt    rank;

src/dm/impls/plex/examples/tests/ex5.c

 #define __FUNCT__ "CreateSimplex_2D"
 PetscErrorCode CreateSimplex_2D(MPI_Comm comm, PetscInt testNum, DM *dm)
 {
-  DM             idm;
+  DM             idm = NULL;
   PetscInt       p;
   PetscMPIInt    rank;
   PetscErrorCode ierr;
 #define __FUNCT__ "CreateQuad_2D"
 PetscErrorCode CreateQuad_2D(MPI_Comm comm, PetscInt testNum, DM *dm)
 {
-  DM             idm;
+  DM             idm = NULL;
   PetscInt       p;
   PetscMPIInt    rank;
   PetscErrorCode ierr;
 #define __FUNCT__ "CreateHex_3D"
 PetscErrorCode CreateHex_3D(MPI_Comm comm, PetscInt testNum, DM *dm)
 {
-  DM             idm;
+  DM             idm   = NULL;
   PetscInt       depth = 3, p;
   PetscMPIInt    rank;
   PetscErrorCode ierr;

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

 
     ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
     if (interpolate) {
-      DM idm;
+      DM idm = NULL;
 
       ierr = DMPlexInterpolate(dm, &idm);CHKERRQ(ierr);
       ierr = PetscObjectSetName((PetscObject) idm, "triangle");CHKERRQ(ierr);
 
     ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
     if (interpolate) {
-      DM idm;
+      DM idm = NULL;
 
       ierr = DMPlexInterpolate(dm, &idm);CHKERRQ(ierr);
       ierr = PetscObjectSetName((PetscObject) idm, "quadrilateral");CHKERRQ(ierr);
 
     ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
     if (interpolate) {
-      DM idm;
+      DM idm = NULL;
 
       ierr = DMPlexInterpolate(dm, &idm);CHKERRQ(ierr);
       ierr = PetscObjectSetName((PetscObject) idm, "tetrahedron");CHKERRQ(ierr);
 
     ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
     if (interpolate) {
-      DM idm;
+      DM idm = NULL;
 
       ierr = DMPlexInterpolate(dm, &idm);CHKERRQ(ierr);
       ierr = PetscObjectSetName((PetscObject) idm, "hexahedron");CHKERRQ(ierr);

src/dm/impls/plex/examples/tests/ex9.c

     *newdm = rdm;
   }
   if (user->interpolate) {
-    DM idm;
+    DM idm = NULL;
     const char *name;
 
     ierr = DMPlexInterpolate(*newdm, &idm);CHKERRQ(ierr);
     ierr = PetscLogEventRegister("VecClosure", PETSC_OBJECT_CLASSID, &event);CHKERRQ(ierr);
   }
   ierr = PetscLogStagePush(stage);CHKERRQ(ierr);
-  ierr = DMPlexCreateSection(dm, user->dim, user->numFields, user->numComponents, user->numDof, 0, NULL, NULL, &s);CHKERRQ(ierr);
+  ierr = DMPlexCreateSection(dm, user->dim, user->numFields, user->numComponents, user->numDof, 0, NULL, NULL, NULL, &s);CHKERRQ(ierr);
   ierr = DMSetDefaultSection(dm, s);CHKERRQ(ierr);
   if (useIndex) {ierr = DMPlexCreateClosureIndex(dm, s);CHKERRQ(ierr);}
   ierr = PetscSectionDestroy(&s);CHKERRQ(ierr);

src/dm/impls/plex/examples/tutorials/ex1.c

   bcField[0] = 0;
   ierr = DMPlexGetStratumIS(dm, "marker", 1, &bcPointIS[0]);CHKERRQ(ierr);
   /* Create a PetscSection with this data layout */
-  ierr = DMPlexCreateSection(dm, dim, numFields, numComp, numDof, numBC, bcField, bcPointIS, &section);CHKERRQ(ierr);
+  ierr = DMPlexCreateSection(dm, dim, numFields, numComp, numDof, numBC, bcField, bcPointIS, NULL, &section);CHKERRQ(ierr);
   ierr = ISDestroy(&bcPointIS[0]);CHKERRQ(ierr);
   /* Name the Field variables */
   ierr = PetscSectionSetFieldName(section, 0, "u");CHKERRQ(ierr);

src/dm/impls/plex/examples/tutorials/ex1f90.F

       pBcPointIS => bcPointIS
 !     Create a PetscSection with this data layout
       call DMPlexCreateSection(dm, dim, numFields, pNumComp,
-     $     pNumDof, numBC, pBcField, pBcPointIS, section,
-     $     ierr)
+     $     pNumDof, numBC, pBcField, pBcPointIS, PETSC_NULL_OBJECT,
+     $     section, ierr)
       CHKERRQ(ierr)
       call ISDestroy(bcPointIS(1), ierr)
       CHKERRQ(ierr)

src/dm/impls/plex/examples/tutorials/ex2.c

+static char help[] = "Read in a mesh and test whether it is valid\n\n";
+
+#include <petscdmplex.h>
+#if defined(PETSC_HAVE_CGNS)
+#include <cgnslib.h>
+#endif
+#if defined(PETSC_HAVE_EXODUSII)
+#include <exodusII.h>
+#endif
+
+typedef struct {
+  PetscBool interpolate;                  /* Generate intermediate mesh elements */
+  char      filename[PETSC_MAX_PATH_LEN]; /* Mesh filename */
+} AppCtx;
+
+#undef __FUNCT__
+#define __FUNCT__ "ProcessOptions"
+static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
+{
+  PetscErrorCode ierr;
+
+  PetscFunctionBeginUser;
+  options->interpolate = PETSC_FALSE;
+  options->filename[0] = '\0';
+
+  ierr = PetscOptionsBegin(comm, "", "Meshing Problem Options", "DMPLEX");CHKERRQ(ierr);
+  ierr = PetscOptionsBool("-interpolate", "Generate intermediate mesh elements", "ex1.c", options->interpolate, &options->interpolate, 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);
+};
+
+#undef __FUNCT__
+#define __FUNCT__ "CreateMesh"
+static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
+{
+  PetscBool      interpolate = user->interpolate;
+  const char    *filename    = user->filename;
+  const char    *extGmsh     = ".msh";
+  const char    *extCGNS     = ".cgns";
+  const char    *extExodus   = ".exo";
+  size_t         len;
+  PetscBool      isGmsh, isCGNS, isExodus;
+  PetscMPIInt    rank;
+  PetscErrorCode ierr;
+
+  PetscFunctionBeginUser;
+  ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
+  ierr = PetscStrlen(filename, &len);CHKERRQ(ierr);
+  if (!len) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Must provide input file using -filename");
+  ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extGmsh,   4, &isGmsh);CHKERRQ(ierr);
+  ierr = PetscStrncmp(&filename[PetscMax(0,len-5)], extCGNS,   5, &isCGNS);CHKERRQ(ierr);
+  ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extExodus, 4, &isExodus);CHKERRQ(ierr);
+  if (isGmsh) {
+    PetscViewer viewer;
+
+    ierr = PetscViewerCreate(comm, &viewer);CHKERRQ(ierr);
+    ierr = PetscViewerSetType(viewer, PETSCVIEWERASCII);CHKERRQ(ierr);
+    ierr = PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRQ(ierr);
+    ierr = PetscViewerFileSetName(viewer, filename);CHKERRQ(ierr);
+    ierr = DMPlexCreateGmsh(comm, viewer, interpolate, dm);CHKERRQ(ierr);
+    ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
+  } else if (isCGNS) {
+#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 (isExodus) {
+#if defined(PETSC_HAVE_EXODUSII)
+    int   CPU_word_size = 0, IO_word_size = 0, exoid;
+    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);
+    }
+    ierr = DMPlexCreateExodus(comm, exoid, interpolate, dm);CHKERRQ(ierr);
+    if (!rank) {ierr = ex_close(exoid);CHKERRQ(ierr);}
+#else
+    SETERRQ(comm, PETSC_ERR_SUP, "Loading meshes requires Exodus support. Reconfigure using --download-exodusii");
+#endif
+  } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot load file %s: unrecognized extension", filename);
+  ierr = PetscObjectSetName((PetscObject) *dm, "Mesh");CHKERRQ(ierr);
+  ierr = DMSetFromOptions(*dm);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "CheckMeshTopology"
+static PetscErrorCode CheckMeshTopology(DM dm)
+{
+  PetscInt       dim, coneSize, cStart;
+  PetscBool      isSimplex;
+  PetscErrorCode ierr;
+
+  PetscFunctionBeginUser;
+  ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
+  ierr = DMPlexGetHeightStratum(dm, 0, &cStart, NULL);CHKERRQ(ierr);
+  ierr = DMPlexGetConeSize(dm, cStart, &coneSize);CHKERRQ(ierr);
+  isSimplex = coneSize == dim+1 ? PETSC_TRUE : PETSC_FALSE;
+  ierr = DMPlexCheckSymmetry(dm);CHKERRQ(ierr);
+  ierr = DMPlexCheckSkeleton(dm, isSimplex, 0);CHKERRQ(ierr);
+  ierr = DMPlexCheckFaces(dm, isSimplex, 0);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "CheckMeshGeometry"
+static PetscErrorCode CheckMeshGeometry(DM dm)
+{
+  PetscInt       dim, coneSize, cStart, cEnd, c;
+  PetscReal     *v0, *J, *invJ, detJ;
+  PetscBool      isSimplex;
+  PetscErrorCode ierr;
+
+  PetscFunctionBeginUser;
+  ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
+  ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
+  ierr = DMPlexGetConeSize(dm, cStart, &coneSize);CHKERRQ(ierr);
+  isSimplex = coneSize == dim+1 ? PETSC_TRUE : PETSC_FALSE;
+  ierr = PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr);
+  for (c = cStart; c < cEnd; ++c) {
+    ierr = DMPlexComputeCellGeometry(dm, c, v0, J, invJ, &detJ);CHKERRQ(ierr);
+    if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for cell %d", detJ, c);
+  }
+  ierr = PetscFree3(v0,J,invJ);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "main"
+int main(int argc, char **argv)
+{
+  DM             dm;
+  AppCtx         user;
+  PetscErrorCode ierr;
+
+  ierr = PetscInitialize(&argc, &argv, NULL, help);CHKERRQ(ierr);
+  ierr = ProcessOptions(PETSC_COMM_WORLD, &user);CHKERRQ(ierr);
+  ierr = CreateMesh(PETSC_COMM_WORLD, &user, &dm);CHKERRQ(ierr);
+  ierr = CheckMeshTopology(dm);CHKERRQ(ierr);
+  ierr = CheckMeshGeometry(dm);CHKERRQ(ierr);
+  ierr = DMDestroy(&dm);CHKERRQ(ierr);
+  ierr = PetscFinalize();
+  return 0;
+}

src/dm/impls/plex/examples/tutorials/ex3f90.F

+! setting up 3-D DMPlex using DMPlexCreateFromDAG(), DMPlexInterpolate() and
+! DMPlexComputeCellGeometryFVM()
+! Contributed by Adrian Croucher <a.croucher@auckland.ac.nz>
+      program main
+      implicit none
+#include <finclude/petsc.h90>
+      DM :: dm, dmi
+      PetscInt, parameter :: dim = 3
+
+      PetscInt :: depth = 1
+      PetscErrorCode :: ierr
+      PetscInt, dimension(2) :: numPoints
+      PetscInt, dimension(14) :: coneSize
+      PetscInt, dimension(16) :: cones
+      PetscInt, dimension(16) ::                                        &
+     & coneOrientations
+      PetscScalar, dimension(36) :: vertexCoords
+      PetscReal ::  vol = 0.
+      PetscReal, target, dimension(dim)  ::                             &
+     & centroid
+      PetscReal, target, dimension(dim)  ::                             &
+     & normal
+      PetscReal, pointer :: pcentroid(:)
+      PetscReal, pointer :: pnormal(:)
+
+      PetscInt :: i
+
+      numPoints = [12, 2]
+      coneSize  = [8, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
+      cones = [2,5,4,3,6,7,8,9,  3,4,11,10,7,12,13,8]
+      coneOrientations = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,0]
+      vertexCoords = [                                                  &
+     &  -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 ]
+      pcentroid => centroid
+      pnormal => normal
+
+
+      call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
+      CHKERRQ(ierr)
+
+      call DMPlexCreate(PETSC_COMM_WORLD, dm, ierr)
+      CHKERRQ(ierr)
+      call PetscObjectSetName(dm, "testplex", ierr)
+      CHKERRQ(ierr)
+      call DMPlexSetDimension(dm, dim, ierr)
+      CHKERRQ(ierr)
+
+      call DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones,   &
+     &  coneOrientations, vertexCoords, ierr)
+      CHKERRQ(ierr)
+
+      call DMPlexInterpolate(dm, dmi, ierr)
+      CHKERRQ(ierr)
+      call DMPlexCopyCoordinates(dm, dmi, ierr)
+      CHKERRQ(ierr)
+      call DMDestroy(dm, ierr)
+      CHKERRQ(ierr)
+      dm = dmi
+
+      call DMView(dm, PETSC_VIEWER_STDOUT_WORLD, ierr)
+      CHKERRQ(ierr)
+
+      do i = 0, 1
+        call DMPlexComputeCellGeometryFVM(dm, i, vol, pcentroid, pnormal&
+     &  , ierr)
+        CHKERRQ(ierr)
+        write(*, '(a, i2, a, f8.4, a, 3(f8.4, 1x))'),                   &
+     &           'cell: ', i, ' volume: ', vol, ' centroid: ',          &
+     &           pcentroid(1), pcentroid(2), pcentroid(3)
+      end do
+
+      call DMDestroy(dm, ierr)
+      CHKERRQ(ierr)
+      call PetscFinalize(ierr)
+      CHKERRQ(ierr)
+      end program main

src/dm/impls/plex/examples/tutorials/output/ex3f90_0.out

+DM_0x84000000_0 in 3 dimensions:
+  0-cells: 12
+  1-cells: 20
+  2-cells: 11
+  3-cells: 2
+Labels:
+  depth: 4 strata of sizes (12, 20, 11, 2)
+cell:  0 volume:   0.5000 centroid:  -0.2500   0.5000   0.5000
+cell:  1 volume:   0.5000 centroid:   0.2500   0.5000   0.5000

src/dm/impls/plex/f90-custom/zplexf90.c

 #define dmplexrestoremeet_              DMPLEXRESTOREMEET
 #define dmplexcreatesection_            DMPLEXCREATESECTION
 #define dmplexcomputecellgeometry_      DMPLEXCOMPUTECELLGEOMETRY
+#define dmplexcomputecellgeometryfvm_   DMPLEXCOMPUTECELLGEOMETRYFVM
 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) && !defined(FORTRANDOUBLEUNDERSCORE)
 #define dmplexgetcone_                  dmplexgetcone
 #define dmplexrestorecone_              dmplexrestorecone
 #define dmplexrestoremeet_              dmplexrestoremeet
 #define dmplexcreatesection_            dmplexcreatesection
 #define dmplexcomputecellgeometry_      dmplexcomputecellgeometry
+#define dmplexcomputecellgeometryfvm_   dmplexcomputecellgeometryfvm
 #endif
 
 /* Definitions of Fortran Wrapper routines */
   *ierr = F90Array1dDestroy(cptr, PETSC_INT PETSC_F90_2PTR_PARAM(cptrd));if (*ierr) return;
 }
 
-PETSC_EXTERN void PETSC_STDCALL dmplexcreatesection_(DM *dm, PetscInt *dim, PetscInt *numFields, F90Array1d *ptrC, F90Array1d *ptrD, PetscInt *numBC, F90Array1d *ptrF, F90Array1d *ptrP, PetscSection *section, int *ierr PETSC_F90_2PTR_PROTO(ptrCd) PETSC_F90_2PTR_PROTO(ptrDd) PETSC_F90_2PTR_PROTO(ptrFd) PETSC_F90_2PTR_PROTO(ptrPd))
+PETSC_EXTERN void PETSC_STDCALL dmplexcreatesection_(DM *dm, PetscInt *dim, PetscInt *numFields, F90Array1d *ptrC, F90Array1d *ptrD, PetscInt *numBC, F90Array1d *ptrF, F90Array1d *ptrP, IS *perm, PetscSection *section, int *ierr PETSC_F90_2PTR_PROTO(ptrCd) PETSC_F90_2PTR_PROTO(ptrDd) PETSC_F90_2PTR_PROTO(ptrFd) PETSC_F90_2PTR_PROTO(ptrPd))
 {
   PetscInt *numComp;
   PetscInt *numDof;
   PetscInt *bcField;
   IS       *bcPoints;
 
+  CHKFORTRANNULLOBJECT(perm);
   *ierr = F90Array1dAccess(ptrC, PETSC_INT, (void**) &numComp PETSC_F90_2PTR_PARAM(ptrCd));if (*ierr) return;
   *ierr = F90Array1dAccess(ptrD, PETSC_INT, (void**) &numDof PETSC_F90_2PTR_PARAM(ptrDd));if (*ierr) return;
   *ierr = F90Array1dAccess(ptrF, PETSC_INT, (void**) &bcField PETSC_F90_2PTR_PARAM(ptrFd));if (*ierr) return;
   *ierr = F90Array1dAccess(ptrP, PETSC_FORTRANADDR, (void**) &bcPoints PETSC_F90_2PTR_PARAM(ptrPd));if (*ierr) return;
-  *ierr = DMPlexCreateSection(*dm, *dim, *numFields, numComp, numDof, *numBC, bcField, bcPoints, section);
+  *ierr = DMPlexCreateSection(*dm, *dim, *numFields, numComp, numDof, *numBC, bcField, bcPoints, *perm, section);
 }
 
 PETSC_EXTERN void PETSC_STDCALL dmplexcomputecellgeometry_(DM *dm, PetscInt *cell, F90Array1d *ptrV, F90Array1d *ptrJ, F90Array1d *ptrIJ, PetscReal *detJ, int *ierr PETSC_F90_2PTR_PROTO(ptrVd) PETSC_F90_2PTR_PROTO(ptrJd) PETSC_F90_2PTR_PROTO(ptrIJd))
   *ierr = DMPlexComputeCellGeometry(*dm, *cell, v0, J, invJ, detJ);
 }
 
+PETSC_EXTERN void PETSC_STDCALL dmplexcomputecellgeometryfvm_(DM *dm, PetscInt *cell, PetscReal *vol, F90Array1d *ptrCentroid, F90Array1d *ptrNormal, int *ierr PETSC_F90_2PTR_PROTO(ptrCentroidd) PETSC_F90_2PTR_PROTO(ptrNormald))
+{
+  PetscReal *centroid;
+  PetscReal *normal;
+
+  *ierr = F90Array1dAccess(ptrCentroid, PETSC_REAL, (void**) &centroid PETSC_F90_2PTR_PARAM(ptrCentroidd));if (*ierr) return;
+  *ierr = F90Array1dAccess(ptrNormal,   PETSC_REAL, (void**) &normal   PETSC_F90_2PTR_PARAM(ptrNormald));if (*ierr) return;
+  *ierr = DMPlexComputeCellGeometryFVM(*dm, *cell, vol, centroid, normal);
+}
+

src/dm/impls/plex/plex.c

 #include <../src/sys/utils/hash.h>
 #include <petsc-private/isimpl.h>
 #include <petscsf.h>
+#include <petscviewerhdf5.h>
 
 /* Logging support */
 PetscLogEvent DMPLEX_Interpolate, DMPLEX_Partition, DMPLEX_Distribute, DMPLEX_DistributeCones, DMPLEX_DistributeLabels, DMPLEX_DistributeSF, DMPLEX_DistributeField, DMPLEX_DistributeData, DMPLEX_Stratify, DMPLEX_Preallocate, DMPLEX_ResidualFEM, DMPLEX_JacobianFEM;
 
 PETSC_EXTERN PetscErrorCode VecView_Seq(Vec, PetscViewer);
 PETSC_EXTERN PetscErrorCode VecView_MPI(Vec, PetscViewer);
+PETSC_EXTERN PetscErrorCode VecLoad_Default(Vec, PetscViewer);
+PETSC_EXTERN PetscErrorCode DMTSGetTimeStepNumber(DM,PetscInt*);
+
+#undef __FUNCT__
+#define __FUNCT__ "GetFieldType_Static"
+static PetscErrorCode GetFieldType_Static(DM dm, PetscSection section, PetscInt field, PetscInt *sStart, PetscInt *sEnd, PetscViewerVTKFieldType *ft)
+{
+  PetscInt       dim, pStart, pEnd, vStart, vEnd, cStart, cEnd, vdof = 0, cdof = 0;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  *ft  = PETSC_VTK_POINT_FIELD;
+  ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
+  ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
+  ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
+  ierr = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr);
+  if (field >= 0) {
+    if ((vStart >= pStart) && (vStart < pEnd)) {ierr = PetscSectionGetFieldDof(section, vStart, field, &vdof);CHKERRQ(ierr);}
+    if ((cStart >= pStart) && (cStart < pEnd)) {ierr = PetscSectionGetFieldDof(section, cStart, field, &cdof);CHKERRQ(ierr);}
+  } else {
+    if ((vStart >= pStart) && (vStart < pEnd)) {ierr = PetscSectionGetDof(section, vStart, &vdof);CHKERRQ(ierr);}
+    if ((cStart >= pStart) && (cStart < pEnd)) {ierr = PetscSectionGetDof(section, cStart, &cdof);CHKERRQ(ierr);}
+  }
+  if (vdof) {
+    *sStart = vStart;
+    *sEnd   = vEnd;
+    if (vdof == dim) *ft = PETSC_VTK_POINT_VECTOR_FIELD;
+    else             *ft = PETSC_VTK_POINT_FIELD;
+  } else if (cdof) {
+    *sStart = cStart;
+    *sEnd   = cEnd;
+    if (cdof == dim) *ft = PETSC_VTK_CELL_VECTOR_FIELD;
+    else             *ft = PETSC_VTK_CELL_FIELD;
+  } else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Could not classify input Vec for VTK");
+  PetscFunctionReturn(0);
+}
+
+#if defined(PETSC_HAVE_HDF5)
+#undef __FUNCT__
+#define __FUNCT__ "GetField_Static"
+static PetscErrorCode GetField_Static(DM dm, PetscSection section, PetscSection sectionGlobal, Vec v, PetscInt field, PetscInt pStart, PetscInt pEnd, IS *is, Vec *subv)
+{
+  PetscInt      *subIndices;
+  PetscInt       Nc, subSize = 0, subOff = 0, p;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  ierr = PetscSectionGetFieldComponents(section, field, &Nc);CHKERRQ(ierr);
+  for (p = pStart; p < pEnd; ++p) {
+    PetscInt gdof, fdof = 0;
+
+    ierr = PetscSectionGetDof(sectionGlobal, p, &gdof);CHKERRQ(ierr);
+    if (gdof > 0) {ierr = PetscSectionGetFieldDof(section, p, field, &fdof);CHKERRQ(ierr);}
+    subSize += fdof;
+  }
+  ierr = PetscMalloc1(subSize, &subIndices);CHKERRQ(ierr);
+  for (p = pStart; p < pEnd; ++p) {
+    PetscInt gdof, goff;
+
+    ierr = PetscSectionGetDof(sectionGlobal, p, &gdof);CHKERRQ(ierr);
+    if (gdof > 0) {
+      PetscInt fdof, fc, f2, poff = 0;
+
+      ierr = PetscSectionGetOffset(sectionGlobal, p, &goff);CHKERRQ(ierr);
+      /* Can get rid of this loop by storing field information in the global section */
+      for (f2 = 0; f2 < field; ++f2) {
+        ierr  = PetscSectionGetFieldDof(section, p, f2, &fdof);CHKERRQ(ierr);
+        poff += fdof;
+      }
+      ierr = PetscSectionGetFieldDof(section, p, field, &fdof);CHKERRQ(ierr);
+      for (fc = 0; fc < fdof; ++fc, ++subOff) subIndices[subOff] = goff+poff+fc;
+    }
+  }
+  ierr = ISCreateGeneral(PetscObjectComm((PetscObject) dm), subSize, subIndices, PETSC_OWN_POINTER, is);CHKERRQ(ierr);
+  ierr = VecGetSubVector(v, *is, subv);CHKERRQ(ierr);
+  ierr = VecSetBlockSize(*subv, Nc);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "RestoreField_Static"
+static PetscErrorCode RestoreField_Static(DM dm, PetscSection section, PetscSection sectionGlobal, Vec v, PetscInt field, PetscInt pStart, PetscInt pEnd, IS *is, Vec *subv)
+{
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  ierr = VecRestoreSubVector(v, *is, subv);CHKERRQ(ierr);
+  ierr = ISDestroy(is);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "DMSequenceView_HDF5"
+static PetscErrorCode DMSequenceView_HDF5(DM dm, const char *seqname, PetscInt seqnum, PetscScalar value, PetscViewer viewer)
+{
+  Vec            stamp;
+  PetscMPIInt    rank;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) viewer), &rank);CHKERRQ(ierr);
+  ierr = VecCreateMPI(PetscObjectComm((PetscObject) viewer), rank ? 0 : 1, 1, &stamp);CHKERRQ(ierr);
+  ierr = VecSetBlockSize(stamp, 1);CHKERRQ(ierr);
+  ierr = PetscObjectSetName((PetscObject) stamp, seqname);CHKERRQ(ierr);
+  if (!rank) {
+    PetscReal timeScale;
+    PetscBool istime;
+
+    ierr = PetscStrncmp(seqname, "time", 5, &istime);CHKERRQ(ierr);
+    if (istime) {ierr = DMPlexGetScale(dm, PETSC_UNIT_TIME, &timeScale);CHKERRQ(ierr); value *= timeScale;}
+    ierr = VecSetValue(stamp, 0, value, INSERT_VALUES);CHKERRQ(ierr);
+  }
+  ierr = VecAssemblyBegin(stamp);CHKERRQ(ierr);
+  ierr = VecAssemblyEnd(stamp);CHKERRQ(ierr);
+  ierr = PetscViewerHDF5PushGroup(viewer, "/");CHKERRQ(ierr);
+  ierr = PetscViewerHDF5SetTimestep(viewer, seqnum);CHKERRQ(ierr);
+  ierr = VecView(stamp, viewer);CHKERRQ(ierr);
+  ierr = PetscViewerHDF5PopGroup(viewer);CHKERRQ(ierr);
+  ierr = VecDestroy(&stamp);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "VecView_Plex_Local_HDF5"
+static PetscErrorCode VecView_Plex_Local_HDF5(Vec v, PetscViewer viewer)
+{
+  DM                      dm;
+  DM                      dmBC;
+  PetscSection            section, sectionGlobal;
+  Vec                     gv;
+  const char             *name;
+  PetscViewerVTKFieldType ft;
+  PetscViewerFormat       format;
+  PetscInt                seqnum;
+  PetscBool               isseq;
+  PetscErrorCode          ierr;
+
+  PetscFunctionBegin;
+  ierr = PetscObjectTypeCompare((PetscObject) v, VECSEQ, &isseq);CHKERRQ(ierr);
+  ierr = VecGetDM(v, &dm);CHKERRQ(ierr);
+  ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
+  ierr = DMGetOutputSequenceNumber(dm, &seqnum);CHKERRQ(ierr);
+  ierr = PetscViewerHDF5SetTimestep(viewer, seqnum);CHKERRQ(ierr);
+  ierr = DMSequenceView_HDF5(dm, "time", seqnum, (PetscScalar) seqnum, viewer);CHKERRQ(ierr);
+  ierr = PetscViewerGetFormat(viewer, &format);CHKERRQ(ierr);
+  ierr = DMGetOutputDM(dm, &dmBC);CHKERRQ(ierr);
+  ierr = DMGetDefaultGlobalSection(dmBC, &sectionGlobal);CHKERRQ(ierr);
+  ierr = DMGetGlobalVector(dmBC, &gv);CHKERRQ(ierr);
+  ierr = PetscObjectGetName((PetscObject) v, &name);CHKERRQ(ierr);
+  ierr = PetscObjectSetName((PetscObject) gv, name);CHKERRQ(ierr);
+  ierr = DMLocalToGlobalBegin(dmBC, v, INSERT_VALUES, gv);CHKERRQ(ierr);
+  ierr = DMLocalToGlobalEnd(dmBC, v, INSERT_VALUES, gv);CHKERRQ(ierr);
+  ierr = PetscObjectTypeCompare((PetscObject) gv, VECSEQ, &isseq);CHKERRQ(ierr);
+  if (isseq) {ierr = VecView_Seq(gv, viewer);CHKERRQ(ierr);}
+  else       {ierr = VecView_MPI(gv, viewer);CHKERRQ(ierr);}
+  if (format == PETSC_VIEWER_HDF5_VIZ) {
+    /* Output visualization representation */
+    PetscInt numFields, f;
+
+    ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr);
+    for (f = 0; f < numFields; ++f) {
+      Vec         subv;
+      IS          is;
+      const char *fname, *fgroup;
+      char        group[PETSC_MAX_PATH_LEN];
+      PetscInt    pStart, pEnd;
+      PetscBool   flag;
+
+      ierr = GetFieldType_Static(dm, section, f, &pStart, &pEnd, &ft);CHKERRQ(ierr);
+      fgroup = (ft == PETSC_VTK_POINT_VECTOR_FIELD) || (ft == PETSC_VTK_POINT_FIELD) ? "/vertex_fields" : "/cell_fields";
+      ierr = PetscSectionGetFieldName(section, f, &fname);CHKERRQ(ierr);
+      ierr = PetscViewerHDF5PushGroup(viewer, fgroup);CHKERRQ(ierr);
+      ierr = GetField_Static(dmBC, section, sectionGlobal, gv, f, pStart, pEnd, &is, &subv);CHKERRQ(ierr);
+      ierr = PetscObjectSetName((PetscObject) subv, fname);CHKERRQ(ierr);
+      if (isseq) {ierr = VecView_Seq(subv, viewer);CHKERRQ(ierr);}
+      else       {ierr = VecView_MPI(subv, viewer);CHKERRQ(ierr);}
+      ierr = RestoreField_Static(dmBC, section, sectionGlobal, gv, f, pStart, pEnd, &is, &subv);CHKERRQ(ierr);
+      ierr = PetscViewerHDF5PopGroup(viewer);CHKERRQ(ierr);
+      ierr = PetscSNPrintf(group, PETSC_MAX_PATH_LEN, "%s/%s", fgroup, fname);CHKERRQ(ierr);
+      ierr = PetscViewerHDF5HasAttribute(viewer, group, "vector_field_type", &flag);CHKERRQ(ierr);
+      if (!flag) {
+        if ((ft == PETSC_VTK_POINT_VECTOR_FIELD) || (ft == PETSC_VTK_CELL_VECTOR_FIELD)) {
+          ierr = PetscViewerHDF5WriteAttribute(viewer, group, "vector_field_type", PETSC_STRING, "vector");CHKERRQ(ierr);
+        } else {
+          ierr = PetscViewerHDF5WriteAttribute(viewer, group, "vector_field_type", PETSC_STRING, "scalar");CHKERRQ(ierr);
+        }
+      }
+    }
+  }
+  ierr = DMRestoreGlobalVector(dmBC, &gv);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+#endif
 
 #undef __FUNCT__
 #define __FUNCT__ "VecView_Plex_Local"
 PetscErrorCode VecView_Plex_Local(Vec v, PetscViewer viewer)
 {
   DM             dm;
-  PetscBool      isvtk;
+  PetscBool      isvtk, ishdf5, isseq;
   PetscErrorCode ierr;
 
   PetscFunctionBegin;
   ierr = VecGetDM(v, &dm);CHKERRQ(ierr);
   if (!dm) SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_WRONG, "Vector not generated from a DM");
-  ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERVTK, &isvtk);CHKERRQ(ierr);
+  ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERVTK,  &isvtk);CHKERRQ(ierr);
+  ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);CHKERRQ(ierr);
+  ierr = PetscObjectTypeCompare((PetscObject) v, VECSEQ, &isseq);CHKERRQ(ierr);
+  if (isvtk || ishdf5) {ierr = DMPlexInsertBoundaryValues(dm, v);CHKERRQ(ierr);}
   if (isvtk) {
-    PetscViewerVTKFieldType ft = PETSC_VTK_POINT_FIELD;
     PetscSection            section;
-    PetscInt                dim, pStart, pEnd, cStart, fStart, vStart, cdof = 0, fdof = 0, vdof = 0;
+    PetscViewerVTKFieldType ft;
+    PetscInt                pStart, pEnd;
 
-    ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
     ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
-    ierr = DMPlexGetHeightStratum(dm, 0, &cStart, NULL);CHKERRQ(ierr);
-    ierr = DMPlexGetHeightStratum(dm, 1, &fStart, NULL);CHKERRQ(ierr);
-    ierr = DMPlexGetDepthStratum(dm, 0, &vStart, NULL);CHKERRQ(ierr);
-    ierr = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr);
-    /* Assumes that numer of dofs per point of each stratum is constant, natural for VTK */
-    if ((cStart >= pStart) && (cStart < pEnd)) {ierr = PetscSectionGetDof(section, cStart, &cdof);CHKERRQ(ierr);}
-    if ((fStart >= pStart) && (fStart < pEnd)) {ierr = PetscSectionGetDof(section, fStart, &fdof);CHKERRQ(ierr);}
-    if ((vStart >= pStart) && (vStart < pEnd)) {ierr = PetscSectionGetDof(section, vStart, &vdof);CHKERRQ(ierr);}
-    if (cdof && fdof && vdof) { /* Actually Q2 or some such, but visualize as Q1 */
-      ft = (cdof == dim) ? PETSC_VTK_POINT_VECTOR_FIELD : PETSC_VTK_POINT_FIELD;
-    } else if (cdof && vdof) {
-      SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_SUP,"No support for viewing mixed space with dofs at both vertices and cells");
-    } else if (cdof) {
-      /* TODO: This assumption should be removed when there is a way of identifying whether a space is conceptually a
-       * vector or just happens to have the same number of dofs as the dimension. */
-      if (cdof == dim) {
-        ft = PETSC_VTK_CELL_VECTOR_FIELD;
-      } else {
-        ft = PETSC_VTK_CELL_FIELD;
-      }
-    } else if (vdof) {
-      if (vdof == dim) {
-        ft = PETSC_VTK_POINT_VECTOR_FIELD;
-      } else {
-        ft = PETSC_VTK_POINT_FIELD;
-      }
-    } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Could not classify input Vec for VTK");
-
+    ierr = GetFieldType_Static(dm, section, PETSC_DETERMINE, &pStart, &pEnd, &ft);CHKERRQ(ierr);
     ierr = PetscObjectReference((PetscObject) dm);CHKERRQ(ierr); /* viewer drops reference */
     ierr = PetscObjectReference((PetscObject) v);CHKERRQ(ierr);  /* viewer drops reference */
     ierr = PetscViewerVTKAddField(viewer, (PetscObject) dm, DMPlexVTKWriteAll, ft, (PetscObject) v);CHKERRQ(ierr);
+  } else if (ishdf5) {
+#if defined(PETSC_HAVE_HDF5)
+    ierr = VecView_Plex_Local_HDF5(v, viewer);CHKERRQ(ierr);
+#else
+    SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
+#endif
   } else {
-    PetscBool isseq;
-
-    ierr = PetscObjectTypeCompare((PetscObject) v, VECSEQ, &isseq);CHKERRQ(ierr);
-    if (isseq) {
-      ierr = VecView_Seq(v, viewer);CHKERRQ(ierr);
-    } else {
-      ierr = VecView_MPI(v, viewer);CHKERRQ(ierr);
-    }
+    if (isseq) {ierr = VecView_Seq(v, viewer);CHKERRQ(ierr);}
+    else       {ierr = VecView_MPI(v, viewer);CHKERRQ(ierr);}
   }
   PetscFunctionReturn(0);
 }
 
+#if defined(PETSC_HAVE_HDF5)
+#undef __FUNCT__
+#define __FUNCT__ "VecView_Plex_HDF5"
+static PetscErrorCode VecView_Plex_HDF5(Vec v, PetscViewer viewer)
+{
+  DM             dm;
+  Vec            locv;
+  const char    *name;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  ierr = VecGetDM(v, &dm);CHKERRQ(ierr);
+  ierr = DMGetLocalVector(dm, &locv);CHKERRQ(ierr);
+  ierr = PetscObjectGetName((PetscObject) v, &name);CHKERRQ(ierr);
+  ierr = PetscObjectSetName((PetscObject) locv, name);CHKERRQ(ierr);
+  ierr = DMGlobalToLocalBegin(dm, v, INSERT_VALUES, locv);CHKERRQ(ierr);
+  ierr = DMGlobalToLocalEnd(dm, v, INSERT_VALUES, locv);CHKERRQ(ierr);
+  ierr = PetscViewerHDF5PushGroup(viewer, "/fields");CHKERRQ(ierr);
+  ierr = PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_VIZ);CHKERRQ(ierr);
+  ierr = VecView_Plex_Local(locv, viewer);CHKERRQ(ierr);
+  ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);
+  ierr = PetscViewerHDF5PopGroup(viewer);CHKERRQ(ierr);
+  ierr = DMRestoreLocalVector(dm, &locv);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+#endif
+
 #undef __FUNCT__
 #define __FUNCT__ "VecView_Plex"
 PetscErrorCode VecView_Plex(Vec v, PetscViewer viewer)
 {
   DM             dm;
-  PetscBool      isvtk;
+  PetscBool      isvtk, ishdf5, isseq;
   PetscErrorCode ierr;
 
   PetscFunctionBegin;
   ierr = VecGetDM(v, &dm);CHKERRQ(ierr);
   if (!dm) SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_WRONG, "Vector not generated from a DM");
-  ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERVTK, &isvtk);CHKERRQ(ierr);
+  ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERVTK,  &isvtk);CHKERRQ(ierr);
+  ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);CHKERRQ(ierr);
+  ierr = PetscObjectTypeCompare((PetscObject) v, VECSEQ, &isseq);CHKERRQ(ierr);
   if (isvtk) {
     Vec         locv;
     const char *name;
     ierr = DMGlobalToLocalEnd(dm, v, INSERT_VALUES, locv);CHKERRQ(ierr);
     ierr = VecView_Plex_Local(locv, viewer);CHKERRQ(ierr);
     ierr = DMRestoreLocalVector(dm, &locv);CHKERRQ(ierr);
+  } else if (ishdf5) {
+#if defined(PETSC_HAVE_HDF5)
+    ierr = VecView_Plex_HDF5(v, viewer);CHKERRQ(ierr);
+#else
+    SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
+#endif
   } else {
-    PetscBool isseq;
+    if (isseq) {ierr = VecView_Seq(v, viewer);CHKERRQ(ierr);}
+    else       {ierr = VecView_MPI(v, viewer);CHKERRQ(ierr);}
+  }
+  PetscFunctionReturn(0);
+}
 
-    ierr = PetscObjectTypeCompare((PetscObject) v, VECSEQ, &isseq);CHKERRQ(ierr);
-    if (isseq) {
-      ierr = VecView_Seq(v, viewer);CHKERRQ(ierr);
-    } else {
-      ierr = VecView_MPI(v, viewer);CHKERRQ(ierr);
-    }
+#undef __FUNCT__
+#define __FUNCT__ "VecLoad_Plex_Local"
+PetscErrorCode VecLoad_Plex_Local(Vec v, PetscViewer viewer)
+{
+  DM             dm;
+  PetscBool      ishdf5;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  ierr = VecGetDM(v, &dm);CHKERRQ(ierr);
+  if (!dm) SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_WRONG, "Vector not generated from a DM");
+  ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);CHKERRQ(ierr);
+  if (ishdf5) {
+    DM          dmBC;
+    Vec         gv;
+    const char *name;
+
+    ierr = DMGetOutputDM(dm, &dmBC);CHKERRQ(ierr);
+    ierr = DMGetGlobalVector(dmBC, &gv);CHKERRQ(ierr);
+    ierr = PetscObjectGetName((PetscObject) v, &name);CHKERRQ(ierr);
+    ierr = PetscObjectSetName((PetscObject) gv, name);CHKERRQ(ierr);
+    ierr = VecLoad_Default(gv, viewer);CHKERRQ(ierr);
+    ierr = DMGlobalToLocalBegin(dmBC, gv, INSERT_VALUES, v);CHKERRQ(ierr);
+    ierr = DMGlobalToLocalEnd(dmBC, gv, INSERT_VALUES, v);CHKERRQ(ierr);
+    ierr = DMRestoreGlobalVector(dmBC, &gv);CHKERRQ(ierr);
+  } else {
+    ierr = VecLoad_Default(v, viewer);CHKERRQ(ierr);
+  }
+  PetscFunctionReturn(0);
+}
+
+#if defined(PETSC_HAVE_HDF5)
+#undef __FUNCT__
+#define __FUNCT__ "VecLoad_Plex_HDF5"
+static PetscErrorCode VecLoad_Plex_HDF5(Vec v, PetscViewer viewer)
+{
+  DM             dm;
+  Vec            locv;
+  const char    *name;
+  PetscInt       seqnum;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  ierr = VecGetDM(v, &dm);CHKERRQ(ierr);
+  ierr = DMGetLocalVector(dm, &locv);CHKERRQ(ierr);
+  ierr = PetscObjectGetName((PetscObject) v, &name);CHKERRQ(ierr);
+  ierr = PetscObjectSetName((PetscObject) locv, name);CHKERRQ(ierr);
+  ierr = DMGetOutputSequenceNumber(dm, &seqnum);CHKERRQ(ierr);
+  ierr = PetscViewerHDF5SetTimestep(viewer, seqnum);CHKERRQ(ierr);
+  ierr = PetscViewerHDF5PushGroup(viewer, "/fields");CHKERRQ(ierr);
+  ierr = VecLoad_Plex_Local(locv, viewer);CHKERRQ(ierr);
+  ierr = PetscViewerHDF5PopGroup(viewer);CHKERRQ(ierr);
+  ierr = DMLocalToGlobalBegin(dm, locv, INSERT_VALUES, v);CHKERRQ(ierr);
+  ierr = DMLocalToGlobalEnd(dm, locv, INSERT_VALUES, v);CHKERRQ(ierr);
+  ierr = DMRestoreLocalVector(dm, &locv);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+#endif
+
+#undef __FUNCT__
+#define __FUNCT__ "VecLoad_Plex"
+PetscErrorCode VecLoad_Plex(Vec v, PetscViewer viewer)
+{
+  DM             dm;
+  PetscBool      ishdf5;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  ierr = VecGetDM(v, &dm);CHKERRQ(ierr);
+  if (!dm) SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_WRONG, "Vector not generated from a DM");
+  ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);CHKERRQ(ierr);
+  if (ishdf5) {
+#if defined(PETSC_HAVE_HDF5)
+    ierr = VecLoad_Plex_HDF5(v, viewer);CHKERRQ(ierr);
+#else
+    SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
+#endif
+  } else {
+    ierr = VecLoad_Default(v, viewer);CHKERRQ(ierr);
   }
   PetscFunctionReturn(0);
 }
   PetscFunctionReturn(0);
 }
 
+#if defined(PETSC_HAVE_HDF5)
+#undef __FUNCT__
+#define __FUNCT__ "DMPlexView_HDF5"
+/* We only write cells and vertices. Does this screw up parallel reading? */
+static PetscErrorCode DMPlexView_HDF5(DM dm, PetscViewer viewer)
+{
+  DM              cdm;
+  Vec             coordinates, newcoords;
+  Vec             coneVec, cellVec;
+  IS              globalVertexNumbers;
+  const PetscInt *gvertex;
+  PetscScalar    *sizes, *vertices;
+  PetscReal       lengthScale;
+  const char     *label   = NULL;
+  PetscInt        labelId = 0, dim;
+  char            group[PETSC_MAX_PATH_LEN];
+  PetscInt        vStart, vEnd, v, cellHeight, cStart, cEnd, cMax, cell, conesSize = 0, numCornersLocal = 0, numCorners, numLabels, l;
+  hid_t           fileId, groupId;
+  herr_t          status;
+  PetscErrorCode  ierr;
+
+  PetscFunctionBegin;
+  /* Write coordinates */
+  ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
+  ierr = DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);CHKERRQ(ierr);
+  ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr);
+  ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
+  ierr = VecDuplicate(coordinates, &newcoords);CHKERRQ(ierr);
+  ierr = PetscObjectSetName((PetscObject) newcoords, "vertices");CHKERRQ(ierr);
+  ierr = VecCopy(coordinates, newcoords);CHKERRQ(ierr);
+  ierr = VecScale(newcoords, lengthScale);CHKERRQ(ierr);
+  /* Use the local version to bypass the default group setting */
+  ierr = PetscViewerHDF5PushGroup(viewer, "/geometry");CHKERRQ(ierr);
+  ierr = VecView(newcoords, viewer);CHKERRQ(ierr);
+  ierr = PetscViewerHDF5PopGroup(viewer);CHKERRQ(ierr);
+  ierr = VecDestroy(&newcoords);CHKERRQ(ierr);
+  /* Write toplogy */
+  ierr = DMPlexGetVTKCellHeight(dm, &cellHeight);CHKERRQ(ierr);
+  ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
+  ierr = DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);CHKERRQ(ierr);
+  ierr = DMPlexGetHybridBounds(dm, &cMax, PETSC_NULL, PETSC_NULL, PETSC_NULL);CHKERRQ(ierr);
+  if (cMax >= 0) cEnd = PetscMin(cEnd, cMax);
+
+  ierr = VecCreate(PetscObjectComm((PetscObject) dm), &coneVec);CHKERRQ(ierr);
+  ierr = VecSetSizes(coneVec, cEnd-cStart, PETSC_DETERMINE);CHKERRQ(ierr);
+  ierr = VecSetBlockSize(coneVec, 1);CHKERRQ(ierr);
+  ierr = VecSetFromOptions(coneVec);CHKERRQ(ierr);
+  ierr = PetscObjectSetName((PetscObject) coneVec, "coneSize");CHKERRQ(ierr);
+  ierr = VecGetArray(coneVec, &sizes);CHKERRQ(ierr);
+  for (cell = cStart; cell < cEnd; ++cell) {
+    PetscInt *closure = NULL;
+    PetscInt  closureSize, v, Nc = 0;
+
+    if (label) {
+      PetscInt value;
+      ierr = DMPlexGetLabelValue(dm, label, cell, &value);CHKERRQ(ierr);
+      if (value == labelId) continue;
+    }
+    ierr = DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
+    for (v = 0; v < closureSize*2; v += 2) {
+      if ((closure[v] >= vStart) && (closure[v] < vEnd)) ++Nc;
+    }
+    ierr = DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
+    conesSize += Nc;
+    if (!numCornersLocal)           numCornersLocal = Nc;
+    else if (numCornersLocal != Nc) numCornersLocal = 1;
+  }
+  ierr = VecRestoreArray(coneVec, &sizes);CHKERRQ(ierr);
+  ierr = MPI_Allreduce(&numCornersLocal, &numCorners, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject) dm));CHKERRQ(ierr);
+  if (numCornersLocal && numCornersLocal != numCorners) numCorners = 1;
+
+  ierr = DMPlexGetVertexNumbering(dm, &globalVertexNumbers);CHKERRQ(ierr);
+  ierr = ISGetIndices(globalVertexNumbers, &gvertex);CHKERRQ(ierr);
+  ierr = VecCreate(PetscObjectComm((PetscObject) dm), &cellVec);CHKERRQ(ierr);
+  ierr = VecSetSizes(cellVec, conesSize, PETSC_DETERMINE);CHKERRQ(ierr);
+  ierr = VecSetBlockSize(cellVec, numCorners);CHKERRQ(ierr);
+  ierr = VecSetFromOptions(cellVec);CHKERRQ(ierr);
+  ierr = PetscObjectSetName((PetscObject) cellVec, "cells");CHKERRQ(ierr);
+  ierr = VecGetArray(cellVec, &vertices);CHKERRQ(ierr);
+  for (cell = cStart, v = 0; cell < cEnd; ++cell) {
+    PetscInt *closure = NULL;
+    PetscInt  closureSize, Nc = 0, p;
+
+    if (label) {
+      PetscInt value;
+      ierr = DMPlexGetLabelValue(dm, label, cell, &value);CHKERRQ(ierr);
+      if (value == labelId) continue;
+    }
+    ierr = DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
+    for (p = 0; p < closureSize*2; p += 2) {
+      if ((closure[p] >= vStart) && (closure[p] < vEnd)) {
+        closure[Nc++] = closure[p];
+        }
+    }
+    ierr = DMPlexInvertCell(dim, Nc, closure);CHKERRQ(ierr);
+    for (p = 0; p < Nc; ++p) {
+      const PetscInt gv = gvertex[closure[p] - vStart];
+      vertices[v++] = gv < 0 ? -(gv+1) : gv;
+    }
+    ierr = DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
+  }
+  if (v != conesSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of cell vertices %d != %d", v, conesSize);
+  ierr = VecRestoreArray(cellVec, &vertices);CHKERRQ(ierr);
+  ierr = PetscViewerHDF5PushGroup(viewer, "/topology");CHKERRQ(ierr);
+  ierr = VecView(cellVec, viewer);CHKERRQ(ierr);
+  if (numCorners == 1) {
+    ierr = VecView(coneVec, viewer);CHKERRQ(ierr);
+  } else {
+    ierr = PetscViewerHDF5WriteAttribute(viewer, "/topology/cells", "cell_corners", PETSC_INT, (void *) &numCorners);CHKERRQ(ierr);
+  }
+  ierr = PetscViewerHDF5PopGroup(viewer);CHKERRQ(ierr);
+  ierr = VecDestroy(&cellVec);CHKERRQ(ierr);
+  ierr = VecDestroy(&coneVec);CHKERRQ(ierr);
+  ierr = ISRestoreIndices(globalVertexNumbers, &gvertex);CHKERRQ(ierr);
+
+  ierr = PetscViewerHDF5WriteAttribute(viewer, "/topology/cells", "cell_dim", PETSC_INT, (void *) &dim);CHKERRQ(ierr);
+  /* Write Labels*/
+  ierr = PetscViewerHDF5PushGroup(viewer, "/labels");CHKERRQ(ierr);
+  ierr = PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);CHKERRQ(ierr);
+  if (groupId != fileId) {status = H5Gclose(groupId);CHKERRQ(status);}
+  ierr = PetscViewerHDF5PopGroup(viewer);CHKERRQ(ierr);
+  ierr = DMPlexGetNumLabels(dm, &numLabels);CHKERRQ(ierr);
+  for (l = 0; l < numLabels; ++l) {
+    DMLabel         label;
+    const char     *name;
+    IS              valueIS;
+    const PetscInt *values;
+    PetscInt        numValues, v;
+    PetscBool       isDepth;
+
+    ierr = DMPlexGetLabelName(dm, l, &name);CHKERRQ(ierr);
+    ierr = PetscStrncmp(name, "depth", 10, &isDepth);CHKERRQ(ierr);
+    if (isDepth) continue;
+    ierr = DMPlexGetLabel(dm, name, &label);CHKERRQ(ierr);
+    ierr = DMLabelGetNumValues(label, &numValues);CHKERRQ(ierr);
+    ierr = DMLabelGetValueIS(label, &valueIS);CHKERRQ(ierr);
+    ierr = ISGetIndices(valueIS, &values);CHKERRQ(ierr);
+    ierr = PetscSNPrintf(group, PETSC_MAX_PATH_LEN, "/labels/%s", name);CHKERRQ(ierr);
+    ierr = PetscViewerHDF5PushGroup(viewer, group);CHKERRQ(ierr);
+    ierr = PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);CHKERRQ(ierr);
+    if (groupId != fileId) {status = H5Gclose(groupId);CHKERRQ(status);}
+    ierr = PetscViewerHDF5PopGroup(viewer);CHKERRQ(ierr);
+    /* TODO: Need to actually loop over the union of label values, ISAllGather() */
+    for (v = 0; v < numValues; ++v) {
+      IS   stratumIS;
+
+      ierr = PetscSNPrintf(group, PETSC_MAX_PATH_LEN, "/labels/%s/%d", name, values[v]);CHKERRQ(ierr);
+      ierr = DMLabelGetStratumIS(label, values[v], &stratumIS);CHKERRQ(ierr);
+      /* TODO: Need to globalize point names and remove unowned points */
+      ierr = PetscViewerHDF5PushGroup(viewer, group);CHKERRQ(ierr);
+      ierr = ISView(stratumIS, viewer);CHKERRQ(ierr);
+      ierr = PetscViewerHDF5PopGroup(viewer);CHKERRQ(ierr);
+      ierr = ISDestroy(&stratumIS);CHKERRQ(ierr);
+    }
+    ierr = ISRestoreIndices(valueIS, &values);CHKERRQ(ierr);
+    ierr = ISDestroy(&valueIS);CHKERRQ(ierr);
+  }
+  PetscFunctionReturn(0);
+}
+#endif
+
 #undef __FUNCT__
 #define __FUNCT__ "DMView_Plex"
 PetscErrorCode DMView_Plex(DM dm, PetscViewer viewer)
 {
-  PetscBool      iascii, isbinary;
+  PetscBool      iascii, ishdf5;
   PetscErrorCode ierr;
 
   PetscFunctionBegin;
   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
   PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2);
   ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr);
-  ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERBINARY, &isbinary);CHKERRQ(ierr);
+  ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5,  &ishdf5);CHKERRQ(ierr);
   if (iascii) {
     ierr = DMPlexView_Ascii(dm, viewer);CHKERRQ(ierr);
-#if 0
-  } else if (isbinary) {
-    ierr = DMPlexView_Binary(dm, viewer);CHKERRQ(ierr);
+  } else if (ishdf5) {
+#if defined(PETSC_HAVE_HDF5)
+    ierr = DMPlexView_HDF5(dm, viewer);CHKERRQ(ierr);
+#else
+    SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
+#endif
+  }
+  PetscFunctionReturn(0);
+}
+
+#if defined(PETSC_HAVE_HDF5)
+typedef struct {
+  DM          dm;
+  PetscViewer viewer;
+  DMLabel     label;
+} LabelCtx;
+
+static herr_t ReadLabelStratumHDF5_Static(hid_t g_id, const char *name, const H5L_info_t *info, void *op_data)
+{
+  PetscViewer     viewer = ((LabelCtx *) op_data)->viewer;
+  DMLabel         label  = ((LabelCtx *) op_data)->label;
+  IS              stratumIS;
+  const PetscInt *ind;
+  PetscInt        value, N, i;
+  const char     *lname;
+  char            group[PETSC_MAX_PATH_LEN];
+  PetscErrorCode  ierr;
+
+  ierr = PetscOptionsStringToInt(name, &value);
+  ierr = ISCreate(PetscObjectComm((PetscObject) viewer), &stratumIS);
+  ierr = PetscObjectSetName((PetscObject) stratumIS, "indices");
+  ierr = DMLabelGetName(label, &lname);
+  ierr = PetscSNPrintf(group, PETSC_MAX_PATH_LEN, "/labels/%s/%s", lname, name);CHKERRQ(ierr);
+  ierr = PetscViewerHDF5PushGroup(viewer, group);CHKERRQ(ierr);
+  ierr = ISLoad(stratumIS, viewer);
+  ierr = PetscViewerHDF5PopGroup(viewer);CHKERRQ(ierr);
+  ierr = ISGetSize(stratumIS, &N);
+  ierr = ISGetIndices(stratumIS, &ind);
+  for (i = 0; i < N; ++i) {ierr = DMLabelSetValue(label, ind[i], value);}
+  ierr = ISRestoreIndices(stratumIS, &ind);
+  ierr = ISDestroy(&stratumIS);
+  return 0;
+}
+
+static herr_t ReadLabelHDF5_Static(hid_t g_id, const char *name, const H5L_info_t *info, void *op_data)
+{
+  DM             dm = ((LabelCtx *) op_data)->dm;
+  hsize_t        idx;
+  herr_t         status;
+  PetscErrorCode ierr;
+
+  ierr = DMPlexCreateLabel(dm, name); if (ierr) return (herr_t) ierr;
+  ierr = DMPlexGetLabel(dm, name, &((LabelCtx *) op_data)->label); if (ierr) return (herr_t) ierr;
+  status = H5Literate_by_name(g_id, name, H5_INDEX_NAME, H5_ITER_NATIVE, &idx, ReadLabelStratumHDF5_Static, op_data, 0);
+  return status;
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "DMPlexLoad_HDF5"
+/* The first version will read everything onto proc 0, letting the user distribute
+   The next will create a naive partition, and then rebalance after reading
+*/
+static PetscErrorCode DMPlexLoad_HDF5(DM dm, PetscViewer viewer)
+{
+  LabelCtx        ctx;
+  PetscSection    coordSection;
+  Vec             coordinates;
+  Vec             cellVec;
+  PetscScalar    *cells;
+  PetscReal       lengthScale;
+  PetscInt       *cone;
+  PetscInt        dim, spatialDim, numVertices, v, numCorners, numCells, cell, c;
+  hid_t           fileId, groupId;
+  hsize_t         idx;
+  herr_t          status;
+  PetscErrorCode  ierr;
+
+  PetscFunctionBegin;
+  /* Read geometry */
+  ierr = PetscViewerHDF5PushGroup(viewer, "/geometry");CHKERRQ(ierr);
+  ierr = VecCreate(PetscObjectComm((PetscObject) dm), &coordinates);CHKERRQ(ierr);
+  ierr = PetscObjectSetName((PetscObject) coordinates, "vertices");CHKERRQ(ierr);
+  ierr = VecLoad(coordinates, viewer);CHKERRQ(ierr);
+  ierr = PetscViewerHDF5PopGroup(viewer);CHKERRQ(ierr);
+  ierr = DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);CHKERRQ(ierr);
+  ierr = VecScale(coordinates, 1.0/lengthScale);CHKERRQ(ierr);
+  ierr = VecGetSize(coordinates, &numVertices);CHKERRQ(ierr);
+  ierr = VecGetBlockSize(coordinates, &spatialDim);CHKERRQ(ierr);
+  numVertices /= spatialDim;
+  /* Read toplogy */
+  ierr = PetscViewerHDF5ReadAttribute(viewer, "/topology/cells", "cell_dim", PETSC_INT, (void *) &dim);CHKERRQ(ierr);
+  ierr = DMPlexSetDimension(dm, dim);CHKERRQ(ierr);
+  /*   TODO Check for coneSize vector rather than this attribute */
+  ierr = PetscViewerHDF5ReadAttribute(viewer, "/topology/cells", "cell_corners", PETSC_INT, (void *) &numCorners);CHKERRQ(ierr);
+  ierr = PetscViewerHDF5PushGroup(viewer, "/topology");CHKERRQ(ierr);
+  ierr = VecCreate(PetscObjectComm((PetscObject) dm), &cellVec);CHKERRQ(ierr);
+  ierr = VecSetBlockSize(cellVec, numCorners);CHKERRQ(ierr);
+  ierr = PetscObjectSetName((PetscObject) cellVec, "cells");CHKERRQ(ierr);
+  ierr = VecLoad(cellVec, viewer);CHKERRQ(ierr);
+  ierr = PetscViewerHDF5PopGroup(viewer);CHKERRQ(ierr);
+  ierr = VecGetSize(cellVec, &numCells);CHKERRQ(ierr);
+  numCells /= numCorners;
+  /* Create Plex */
+  ierr = DMPlexSetChart(dm, 0, numCells+numVertices);CHKERRQ(ierr);
+  for (cell = 0; cell < numCells; ++cell) {ierr = DMPlexSetConeSize(dm, cell, numCorners);CHKERRQ(ierr);}
+  ierr = DMSetUp(dm);CHKERRQ(ierr);
+  ierr = PetscMalloc1(numCorners,&cone);CHKERRQ(ierr);
+  ierr = VecGetArray(cellVec, &cells);CHKERRQ(ierr);
+  for (cell = 0; cell < numCells; ++cell) {
+    for (c = 0; c < numCorners; ++c) {cone[c] = numCells + cells[cell*numCorners+c];}
+    ierr = DMPlexSetCone(dm, cell, cone);CHKERRQ(ierr);
+  }
+  ierr = VecRestoreArray(cellVec, &cells);CHKERRQ(ierr);
+  ierr = PetscFree(cone);CHKERRQ(ierr);
+  ierr = VecDestroy(&cellVec);CHKERRQ(ierr);
+  ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr);
+  ierr = DMPlexStratify(dm);CHKERRQ(ierr);
+  /* Create coordinates */
+  ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
+  ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr);
+  ierr = PetscSectionSetFieldComponents(coordSection, 0, spatialDim);CHKERRQ(ierr);
+  ierr = PetscSectionSetChart(coordSection, numCells, numCells+numVertices);CHKERRQ(ierr);
+  for (v = numCells; v < numCells+numVertices; ++v) {
+    ierr = PetscSectionSetDof(coordSection, v, spatialDim);CHKERRQ(ierr);
+    ierr = PetscSectionSetFieldDof(coordSection, v, 0, spatialDim);CHKERRQ(ierr);
+  }
+  ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr);
+  ierr = DMSetCoordinates(dm, coordinates);CHKERRQ(ierr);
+  ierr = VecDestroy(&coordinates);CHKERRQ(ierr);
+  /* Read Labels*/
+  ctx.dm     = dm;
+  ctx.viewer = viewer;
+  ierr = PetscViewerHDF5PushGroup(viewer, "/labels");CHKERRQ(ierr);
+  ierr = PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);CHKERRQ(ierr);
+  status = H5Literate(groupId, H5_INDEX_NAME, H5_ITER_NATIVE, &idx, ReadLabelHDF5_Static, &ctx);CHKERRQ(status);
+  status = H5Gclose(groupId);CHKERRQ(status);
+  ierr = PetscViewerHDF5PopGroup(viewer);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+#endif
+
+#undef __FUNCT__
+#define __FUNCT__ "DMLoad_Plex"
+PetscErrorCode DMLoad_Plex(DM dm, PetscViewer viewer)
+{
+  PetscBool      isbinary, ishdf5;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
+  PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2);
+  ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERBINARY, &isbinary);CHKERRQ(ierr);
+  ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5,   &ishdf5);CHKERRQ(ierr);
+  if (isbinary) {SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Do not yet support binary viewers");}
+  else if (ishdf5) {
+#if defined(PETSC_HAVE_HDF5)
+    DM odm;
+    ierr = DMCreate(PetscObjectComm((PetscObject) dm), &odm);CHKERRQ(ierr);
+    ierr = DMSetType(odm, DMPLEX);CHKERRQ(ierr);
+    ierr = DMPlexLoad_HDF5(odm, viewer);CHKERRQ(ierr);
+    ierr = DMPlexInterpolate(odm, &dm);CHKERRQ(ierr);
+    ierr = DMPlexCopyCoordinates(odm, dm);CHKERRQ(ierr);
+    ierr = DMPlexCopyLabels(odm, dm);CHKERRQ(ierr);
+    ierr = DMDestroy(&odm);CHKERRQ(ierr);
+#else
+    SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
 #endif
   }
   PetscFunctionReturn(0);
 . numDof    - An array of size numFields*(dim+1) which holds the number of dof for each field on a mesh piece of dimension d
 . numBC     - The number of boundary conditions
 . bcField   - An array of size numBC giving the field number for each boundry condition
-- bcPoints  - An array of size numBC giving an IS holding the sieve points to which each boundary condition applies
+. bcPoints  - An array of size numBC giving an IS holding the sieve points to which each boundary condition applies
+- perm      - Optional permutation of the chart, or NULL
 
   Output Parameter:
 . section - The PetscSection object
 
   Notes: numDof[f*(dim+1)+d] gives the number of dof for field f on sieve points of dimension d. For instance, numDof[1] is the
-  nubmer of dof for field 0 on each edge.
+  number of dof for field 0 on each edge.
+
+  The chart permutation is the same one set using PetscSectionSetPermutation()
 
   Level: developer
 
   A Fortran 90 version is available as DMPlexCreateSectionF90()
 
 .keywords: mesh, elements
-.seealso: DMPlexCreate(), PetscSectionCreate()
+.seealso: DMPlexCreate(), PetscSectionCreate(), PetscSectionSetPermutation()
 @*/
-PetscErrorCode DMPlexCreateSection(DM dm, PetscInt dim, PetscInt numFields,const PetscInt numComp[],const PetscInt numDof[], PetscInt numBC,const PetscInt bcField[],const IS bcPoints[], PetscSection *section)
+PetscErrorCode DMPlexCreateSection(DM dm, PetscInt dim, PetscInt numFields,const PetscInt numComp[],const PetscInt numDof[], PetscInt numBC,const PetscInt bcField[],const IS bcPoints[], IS perm, PetscSection *section)
 {
   PetscErrorCode ierr;
 
   PetscFunctionBegin;
   ierr = DMPlexCreateSectionInitial(dm, dim, numFields, numComp, numDof, section);CHKERRQ(ierr);
   ierr = DMPlexCreateSectionBCDof(dm, numBC, bcField, bcPoints, PETSC_DETERMINE, *section);CHKERRQ(ierr);
+  if (perm) {ierr = PetscSectionSetPermutation(*section, perm);CHKERRQ(ierr);}
   ierr = PetscSectionSetUp(*section);CHKERRQ(ierr);
   if (numBC) {ierr = DMPlexCreateSectionBCIndicesAll(dm, *section);CHKERRQ(ierr);}
   ierr = PetscSectionViewFromOptions(*section,NULL,"-section_view");CHKERRQ(ierr);
   PetscErrorCode ierr;
 
   PetscFunctionBegin;
-  ierr = PetscSectionCreate(s->atlasLayout.comm, gsection);CHKERRQ(ierr);
+  ierr = PetscSectionCreate(PetscObjectComm((PetscObject) s), gsection);CHKERRQ(ierr);
   ierr = PetscSectionGetChart(s, &pStart, &pEnd);CHKERRQ(ierr);
   ierr = PetscSectionSetChart(*gsection, pStart, pEnd);CHKERRQ(ierr);
   ierr = PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);CHKERRQ(ierr);
     (*gsection)->atlasOff[p] = off;
     off += (*gsection)->atlasDof[p] > 0 ? (*gsection)->atlasDof[p]-cdof : 0;
   }
-  ierr       = MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, s->atlasLayout.comm);CHKERRQ(ierr);
+  ierr       = MPI_Scan(&off, &globalOff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) s));CHKERRQ(ierr);
   globalOff -= off;
   for (p = 0, off = 0; p < pEnd-pStart; ++p) {
     (*gsection)->atlasOff[p] += globalOff;
       if ((numDof[f*(dim+1)+d] > 0) && (depth < dim)) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mesh must be interpolated when unknowns are specified on edges or faces.");
     }
   }
-  ierr = DMPlexCreateSection(dm, dim, numFields, numComp, numDof, numBC, bcFields, bcPoints, &section);CHKERRQ(ierr);
+  ierr = DMPlexCreateSection(dm, dim, numFields, numComp, numDof, numBC, bcFields, bcPoints, NULL, &section);CHKERRQ(ierr);
   for (f = 0; f < numFields; ++f) {
     PetscFE     fe;
     const char *name;

src/dm/impls/plex/plexcgns.c

   ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr);
   ierr = DMPlexStratify(*dm);CHKERRQ(ierr);
   if (interpolate) {
-    DM idm;
+    DM idm = NULL;
 
     ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr);
     /* Maintain zone label */

src/dm/impls/plex/plexcreate.c

     *newdm = rdm;
   }
   if (interpolate) {
-    DM idm;
+    DM idm = NULL;
     const char *name;
 
     ierr = DMPlexInterpolate(*newdm, &idm);CHKERRQ(ierr);
   ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr);
   ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr);
   ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr);
+  ierr = VecSetBlockSize(coordinates, 2);CHKERRQ(ierr);
   ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr);
   ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr);
   ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr);
   ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr);
   ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr);
   ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr);
+  ierr = VecSetBlockSize(coordinates, 3);CHKERRQ(ierr);
   ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr);
   ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr);
   ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr);
     ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr);
     ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr);
     ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr);
+    ierr = VecSetBlockSize(coordinates, 2);CHKERRQ(ierr);
     ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr);
     ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr);
     ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr);
 extern PetscErrorCode DMSetUp_Plex(DM dm);
 extern PetscErrorCode DMDestroy_Plex(DM dm);
 extern PetscErrorCode DMView_Plex(DM dm, PetscViewer viewer);
+extern PetscErrorCode DMLoad_Plex(DM dm, PetscViewer viewer);
 extern PetscErrorCode DMCreateSubDM_Plex(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm);
 extern PetscErrorCode DMLocatePoints_Plex(DM dm, Vec v, IS *cellIS);
 
   PetscFunctionBegin;
   ierr = DMCreateGlobalVector_Section_Private(dm,vec);CHKERRQ(ierr);
   /* ierr = VecSetOperation(*vec, VECOP_DUPLICATE, (void(*)(void)) VecDuplicate_MPI_DM);CHKERRQ(ierr); */
-  ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void))VecView_Plex);CHKERRQ(ierr);
+  ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex);CHKERRQ(ierr);
+  ierr = VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
 
 
   PetscFunctionBegin;
   ierr = DMCreateLocalVector_Section_Private(dm,vec);CHKERRQ(ierr);
-  ierr = VecSetOperation(*vec, VECOP_VIEW, (void(*)(void)) VecView_Plex_Local);CHKERRQ(ierr);
+  ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex_Local);CHKERRQ(ierr);
+  ierr = VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex_Local);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
 
 {
   PetscFunctionBegin;
   dm->ops->view                            = DMView_Plex;
+  dm->ops->load                            = DMLoad_Plex;
   dm->ops->setfromoptions                  = DMSetFromOptions_Plex;
   dm->ops->clone                           = DMClone_Plex;
   dm->ops->setup                           = DMSetUp_Plex;
   ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr);
   ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr);
   ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr);
+  ierr = VecSetBlockSize(coordinates, spaceDim);CHKERRQ(ierr);
   ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr);
   ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr);
   ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr);
   ierr = DMPlexSetDimension(*dm, dim);CHKERRQ(ierr);
   ierr = DMPlexBuildFromCellList_Private(*dm, numCells, numVertices, numCorners, cells);CHKERRQ(ierr);
   if (interpolate) {
-    DM idm;
+    DM idm = NULL;
 
     ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr);
     ierr = DMDestroy(dm);CHKERRQ(ierr);
   ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr);
   ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr);
   ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr);
+  ierr = VecSetBlockSize(coordinates, dim);CHKERRQ(ierr);
   ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr);
   ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr);
   ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr);

src/dm/impls/plex/plexdistribute.c

 
 #undef __FUNCT__
 #define __FUNCT__ "DMPlexGetAdjacency_Internal"
-PetscErrorCode DMPlexGetAdjacency_Internal(DM dm, PetscInt p, PetscBool useCone, PetscBool useTransitiveClosure, PetscInt *adjSize, PetscInt adj[])
+PetscErrorCode DMPlexGetAdjacency_Internal(DM dm, PetscInt p, PetscBool useCone, PetscBool useTransitiveClosure, PetscInt *adjSize, PetscInt *adj[])
 {
+  static PetscInt asiz = 0;
   PetscErrorCode  ierr;
 
   PetscFunctionBeginHot;
+  if (!*adj) {
+    PetscInt depth, maxConeSize, maxSupportSize;
+
+    ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr);
+    ierr = DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);CHKERRQ(ierr);
+    asiz = PetscPowInt(maxConeSize, depth+1) * PetscPowInt(maxSupportSize, depth+1) + 1;
+    ierr = PetscMalloc1(asiz,adj);CHKERRQ(ierr);
+  }
+  if (*adjSize < 0) *adjSize = asiz;
   if (useTransitiveClosure) {
-    ierr = DMPlexGetAdjacency_Transitive_Internal(dm, p, useCone, adjSize, adj);CHKERRQ(ierr);
+    ierr = DMPlexGetAdjacency_Transitive_Internal(dm, p, useCone, adjSize, *adj);CHKERRQ(ierr);
   } else if (useCone) {
-    ierr = DMPlexGetAdjacency_Cone_Internal(dm, p, adjSize, adj);CHKERRQ(ierr);
+    ierr = DMPlexGetAdjacency_Cone_Internal(dm, p, adjSize, *adj);CHKERRQ(ierr);
   } else {
-    ierr = DMPlexGetAdjacency_Support_Internal(dm, p, adjSize, adj);CHKERRQ(ierr);
+    ierr = DMPlexGetAdjacency_Support_Internal(dm, p, adjSize, *adj);CHKERRQ(ierr);
   }
   PetscFunctionReturn(0);
 }
 @*/
 PetscErrorCode DMPlexGetAdjacency(DM dm, PetscInt p, PetscInt *adjSize, PetscInt *adj[])
 {
-  DM_Plex        *mesh = (DM_Plex *) dm->data;
-  static PetscInt asiz = 0;
-  PetscErrorCode  ierr;
+  DM_Plex       *mesh = (DM_Plex *) dm->data;
+  PetscErrorCode ierr;
 
   PetscFunctionBeginHot;
   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
   PetscValidPointer(adjSize,3);
   PetscValidPointer(adj,4);
-  if (!*adj) {
-    PetscInt depth, maxConeSize, maxSupportSize;
-
-    ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr);
-    ierr = DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);CHKERRQ(ierr);
-    asiz = PetscPowInt(maxConeSize, depth+1) * PetscPowInt(maxSupportSize, depth+1) + 1;
-    ierr = PetscMalloc1(asiz,adj);CHKERRQ(ierr);
-  }
-  if (*adjSize < 0) *adjSize = asiz;
-  ierr = DMPlexGetAdjacency_Internal(dm, p, mesh->useCone, mesh->useClosure, adjSize, *adj);CHKERRQ(ierr);
+  ierr = DMPlexGetAdjacency_Internal(dm, p, mesh->useCone, mesh->useClosure, adjSize, adj);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
 

src/dm/impls/plex/plexexodusii.c

   ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr);
   ierr = DMPlexStratify(*dm);CHKERRQ(ierr);
   if (interpolate) {
-    DM idm;
+    DM idm = NULL;
 
     ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr);
     /* Maintain Cell Sets label */