Commits

Matt Knepley committed fd0828a Merge

Merge branch 'knepley/feature-plex-hdf5-parallel-load'

* knepley/feature-plex-hdf5-parallel-load:
HDF5: Fixed Xdmf script
Plex+HDF5: Do not use vectors for periodic visualization
Plex+HDF5: Ignore negative sequence numbers
DMPlex: Now Plex output is parallel - We now write visualization specific topology in /viz/topology - We output a point reordering, coneSizes, cones, and orientations - Now longer need to interpolate on load
DMPlex: Use the presence of faceGeometry in the DM to signal we are using FVM - This will become a PetscFVM object soon
DMPlex: Preserve the block size of the coordinate vector after distribution
Plex+HDF5: Moved all HDF5 to a separate file, and mapped 2D periodic mesh to the cylinder - Now visualization specific things are in /viz
HDF5: Made petsc_gen_xdmf.py executable
Viewer+Options: Added Fortran interface for PetscObjectViewFromOptions()
IS: Stupid mistake. Damn you compiler
DMPlex: Now parallel HDF5 label output does not fail - However, it is also now clear that we will have to write the full interpolation connectivity in order for these to be meaningful
DMPlex: Added DMPlexCreatePointNumbering() - Made DMPlexCreateNumbering_Private() more flexible
IS: Added ISSortRemoveDups()
DMPlex: Added DMPlexInvertCell_Internal() - Stupid type matching
HDF5: Added petsc_gen_xdmf.py which processes PETSc HDF5 output and produces an Xdmf file
DMPlex: Force a serial load in DMLoad_HDF5(), after which we call DMPlexDistribute() - Eventually we would load into a naive partition
IS+HDF5: Corrected code for new block size interface where it is never negative
Viewer+HDF5: Added PetscViewerHDF5ReadSizes() - This allows me to check the size of a Vec or IS to be loaded

Comments (0)

Files changed (18)

bin/pythonscripts/petsc_gen_xdmf.py

+#!/usr/bin/env python
+import h5py
+import numpy as np
+import os, sys
+
+class Xdmf:
+  def __init__(self, filename):
+    self.filename = filename
+    self.cellMap  = {1 : {1 : 'Polyvertex', 2 : 'Polyline'}, 2 : {3 : 'Triangle', 4 : 'Quadrilateral'}, 3 : {4 : 'Tetrahedron', 8 : 'Hexahedron'}}
+    self.typeMap  = {'scalar' : 'Scalar', 'vector' : 'Vector', 'tensor' : 'Tensor6', 'matrix' : 'Matrix'}
+    self.typeExt  = {2 : {'vector' : ['x', 'y'], 'tensor' : ['xx', 'yy', 'xy']}, 3 : {'vector' : ['x', 'y', 'z'], 'tensor' : ['xx', 'yy', 'zz', 'xy', 'yz', 'xz']}}
+    return
+
+  def writeHeader(self, fp, hdfFilename):
+    fp.write('''\
+<?xml version="1.0" ?>
+<!DOCTYPE Xdmf SYSTEM "Xdmf.dtd" [
+<!ENTITY HeavyData "%s">
+]>
+''' % os.path.basename(hdfFilename))
+    fp.write('\n<Xdmf>\n  <Domain Name="domain">\n')
+    return
+
+  def writeCells(self, fp, topologyPath, numCells, numCorners):
+    fp.write('''\
+    <DataItem Name="cells"
+	      ItemType="Uniform"
+	      Format="HDF"
+	      NumberType="Float" Precision="8"
+	      Dimensions="%d %d">
+      &HeavyData;:/%s/cells
+    </DataItem>
+''' % (numCells, numCorners, topologyPath))
+    return
+
+  def writeVertices(self, fp, geometryPath, numVertices, spaceDim):
+    fp.write('''\
+    <DataItem Name="vertices"
+	      Format="HDF"
+	      Dimensions="%d %d">
+      &HeavyData;:/%s/vertices
+    </DataItem>
+    <!-- ============================================================ -->
+''' % (numVertices, spaceDim, geometryPath))
+    return
+
+  def writeTimeGridHeader(self, fp, time):
+    fp.write('''\
+    <Grid Name="TimeSeries" GridType="Collection" CollectionType="Temporal">
+      <Time TimeType="List">
+        <DataItem Format="XML" NumberType="Float" Dimensions="%d">
+          ''' % (len(time)))
+    fp.write(' '.join(map(str, map(int, time))))
+    fp.write('''
+        </DataItem>
+      </Time>
+''')
+    return
+
+  def writeSpaceGridHeader(self, fp, numCells, numCorners, cellDim, spaceDim):
+    fp.write('''\
+      <Grid Name="domain" GridType="Uniform">
+	<Topology
+	   TopologyType="%s"
+	   NumberOfElements="%d">
+	  <DataItem Reference="XML">
+	    /Xdmf/Domain/DataItem[@Name="cells"]
+	  </DataItem>
+	</Topology>
+	<Geometry GeometryType="%s">
+	  <DataItem Reference="XML">
+	    /Xdmf/Domain/DataItem[@Name="vertices"]
+	  </DataItem>
+	</Geometry>
+''' % (self.cellMap[cellDim][numCorners], numCells, "XYZ" if spaceDim > 2 else "XY"))
+    return
+
+  def writeFieldSingle(self, fp, numSteps, timestep, spaceDim, name, f, domain):
+    if len(f[1].shape) > 2:
+      dof = f[1].shape[1]
+      bs  = f[1].shape[2]
+    else:
+      dof = f[1].shape[0]
+      bs  = f[1].shape[1]
+    fp.write('''\
+	<Attribute
+	   Name="%s"
+	   Type="%s"
+	   Center="%s">
+          <DataItem ItemType="HyperSlab"
+		    Dimensions="1 %d %d"
+		    Type="HyperSlab">
+            <DataItem
+	       Dimensions="3 3"
+	       Format="XML">
+              %d 0 0
+              1 1 1
+              1 %d %d
+	    </DataItem>
+	    <DataItem
+	       DataType="Float" Precision="8"
+	       Dimensions="%d %d %d"
+	       Format="HDF">
+	      &HeavyData;:%s
+	    </DataItem>
+	  </DataItem>
+	</Attribute>
+''' % (f[0], self.typeMap[f[1].attrs['vector_field_type']], domain, dof, bs, timestep, dof, bs, numSteps, dof, bs, name))
+    return
+
+  def writeFieldComponents(self, fp, numSteps, timestep, spaceDim, name, f, domain):
+    vtype = f[1].attrs['vector_field_type']
+    if len(f[1].shape) > 2:
+      dof = f[1].shape[1]
+      bs  = f[1].shape[2]
+    else:
+      dof = f[1].shape[0]
+      bs  = f[1].shape[1]
+    for c in range(bs):
+      ext = self.typeExt[spaceDim][vtype][c]
+      fp.write('''\
+	<Attribute
+	   Name="%s"
+	   Type="Scalar"
+	   Center="%s">
+          <DataItem ItemType="HyperSlab"
+		    Dimensions="1 %d 1"
+		    Type="HyperSlab">
+            <DataItem
+	       Dimensions="3 3"
+	       Format="XML">
+              %d 0 %d
+              1 1 1
+              1 %d 1
+	    </DataItem>
+	    <DataItem
+	       DataType="Float" Precision="8"
+	       Dimensions="%d %d %d"
+	       Format="HDF">
+	      &HeavyData;:%s
+	    </DataItem>
+	  </DataItem>
+	</Attribute>
+''' % (f[0]+'_'+ext, domain, dof, timestep, c, dof, numSteps, dof, bs, name))
+    return
+
+  def writeField(self, fp, numSteps, timestep, cellDim, spaceDim, name, f, domain):
+    ctypes = ['tensor', 'matrix']
+    if spaceDim == 2 or cellDim != spaceDim: ctypes.append('vector')
+    if f[1].attrs['vector_field_type'] in ctypes:
+      self.writeFieldComponents(fp, numSteps, timestep, spaceDim, name, f, domain)
+    else:
+      self.writeFieldSingle(fp, numSteps, timestep, spaceDim, name, f, domain)
+    return
+
+  def writeSpaceGridFooter(self, fp):
+    fp.write('      </Grid>\n')
+    return
+
+  def writeTimeGridFooter(self, fp):
+    fp.write('    </Grid>\n')
+    return
+
+  def writeFooter(self, fp):
+    fp.write('  </Domain>\n</Xdmf>\n')
+    return
+
+  def write(self, hdfFilename, topologyPath, numCells, numCorners, cellDim, geometryPath, numVertices, spaceDim, time, vfields, cfields):
+    useTime = not (len(time) < 2 and time[0] == -1)
+    with file(self.filename, 'w') as fp:
+      self.writeHeader(fp, hdfFilename)
+      self.writeCells(fp, topologyPath, numCells, numCorners)
+      self.writeVertices(fp, geometryPath, numVertices, spaceDim)
+      if useTime: self.writeTimeGridHeader(fp, time)
+      for t in range(len(time)):
+        self.writeSpaceGridHeader(fp, numCells, numCorners, cellDim, spaceDim)
+        for vf in vfields: self.writeField(fp, len(time), t, cellDim, spaceDim, '/vertex_fields/'+vf[0], vf, 'Node')
+        for cf in cfields: self.writeField(fp, len(time), t, cellDim, spaceDim, '/cell_fields/'+cf[0], cf, 'Cell')
+        self.writeSpaceGridFooter(fp)
+      if useTime: self.writeTimeGridFooter(fp)
+      self.writeFooter(fp)
+    return
+
+def generateXdmf(hdfFilename, xdmfFilename = None):
+  if xdmfFilename is None:
+    xdmfFilename = os.path.splitext(hdfFilename)[0] + '.xmf'
+  # Read mesh
+  h5          = h5py.File(hdfFilename, 'r')
+  if 'viz' in h5 and 'geometry' in h5['viz']:
+    geomPath  = 'viz/geometry'
+    geom      = h5['viz']['geometry']
+  else:
+    geomPath  = 'geometry'
+    geom      = h5['geometry']
+  if 'viz' in h5 and 'topology' in h5['viz']:
+    topoPath  = 'viz/topology'
+    topo      = h5['viz']['topology']
+  else:
+    topoPath  = 'topology'
+    topo      = h5['topology']
+  vertices    = geom['vertices']
+  numVertices = vertices.shape[0]
+  spaceDim    = vertices.shape[1]
+  cells       = topo['cells']
+  numCells    = cells.shape[0]
+  numCorners  = cells.shape[1]
+  cellDim     = topo['cells'].attrs['cell_dim']
+  time        = np.array(h5['time']).flatten()
+  vfields     = []
+  cfields     = []
+  if 'vertex_fields' in h5: vfields = h5['vertex_fields'].items()
+  if 'cell_fields' in h5: cfields = h5['cell_fields'].items()
+
+  # Write Xdmf
+  Xdmf(xdmfFilename).write(hdfFilename, topoPath, numCells, numCorners, cellDim, geomPath, numVertices, spaceDim, time, vfields, cfields)
+  h5.close()
+  return
+
+if __name__ == '__main__':
+  generateXdmf(sys.argv[1])

include/petsc-private/dmpleximpl.h

 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 DMPlexGetFieldType_Internal(DM, PetscSection, PetscInt, PetscInt *, PetscInt *, PetscViewerVTKFieldType *);
+#if defined(PETSC_HAVE_HDF5)
+PETSC_EXTERN PetscErrorCode VecView_Plex_Local_HDF5(Vec, PetscViewer);
+PETSC_EXTERN PetscErrorCode VecView_Plex_HDF5(Vec, PetscViewer);
+PETSC_EXTERN PetscErrorCode VecLoad_Plex_HDF5(Vec, PetscViewer);
+PETSC_EXTERN PetscErrorCode DMPlexView_HDF5(DM, PetscViewer);
+PETSC_EXTERN PetscErrorCode DMPlexLoad_HDF5(DM, PetscViewer);
+#endif
 
 PETSC_EXTERN PetscErrorCode DMPlexGetAdjacency_Internal(DM,PetscInt,PetscBool,PetscBool,PetscInt*,PetscInt*[]);
 PETSC_EXTERN PetscErrorCode DMPlexGetFaces_Internal(DM,PetscInt,PetscInt,PetscInt*,PetscInt*,const PetscInt*[]);

include/petsc-private/isimpl.h

   PetscErrorCode (*restoreindices)(IS,const PetscInt*[]);
   PetscErrorCode (*invertpermutation)(IS,PetscInt,IS*);
   PetscErrorCode (*sort)(IS);
+  PetscErrorCode (*sortremovedups)(IS);
   PetscErrorCode (*sorted)(IS,PetscBool*);
   PetscErrorCode (*duplicate)(IS,IS*);
   PetscErrorCode (*destroy)(IS);

include/petscdmplex.h

 PETSC_EXTERN PetscErrorCode DMPlexRemoveLabel(DM, const char [], DMLabel *);
 PETSC_EXTERN PetscErrorCode DMPlexGetCellNumbering(DM, IS *);
 PETSC_EXTERN PetscErrorCode DMPlexGetVertexNumbering(DM, IS *);
+PETSC_EXTERN PetscErrorCode DMPlexCreatePointNumbering(DM, IS *);
 
 PETSC_EXTERN PetscErrorCode DMPlexGetDepth(DM, PetscInt *);
 PETSC_EXTERN PetscErrorCode DMPlexGetDepthLabel(DM, DMLabel *);

include/petscis.h

 PETSC_EXTERN PetscErrorCode ISLoad(IS,PetscViewer);
 PETSC_EXTERN PetscErrorCode ISEqual(IS,IS,PetscBool  *);
 PETSC_EXTERN PetscErrorCode ISSort(IS);
+PETSC_EXTERN PetscErrorCode ISSortRemoveDups(IS);
 PETSC_EXTERN PetscErrorCode ISSorted(IS,PetscBool  *);
 PETSC_EXTERN PetscErrorCode ISDifference(IS,IS,IS*);
 PETSC_EXTERN PetscErrorCode ISSum(IS,IS,IS*);

include/petscviewerhdf5.h

 #include <hdf5.h>
 PETSC_EXTERN PetscErrorCode PetscViewerHDF5GetFileId(PetscViewer,hid_t*);
 PETSC_EXTERN PetscErrorCode PetscViewerHDF5OpenGroup(PetscViewer, hid_t *, hid_t *);
+PETSC_EXTERN PetscErrorCode PetscViewerHDF5ReadSizes(PetscViewer, const char[], PetscInt *, PetscInt *);
 
 /* On 32 bit systems HDF5 is limited by size of integer, because hsize_t is defined as size_t */
 #define PETSC_HDF5_INT_MAX  2147483647

src/dm/impls/plex/makefile

 CPPFLAGS =
 CFLAGS   =
 FFLAGS   =
-SOURCEC  = plexcreate.c plex.c plexpartition.c plexdistribute.c plexrefine.c plexinterpolate.c plexpreallocate.c plexreorder.c plexgeometry.c plexlabel.c plexsubmesh.c plexexodusii.c plexgmsh.c plexcgns.c plexvtk.c plexpoint.c plexvtu.c plexfem.c plexindices.c
+SOURCEC  = plexcreate.c plex.c plexpartition.c plexdistribute.c plexrefine.c plexinterpolate.c plexpreallocate.c plexreorder.c plexgeometry.c plexlabel.c plexsubmesh.c plexhdf5.c plexexodusii.c plexgmsh.c plexcgns.c plexvtk.c plexpoint.c plexvtu.c plexfem.c plexindices.c
 SOURCEF  =
 SOURCEH  =
 DIRS     = examples

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)
+#define __FUNCT__ "DMPlexGetFieldType_Internal"
+PetscErrorCode DMPlexGetFieldType_Internal(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;
   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)
     if (fem) {
       ierr = DMPlexInsertBoundaryValuesFEM(dm, v);CHKERRQ(ierr);
     } else {
+      PetscObject faceGeometry;
+
       /* TODO Fix the time, and add FVM objects */
-      ierr = DMPlexInsertBoundaryValuesFVM(dm, 0.0, v);CHKERRQ(ierr);
+      ierr = PetscObjectQuery((PetscObject) dm, "FaceGeometry", &faceGeometry);CHKERRQ(ierr);
+      if (faceGeometry) {ierr = DMPlexInsertBoundaryValuesFVM(dm, 0.0, v);CHKERRQ(ierr);}
     }
   }
   if (isvtk) {
     PetscInt                pStart, pEnd;
 
     ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
-    ierr = GetFieldType_Static(dm, section, PETSC_DETERMINE, &pStart, &pEnd, &ft);CHKERRQ(ierr);
+    ierr = DMPlexGetFieldType_Internal(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);
   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)
   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)
   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)
     ierr = DMPlexView_Ascii(dm, viewer);CHKERRQ(ierr);
   } else if (ishdf5) {
 #if defined(PETSC_HAVE_HDF5)
+    ierr = PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_VIZ);CHKERRQ(ierr);
     ierr = DMPlexView_HDF5(dm, viewer);CHKERRQ(ierr);
+    ierr = PetscViewerPopFormat(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 = 0;
-  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)
   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);
+    ierr = DMPlexLoad_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
 }
 
 #undef __FUNCT__
+#define __FUNCT__ "DMPlexInvertCell_Internal"
+PetscErrorCode DMPlexInvertCell_Internal(PetscInt dim, PetscInt numCorners, PetscInt cone[])
+{
+  int tmpc;
+
+  PetscFunctionBegin;
+  if (dim != 3) PetscFunctionReturn(0);
+  switch (numCorners) {
+  case 4:
+    tmpc    = cone[0];
+    cone[0] = cone[1];
+    cone[1] = tmpc;
+    break;
+  case 8:
+    tmpc    = cone[1];
+    cone[1] = cone[3];
+    cone[3] = tmpc;
+    break;
+  default: break;
+  }
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
 #define __FUNCT__ "DMPlexInvertCell"
 /*@C
   DMPlexInvertCell - This flips tetrahedron and hexahedron orientation since Plex stores them internally with outward normals. Other cells are left untouched.
 #undef __FUNCT__
 #define __FUNCT__ "DMPlexCreateNumbering_Private"
 /* We can easily have a form that takes an IS instead */
-PetscErrorCode DMPlexCreateNumbering_Private(DM dm, PetscInt pStart, PetscInt pEnd, PetscSF sf, IS *numbering)
+static PetscErrorCode DMPlexCreateNumbering_Private(DM dm, PetscInt pStart, PetscInt pEnd, PetscInt shift, PetscInt *globalSize, PetscSF sf, IS *numbering)
 {
   PetscSection   section, globalSection;
   PetscInt      *numbers, p;
   ierr = PetscMalloc1((pEnd - pStart), &numbers);CHKERRQ(ierr);
   for (p = pStart; p < pEnd; ++p) {
     ierr = PetscSectionGetOffset(globalSection, p, &numbers[p-pStart]);CHKERRQ(ierr);
+    if (numbers[p-pStart] < 0) numbers[p-pStart] -= shift;
+    else                       numbers[p-pStart] += shift;
+  }
+  ierr = ISCreateGeneral(PetscObjectComm((PetscObject) dm), pEnd - pStart, numbers, PETSC_OWN_POINTER, numbering);CHKERRQ(ierr);
+  if (globalSize) {
+    PetscLayout layout;
+    ierr = PetscSectionGetPointLayout(PetscObjectComm((PetscObject) dm), globalSection, &layout);CHKERRQ(ierr);
+    ierr = PetscLayoutGetSize(layout, globalSize);CHKERRQ(ierr);
+    ierr = PetscLayoutDestroy(&layout);CHKERRQ(ierr);
   }
-  ierr = ISCreateGeneral(PetscObjectComm((PetscObject)dm), pEnd - pStart, numbers, PETSC_OWN_POINTER, numbering);CHKERRQ(ierr);
   ierr = PetscSectionDestroy(&section);CHKERRQ(ierr);
   ierr = PetscSectionDestroy(&globalSection);CHKERRQ(ierr);
   PetscFunctionReturn(0);
     ierr = DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);CHKERRQ(ierr);
     ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr);
     if (cMax >= 0) cEnd = PetscMin(cEnd, cMax);
-    ierr = DMPlexCreateNumbering_Private(dm, cStart, cEnd, dm->sf, &mesh->globalCellNumbers);CHKERRQ(ierr);
+    ierr = DMPlexCreateNumbering_Private(dm, cStart, cEnd, 0, NULL, dm->sf, &mesh->globalCellNumbers);CHKERRQ(ierr);
   }
   *globalCellNumbers = mesh->globalCellNumbers;
   PetscFunctionReturn(0);
     ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
     ierr = DMPlexGetHybridBounds(dm, NULL, NULL, NULL, &vMax);CHKERRQ(ierr);
     if (vMax >= 0) vEnd = PetscMin(vEnd, vMax);
-    ierr = DMPlexCreateNumbering_Private(dm, vStart, vEnd, dm->sf, &mesh->globalVertexNumbers);CHKERRQ(ierr);
+    ierr = DMPlexCreateNumbering_Private(dm, vStart, vEnd, 0, NULL, dm->sf, &mesh->globalVertexNumbers);CHKERRQ(ierr);
   }
   *globalVertexNumbers = mesh->globalVertexNumbers;
   PetscFunctionReturn(0);
 }
 
+#undef __FUNCT__
+#define __FUNCT__ "DMPlexCreatePointNumbering"
+PetscErrorCode DMPlexCreatePointNumbering(DM dm, IS *globalPointNumbers)
+{
+  IS             nums[4];
+  PetscInt       depths[4];
+  PetscInt       depth, d, shift = 0;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
+  ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr);
+  depths[0] = depth; depths[1] = 0;
+  for (d = 2; d <= depth; ++d) depths[d] = depth-d+1;
+  for (d = 0; d <= depth; ++d) {
+    PetscInt pStart, pEnd, gsize;
+
+    ierr = DMPlexGetDepthStratum(dm, depths[d], &pStart, &pEnd);CHKERRQ(ierr);
+    ierr = DMPlexCreateNumbering_Private(dm, pStart, pEnd, shift, &gsize, dm->sf, &nums[d]);CHKERRQ(ierr);
+    shift += gsize;
+  }
+  ierr = ISConcatenate(PetscObjectComm((PetscObject) dm), depth+1, nums, globalPointNumbers);
+  for (d = 0; d <= depth; ++d) {ierr = ISDestroy(&nums[d]);CHKERRQ(ierr);}
+  PetscFunctionReturn(0);
+}
+
 
 #undef __FUNCT__
 #define __FUNCT__ "PetscSectionCreateGlobalSectionLabel"

src/dm/impls/plex/plexdistribute.c

   {
     PetscSection originalCoordSection, newCoordSection;
     Vec          originalCoordinates, newCoordinates;
+    PetscInt     bs;
     const char  *name;
 
     ierr = DMGetCoordinateSection(dm, &originalCoordSection);CHKERRQ(ierr);
 
     ierr = DMPlexDistributeField(dm, pointSF, originalCoordSection, originalCoordinates, newCoordSection, newCoordinates);CHKERRQ(ierr);
     ierr = DMSetCoordinatesLocal(*dmParallel, newCoordinates);CHKERRQ(ierr);
+    ierr = VecGetBlockSize(originalCoordinates, &bs);CHKERRQ(ierr);
+    ierr = VecSetBlockSize(newCoordinates, bs);CHKERRQ(ierr);
     ierr = VecDestroy(&newCoordinates);CHKERRQ(ierr);
   }
   /* Distribute labels */

src/dm/impls/plex/plexhdf5.c

+#include <petsc-private/dmpleximpl.h>   /*I      "petscdmplex.h"   I*/
+#include <petsc-private/isimpl.h>
+#include <petscviewerhdf5.h>
+
+PETSC_EXTERN PetscErrorCode VecView_Seq(Vec, PetscViewer);
+PETSC_EXTERN PetscErrorCode VecView_MPI(Vec, PetscViewer);
+
+#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;
+  if (seqnum < 0) PetscFunctionReturn(0);
+  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"
+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 = DMPlexGetFieldType_Internal(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);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "VecView_Plex_HDF5"
+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);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "VecLoad_Plex_HDF5"
+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);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "DMPlexWriteTopology_HDF5_Static"
+static PetscErrorCode DMPlexWriteTopology_HDF5_Static(DM dm, IS globalPointNumbers, PetscViewer viewer)
+{
+  IS              orderIS, conesIS, cellsIS, orntsIS;
+  const PetscInt *gpoint;
+  PetscInt       *order, *sizes, *cones, *ornts;
+  PetscInt        dim, pStart, pEnd, p, conesSize = 0, cellsSize = 0, c = 0, s = 0;
+  PetscErrorCode  ierr;
+
+  PetscFunctionBegin;
+  ierr = ISGetIndices(globalPointNumbers, &gpoint);CHKERRQ(ierr);
+  ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
+  ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr);
+  for (p = pStart; p < pEnd; ++p) {
+    if (gpoint[p] >= 0) {
+      PetscInt coneSize;
+
+      ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr);
+      conesSize += 1;
+      cellsSize += coneSize;
+    }
+  }
+  ierr = PetscMalloc1(conesSize, &order);CHKERRQ(ierr);
+  ierr = PetscMalloc1(conesSize, &sizes);CHKERRQ(ierr);
+  ierr = PetscMalloc1(cellsSize, &cones);CHKERRQ(ierr);
+  ierr = PetscMalloc1(cellsSize, &ornts);CHKERRQ(ierr);
+  for (p = pStart; p < pEnd; ++p) {
+    if (gpoint[p] >= 0) {
+      const PetscInt *cone, *ornt;
+      PetscInt        coneSize, cp;
+
+      ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr);
+      ierr = DMPlexGetCone(dm, p, &cone);CHKERRQ(ierr);
+      ierr = DMPlexGetConeOrientation(dm, p, &ornt);CHKERRQ(ierr);
+      order[s]   = gpoint[p];
+      sizes[s++] = coneSize;
+      for (cp = 0; cp < coneSize; ++cp, ++c) {cones[c] = gpoint[cone[cp]] < 0 ? -(gpoint[cone[cp]]+1) : gpoint[cone[cp]]; ornts[c] = ornt[cp];}
+    }
+  }
+  if (s != conesSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of points %d != %d", s, conesSize);
+  if (c != cellsSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of cone points %d != %d", c, cellsSize);
+  ierr = ISCreateGeneral(PetscObjectComm((PetscObject) dm), conesSize, order, PETSC_OWN_POINTER, &orderIS);CHKERRQ(ierr);
+  ierr = PetscObjectSetName((PetscObject) orderIS, "order");CHKERRQ(ierr);
+  ierr = ISCreateGeneral(PetscObjectComm((PetscObject) dm), conesSize, sizes, PETSC_OWN_POINTER, &conesIS);CHKERRQ(ierr);
+  ierr = PetscObjectSetName((PetscObject) conesIS, "cones");CHKERRQ(ierr);
+  ierr = ISCreateGeneral(PetscObjectComm((PetscObject) dm), cellsSize, cones, PETSC_OWN_POINTER, &cellsIS);CHKERRQ(ierr);
+  ierr = PetscObjectSetName((PetscObject) cellsIS, "cells");CHKERRQ(ierr);
+  ierr = ISCreateGeneral(PetscObjectComm((PetscObject) dm), cellsSize, ornts, PETSC_OWN_POINTER, &orntsIS);CHKERRQ(ierr);
+  ierr = PetscObjectSetName((PetscObject) orntsIS, "orientation");CHKERRQ(ierr);
+  ierr = PetscViewerHDF5PushGroup(viewer, "/topology");CHKERRQ(ierr);
+  ierr = ISView(orderIS, viewer);CHKERRQ(ierr);
+  ierr = ISView(conesIS, viewer);CHKERRQ(ierr);
+  ierr = ISView(cellsIS, viewer);CHKERRQ(ierr);
+  ierr = ISView(orntsIS, viewer);CHKERRQ(ierr);
+  ierr = PetscViewerHDF5PopGroup(viewer);CHKERRQ(ierr);
+  ierr = ISDestroy(&orderIS);CHKERRQ(ierr);
+  ierr = ISDestroy(&conesIS);CHKERRQ(ierr);
+  ierr = ISDestroy(&cellsIS);CHKERRQ(ierr);
+  ierr = ISDestroy(&orntsIS);CHKERRQ(ierr);
+  ierr = ISRestoreIndices(globalPointNumbers, &gpoint);CHKERRQ(ierr);
+
+  ierr = PetscViewerHDF5WriteAttribute(viewer, "/topology/cells", "cell_dim", PETSC_INT, (void *) &dim);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "DMPlexWriteTopology_Vertices_HDF5_Static"
+static PetscErrorCode DMPlexWriteTopology_Vertices_HDF5_Static(DM dm, DMLabel label, PetscInt labelId, PetscViewer viewer)
+{
+  IS              cellIS, globalVertexNumbers;
+  const PetscInt *gvertex;
+  PetscInt       *vertices;
+  PetscInt        dim, depth, vStart, vEnd, v, cellHeight, cStart, cMax, cEnd, cell, conesSize = 0, numCornersLocal = 0, numCorners;
+  hid_t           fileId, groupId;
+  herr_t          status;
+  PetscErrorCode  ierr;
+
+  PetscFunctionBegin;
+  ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
+  ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr);
+  ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
+  ierr = DMPlexGetVTKCellHeight(dm, &cellHeight);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);
+  if (depth == 1) PetscFunctionReturn(0);
+  for (cell = cStart; cell < cEnd; ++cell) {
+    PetscInt *closure = NULL;
+    PetscInt  closureSize, v, Nc = 0;
+
+    if (label) {
+      PetscInt value;
+      ierr = DMLabelGetValue(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 = MPI_Allreduce(&numCornersLocal, &numCorners, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject) dm));CHKERRQ(ierr);
+  if (numCornersLocal && (numCornersLocal != numCorners || numCorners == 1)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Visualization topology currently only supports identical cell shapes");
+
+  ierr = PetscViewerHDF5PushGroup(viewer, "/viz");CHKERRQ(ierr);
+  ierr = PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);CHKERRQ(ierr);
+  status = H5Gclose(groupId);CHKERRQ(status);
+  ierr = PetscViewerHDF5PopGroup(viewer);CHKERRQ(ierr);
+
+  ierr = DMPlexGetVertexNumbering(dm, &globalVertexNumbers);CHKERRQ(ierr);
+  ierr = ISGetIndices(globalVertexNumbers, &gvertex);CHKERRQ(ierr);
+  ierr = PetscMalloc1(conesSize, &vertices);CHKERRQ(ierr);
+  for (cell = cStart, v = 0; cell < cEnd; ++cell) {
+    PetscInt *closure = NULL;
+    PetscInt  closureSize, Nc = 0, p;
+
+    if (label) {
+      PetscInt value;
+      ierr = DMLabelGetValue(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_Internal(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 = ISCreateGeneral(PetscObjectComm((PetscObject) dm), conesSize, vertices, PETSC_OWN_POINTER, &cellIS);CHKERRQ(ierr);
+  ierr = PetscLayoutSetBlockSize(cellIS->map, numCorners);CHKERRQ(ierr);
+  ierr = PetscObjectSetName((PetscObject) cellIS, "cells");CHKERRQ(ierr);
+  ierr = PetscViewerHDF5PushGroup(viewer, "/viz/topology");CHKERRQ(ierr);
+  ierr = ISView(cellIS, viewer);CHKERRQ(ierr);
+#if 0
+  if (numCorners == 1) {
+    ierr = ISView(coneIS, viewer);CHKERRQ(ierr);
+  } else {
+    ierr = PetscViewerHDF5WriteAttribute(viewer, "/viz/topology/cells", "cell_corners", PETSC_INT, (void *) &numCorners);CHKERRQ(ierr);
+  }
+  ierr = ISDestroy(&coneIS);CHKERRQ(ierr);
+#else
+  ierr = PetscViewerHDF5WriteAttribute(viewer, "/viz/topology/cells", "cell_corners", PETSC_INT, (void *) &numCorners);CHKERRQ(ierr);
+#endif
+  ierr = ISDestroy(&cellIS);CHKERRQ(ierr);
+
+  ierr = PetscViewerHDF5WriteAttribute(viewer, "/viz/topology/cells", "cell_dim", PETSC_INT, (void *) &dim);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "DMPlexWriteCoordinates_HDF5_Static"
+static PetscErrorCode DMPlexWriteCoordinates_HDF5_Static(DM dm, PetscViewer viewer)
+{
+  DM             cdm;
+  Vec            coordinates, newcoords;
+  PetscReal      lengthScale;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  ierr = DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);CHKERRQ(ierr);
+  ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr);
+  ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
+  ierr = DMGetLocalVector(cdm, &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 = PetscViewerPushFormat(viewer, PETSC_VIEWER_NATIVE);CHKERRQ(ierr);
+  ierr = VecView(newcoords, viewer);CHKERRQ(ierr);
+  ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);
+  ierr = PetscViewerHDF5PopGroup(viewer);CHKERRQ(ierr);
+  ierr = DMRestoreLocalVector(cdm, &newcoords);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "DMPlexWriteCoordinates_Vertices_HDF5_Static"
+static PetscErrorCode DMPlexWriteCoordinates_Vertices_HDF5_Static(DM dm, PetscViewer viewer)
+{
+  Vec              coordinates, newcoords;
+  PetscSection     cSection;
+  PetscScalar     *coords, *ncoords;
+  const PetscReal *L;
+  PetscReal        lengthScale;
+  PetscInt         vStart, vEnd, v, bs, coordSize, dof, off, d;
+  hid_t            fileId, groupId;
+  herr_t           status;
+  PetscErrorCode   ierr;
+
+  PetscFunctionBegin;
+  ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
+  ierr = DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);CHKERRQ(ierr);
+  ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
+  ierr = VecGetBlockSize(coordinates, &bs);CHKERRQ(ierr);
+  ierr = VecGetLocalSize(coordinates, &coordSize);CHKERRQ(ierr);
+  if (coordSize == (vEnd - vStart)*bs) PetscFunctionReturn(0);
+  ierr = DMGetPeriodicity(dm, NULL, &L);CHKERRQ(ierr);
+  ierr = DMGetCoordinateSection(dm, &cSection);CHKERRQ(ierr);
+  ierr = VecCreate(PetscObjectComm((PetscObject) coordinates), &newcoords);CHKERRQ(ierr);
+  coordSize = 0;
+  for (v = vStart; v < vEnd; ++v) {
+    ierr = PetscSectionGetDof(cSection, v, &dof);CHKERRQ(ierr);
+    if (L && dof == 2) coordSize += dof+1;
+    else               coordSize += dof;
+  }
+  if (L && bs == 2) {ierr = VecSetBlockSize(newcoords, bs+1);CHKERRQ(ierr);}
+  else              {ierr = VecSetBlockSize(newcoords, bs);CHKERRQ(ierr);}
+  ierr = VecSetSizes(newcoords, coordSize, PETSC_DETERMINE);CHKERRQ(ierr);
+  ierr = VecSetType(newcoords,VECSTANDARD);CHKERRQ(ierr);
+  ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr);
+  ierr = VecGetArray(newcoords,   &ncoords);CHKERRQ(ierr);
+  coordSize = 0;
+  for (v = vStart; v < vEnd; ++v) {
+    ierr = PetscSectionGetDof(cSection, v, &dof);CHKERRQ(ierr);
+    ierr = PetscSectionGetOffset(cSection, v, &off);CHKERRQ(ierr);
+    if (L && dof == 2) {
+      ncoords[coordSize++] = cos(2.0*PETSC_PI*coords[off+0]/L[0]);
+      ncoords[coordSize++] = coords[off+1];
+      ncoords[coordSize++] = sin(2.0*PETSC_PI*coords[off+0]/L[0]);
+    } else {
+      for (d = 0; d < dof; ++d, ++coordSize) ncoords[coordSize] = coords[off+d];
+    }
+  }
+  ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr);
+  ierr = VecRestoreArray(newcoords,   &ncoords);CHKERRQ(ierr);
+  ierr = PetscObjectSetName((PetscObject) newcoords, "vertices");CHKERRQ(ierr);
+  ierr = VecScale(newcoords, lengthScale);CHKERRQ(ierr);
+  /* Use the local version to bypass the default group setting */
+  ierr = PetscViewerHDF5PushGroup(viewer, "/viz");CHKERRQ(ierr);
+  ierr = PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);CHKERRQ(ierr);
+  status = H5Gclose(groupId);CHKERRQ(status);
+  ierr = PetscViewerHDF5PopGroup(viewer);CHKERRQ(ierr);
+  ierr = PetscViewerHDF5PushGroup(viewer, "/viz/geometry");CHKERRQ(ierr);
+  ierr = VecView(newcoords, viewer);CHKERRQ(ierr);
+  ierr = PetscViewerHDF5PopGroup(viewer);CHKERRQ(ierr);
+  ierr = VecDestroy(&newcoords);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "DMPlexView_HDF5"
+/* We only write cells and vertices. Does this screw up parallel reading? */
+PetscErrorCode DMPlexView_HDF5(DM dm, PetscViewer viewer)
+{
+  DMLabel           label   = NULL;
+  PetscInt          labelId = 0;
+  IS                globalPointNumbers;
+  const PetscInt   *gpoint;
+  PetscInt          numLabels, l;
+  hid_t             fileId, groupId;
+  herr_t            status;
+  PetscViewerFormat format;
+  PetscErrorCode    ierr;
+
+  PetscFunctionBegin;
+  ierr = DMPlexCreatePointNumbering(dm, &globalPointNumbers);CHKERRQ(ierr);
+  ierr = PetscViewerGetFormat(viewer, &format);CHKERRQ(ierr);
+  ierr = DMPlexWriteCoordinates_HDF5_Static(dm, viewer);CHKERRQ(ierr);
+  if (format == PETSC_VIEWER_HDF5_VIZ) {ierr = DMPlexWriteCoordinates_Vertices_HDF5_Static(dm, viewer);CHKERRQ(ierr);}
+  ierr = DMPlexWriteTopology_HDF5_Static(dm, globalPointNumbers, viewer);CHKERRQ(ierr);
+  if (format == PETSC_VIEWER_HDF5_VIZ) {ierr = DMPlexWriteTopology_Vertices_HDF5_Static(dm, label, labelId, viewer);CHKERRQ(ierr);}
+  /* Write Labels*/
+  ierr = ISGetIndices(globalPointNumbers, &gpoint);CHKERRQ(ierr);
+  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, globalValueIS;
+    const PetscInt *values;
+    PetscInt        numValues, v;
+    PetscBool       isDepth;
+    char            group[PETSC_MAX_PATH_LEN];
+
+    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 = 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);
+    ierr = DMLabelGetValueIS(label, &valueIS);CHKERRQ(ierr);
+    ierr = ISAllGather(valueIS, &globalValueIS);CHKERRQ(ierr);
+    ierr = ISSortRemoveDups(globalValueIS);CHKERRQ(ierr);
+    ierr = ISGetLocalSize(globalValueIS, &numValues);CHKERRQ(ierr);
+    ierr = ISGetIndices(globalValueIS, &values);CHKERRQ(ierr);
+    for (v = 0; v < numValues; ++v) {
+      IS              stratumIS, globalStratumIS;
+      const PetscInt *spoints;
+      PetscInt       *gspoints, n, gn, p;
+      const char     *iname;
+
+      ierr = PetscSNPrintf(group, PETSC_MAX_PATH_LEN, "/labels/%s/%d", name, values[v]);CHKERRQ(ierr);
+      ierr = DMLabelGetStratumIS(label, values[v], &stratumIS);CHKERRQ(ierr);
+
+      ierr = ISGetLocalSize(stratumIS, &n);CHKERRQ(ierr);
+      ierr = ISGetIndices(stratumIS, &spoints);CHKERRQ(ierr);
+      for (gn = 0, p = 0; p < n; ++p) if (gpoint[spoints[p]] >= 0) ++gn;
+      ierr = PetscMalloc1(gn,&gspoints);CHKERRQ(ierr);
+      for (gn = 0, p = 0; p < n; ++p) if (gpoint[spoints[p]] >= 0) gspoints[gn++] = gpoint[spoints[p]];
+      ierr = ISRestoreIndices(stratumIS, &spoints);CHKERRQ(ierr);
+      ierr = ISCreateGeneral(PetscObjectComm((PetscObject) dm), gn, gspoints, PETSC_OWN_POINTER, &globalStratumIS);CHKERRQ(ierr);
+      ierr = PetscObjectGetName((PetscObject) stratumIS, &iname);CHKERRQ(ierr);
+      ierr = PetscObjectSetName((PetscObject) globalStratumIS, iname);CHKERRQ(ierr);
+
+      ierr = PetscViewerHDF5PushGroup(viewer, group);CHKERRQ(ierr);
+      ierr = ISView(globalStratumIS, viewer);CHKERRQ(ierr);
+      ierr = PetscViewerHDF5PopGroup(viewer);CHKERRQ(ierr);
+      ierr = ISDestroy(&globalStratumIS);CHKERRQ(ierr);
+      ierr = ISDestroy(&stratumIS);CHKERRQ(ierr);
+    }
+    ierr = ISRestoreIndices(globalValueIS, &values);CHKERRQ(ierr);
+    ierr = ISDestroy(&globalValueIS);CHKERRQ(ierr);
+    ierr = ISDestroy(&valueIS);CHKERRQ(ierr);
+  }
+  ierr = ISRestoreIndices(globalPointNumbers, &gpoint);CHKERRQ(ierr);
+  ierr = ISDestroy(&globalPointNumbers);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+typedef struct {
+  PetscMPIInt rank;
+  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);
+  {
+    /* Force serial load */
+    ierr = PetscViewerHDF5ReadSizes(viewer, "indices", NULL, &N);CHKERRQ(ierr);
+    ierr = PetscLayoutSetLocalSize(stratumIS->map, !((LabelCtx *) op_data)->rank ? N : 0);CHKERRQ(ierr);
+    ierr = PetscLayoutSetSize(stratumIS->map, N);CHKERRQ(ierr);
+  }
+  ierr = ISLoad(stratumIS, viewer);
+  ierr = PetscViewerHDF5PopGroup(viewer);CHKERRQ(ierr);
+  ierr = ISGetLocalSize(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
+*/
+PetscErrorCode DMPlexLoad_HDF5(DM dm, PetscViewer viewer)
+{
+  LabelCtx        ctx;
+  PetscSection    coordSection;
+  Vec             coordinates;
+  IS              orderIS, conesIS, cellsIS, orntsIS;
+  const PetscInt *order, *cones, *cells, *ornts;
+  PetscReal       lengthScale;
+  PetscInt       *cone, *ornt;
+  PetscInt        dim, spatialDim, N, numVertices, vStart, vEnd, v, pEnd, p, q, maxConeSize = 0, c;
+  hid_t           fileId, groupId;
+  hsize_t         idx = 0;
+  herr_t          status;
+  PetscMPIInt     rank;
+  PetscErrorCode  ierr;
+
+  PetscFunctionBegin;
+  ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRQ(ierr);
+  /* Read toplogy */
+  ierr = PetscViewerHDF5ReadAttribute(viewer, "/topology/cells", "cell_dim", PETSC_INT, (void *) &dim);CHKERRQ(ierr);
+  ierr = DMPlexSetDimension(dm, dim);CHKERRQ(ierr);
+  ierr = PetscViewerHDF5PushGroup(viewer, "/topology");CHKERRQ(ierr);
+
+  ierr = ISCreate(PetscObjectComm((PetscObject) dm), &orderIS);CHKERRQ(ierr);
+  ierr = PetscObjectSetName((PetscObject) orderIS, "order");CHKERRQ(ierr);
+  ierr = ISCreate(PetscObjectComm((PetscObject) dm), &conesIS);CHKERRQ(ierr);
+  ierr = PetscObjectSetName((PetscObject) conesIS, "cones");CHKERRQ(ierr);
+  ierr = ISCreate(PetscObjectComm((PetscObject) dm), &cellsIS);CHKERRQ(ierr);
+  ierr = PetscObjectSetName((PetscObject) cellsIS, "cells");CHKERRQ(ierr);
+  ierr = ISCreate(PetscObjectComm((PetscObject) dm), &orntsIS);CHKERRQ(ierr);
+  ierr = PetscObjectSetName((PetscObject) orntsIS, "orientation");CHKERRQ(ierr);
+  {
+    /* Force serial load */
+    ierr = PetscViewerHDF5ReadSizes(viewer, "order", NULL, &pEnd);CHKERRQ(ierr);
+    ierr = PetscLayoutSetLocalSize(orderIS->map, !rank ? pEnd : 0);CHKERRQ(ierr);
+    ierr = PetscLayoutSetSize(orderIS->map, pEnd);CHKERRQ(ierr);
+    ierr = PetscViewerHDF5ReadSizes(viewer, "cones", NULL, &pEnd);CHKERRQ(ierr);
+    ierr = PetscLayoutSetLocalSize(conesIS->map, !rank ? pEnd : 0);CHKERRQ(ierr);
+    ierr = PetscLayoutSetSize(conesIS->map, pEnd);CHKERRQ(ierr);
+    ierr = PetscViewerHDF5ReadSizes(viewer, "cells", NULL, &N);CHKERRQ(ierr);
+    ierr = PetscLayoutSetLocalSize(cellsIS->map, !rank ? N : 0);CHKERRQ(ierr);
+    ierr = PetscLayoutSetSize(cellsIS->map, N);CHKERRQ(ierr);
+    ierr = PetscViewerHDF5ReadSizes(viewer, "orientation", NULL, &N);CHKERRQ(ierr);
+    ierr = PetscLayoutSetLocalSize(orntsIS->map, !rank ? N : 0);CHKERRQ(ierr);
+    ierr = PetscLayoutSetSize(orntsIS->map, N);CHKERRQ(ierr);
+  }
+  ierr = ISLoad(orderIS, viewer);CHKERRQ(ierr);
+  ierr = ISLoad(conesIS, viewer);CHKERRQ(ierr);
+  ierr = ISLoad(cellsIS, viewer);CHKERRQ(ierr);
+  ierr = ISLoad(orntsIS, viewer);CHKERRQ(ierr);
+  ierr = PetscViewerHDF5PopGroup(viewer);CHKERRQ(ierr);
+  /* Read geometry */
+  ierr = PetscViewerHDF5PushGroup(viewer, "/geometry");CHKERRQ(ierr);
+  ierr = VecCreate(PetscObjectComm((PetscObject) dm), &coordinates);CHKERRQ(ierr);
+  ierr = PetscObjectSetName((PetscObject) coordinates, "vertices");CHKERRQ(ierr);
+  {
+    /* Force serial load */
+    ierr = PetscViewerHDF5ReadSizes(viewer, "vertices", &spatialDim, &N);CHKERRQ(ierr);
+    ierr = VecSetSizes(coordinates, !rank ? N : 0, N);CHKERRQ(ierr);
+    ierr = VecSetBlockSize(coordinates, spatialDim);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 = VecGetLocalSize(coordinates, &numVertices);CHKERRQ(ierr);
+  ierr = VecGetBlockSize(coordinates, &spatialDim);CHKERRQ(ierr);
+  numVertices /= spatialDim;
+  /* Create Plex */
+  ierr = DMPlexSetChart(dm, 0, pEnd);CHKERRQ(ierr);
+  ierr = ISGetIndices(orderIS, &order);CHKERRQ(ierr);
+  ierr = ISGetIndices(conesIS, &cones);CHKERRQ(ierr);
+  for (p = 0; p < pEnd; ++p) {
+    ierr = DMPlexSetConeSize(dm, order[p], cones[p]);CHKERRQ(ierr);
+    maxConeSize = PetscMax(maxConeSize, cones[p]);
+  }
+  ierr = DMSetUp(dm);CHKERRQ(ierr);
+  ierr = ISGetIndices(cellsIS, &cells);CHKERRQ(ierr);
+  ierr = ISGetIndices(orntsIS, &ornts);CHKERRQ(ierr);
+  ierr = PetscMalloc2(maxConeSize,&cone,maxConeSize,&ornt);CHKERRQ(ierr);
+  for (p = 0, q = 0; p < pEnd; ++p) {
+    for (c = 0; c < cones[p]; ++c, ++q) {cone[c] = cells[q]; ornt[c] = ornts[q];}
+    ierr = DMPlexSetCone(dm, order[p], cone);CHKERRQ(ierr);
+    ierr = DMPlexSetConeOrientation(dm, order[p], ornt);CHKERRQ(ierr);
+  }
+  ierr = PetscFree2(cone,ornt);CHKERRQ(ierr);
+  ierr = ISRestoreIndices(orderIS, &order);CHKERRQ(ierr);
+  ierr = ISRestoreIndices(conesIS, &cones);CHKERRQ(ierr);
+  ierr = ISRestoreIndices(cellsIS, &cells);CHKERRQ(ierr);
+  ierr = ISRestoreIndices(orntsIS, &ornts);CHKERRQ(ierr);
+  ierr = ISDestroy(&orderIS);CHKERRQ(ierr);
+  ierr = ISDestroy(&conesIS);CHKERRQ(ierr);
+  ierr = ISDestroy(&cellsIS);CHKERRQ(ierr);
+  ierr = ISDestroy(&orntsIS);CHKERRQ(ierr);
+  ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr);
+  ierr = DMPlexStratify(dm);CHKERRQ(ierr);
+  /* Create coordinates */
+  ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
+  if (numVertices != vEnd - vStart) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of coordinates loaded %d does not match number of vertices %d", numVertices, vEnd - vStart);
+  ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
+  ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr);
+  ierr = PetscSectionSetFieldComponents(coordSection, 0, spatialDim);CHKERRQ(ierr);
+  ierr = PetscSectionSetChart(coordSection, vStart, vEnd);CHKERRQ(ierr);
+  for (v = vStart; v < vEnd; ++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.rank   = rank;
+  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

src/sys/objects/ftn-custom/zoptionsf.c

 #define petscoptionsclear_                 PETSCOPTIONSCLEAR
 #define petscoptionsinsertstring_          PETSCOPTIONSINSERTSTRING
 #define petscoptionsview_                  PETSCOPTIONSVIEW
+#define petscobjectviewfromoptions_        PETSCOBJECTVIEWFROMOPTIONS
 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE)
 #define petscoptionsgetenumprivate_        petscoptionsgetenumprivate
 #define petscoptionsgetbool_               petscoptionsgetbool
 #define petscoptionsclear_                 petscoptionsclear
 #define petscoptionsinsertstring_          petscoptionsinsertstring
 #define petscoptionsview_                  petscoptionsview
+#define petscobjectviewfromoptions_        petscobjectviewfromoptions
 #endif
 
 /* ---------------------------------------------------------------------*/
   *ierr = PetscOptionsView(v);
 }
 
+PETSC_EXTERN void PETSC_STDCALL petscobjectviewfromoptions_(PetscObject *obj,CHAR prefix PETSC_MIXED_LEN(lprefix),CHAR option PETSC_MIXED_LEN(loption),PetscErrorCode *ierr PETSC_END_LEN(lprefix) PETSC_END_LEN(loption))
+{
+  char *p, *o;
+
+  FIXCHAR(prefix, lprefix, p);
+  FIXCHAR(option, loption, o);
+  *ierr = PetscObjectViewFromOptions(*obj, p, o);
+  FREECHAR(prefix, p);
+  FREECHAR(option, o);
+}

src/sys/objects/options.c

 
 #undef __FUNCT__
 #define __FUNCT__ "PetscObjectViewFromOptions"
-/*
+/*@C
   PetscObjectViewFromOptions - Processes command line options to determine if/how a PetscObject is to be viewed. 
 
   Collective on PetscObject
 
   Level: intermediate
 
-*/
+@*/
 PetscErrorCode PetscObjectViewFromOptions(PetscObject obj,const char prefix[],const char optionname[])
 {
   PetscErrorCode    ierr;

src/vec/is/is/impls/block/block.c

 }
 
 #undef __FUNCT__
+#define __FUNCT__ "ISSortRemoveDups_Block"
+PetscErrorCode ISSortRemoveDups_Block(IS is)
+{
+  IS_Block       *sub = (IS_Block*)is->data;
+  PetscInt       bs, n, nb;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  if (sub->sorted) PetscFunctionReturn(0);
+  ierr = PetscLayoutGetBlockSize(is->map, &bs);CHKERRQ(ierr);
+  ierr = PetscLayoutGetLocalSize(is->map, &n);CHKERRQ(ierr);
+  nb   = n/bs;
+  ierr = PetscSortRemoveDupsInt(&nb,sub->idx);CHKERRQ(ierr);
+  ierr = PetscLayoutSetLocalSize(is->map, nb*bs);CHKERRQ(ierr);
+  sub->sorted = PETSC_TRUE;
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
 #define __FUNCT__ "ISSorted_Block"
 PetscErrorCode ISSorted_Block(IS is,PetscBool  *flg)
 {
                                ISRestoreIndices_Block,
                                ISInvertPermutation_Block,
                                ISSort_Block,
+                               ISSortRemoveDups_Block,
                                ISSorted_Block,
                                ISDuplicate_Block,
                                ISDestroy_Block,

src/vec/is/is/impls/general/general.c

 }
 
 #undef __FUNCT__
+#define __FUNCT__ "ISSortRemoveDups_General"
+PetscErrorCode ISSortRemoveDups_General(IS is)
+{
+  IS_General     *sub = (IS_General*)is->data;
+  PetscInt       n;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  if (sub->sorted) PetscFunctionReturn(0);
+  ierr = PetscLayoutGetLocalSize(is->map, &n);CHKERRQ(ierr);
+  ierr = PetscSortRemoveDupsInt(&n,sub->idx);CHKERRQ(ierr);
+  ierr = PetscLayoutSetLocalSize(is->map, n);CHKERRQ(ierr);
+  sub->sorted = PETSC_TRUE;
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
 #define __FUNCT__ "ISSorted_General"
 PetscErrorCode ISSorted_General(IS is,PetscBool  *flg)
 {
                                ISRestoreIndices_General,
                                ISInvertPermutation_General,
                                ISSort_General,
+                               ISSortRemoveDups_General,
                                ISSorted_General,
                                ISDuplicate_General,
                                ISDestroy_General,

src/vec/is/is/impls/stride/stride.c

                                ISRestoreIndices_Stride,
                                ISInvertPermutation_Stride,
                                ISSort_Stride,
+                               ISSort_Stride,
                                ISSorted_Stride,
                                ISDuplicate_Stride,
                                ISDestroy_Stride,

src/vec/is/is/interface/index.c

    Concepts: index sets^sorting
    Concepts: sorting^index set
 
-.seealso: ISSorted()
+.seealso: ISSortRemoveDups(), ISSorted()
 @*/
 PetscErrorCode  ISSort(IS is)
 {
 }
 
 #undef __FUNCT__
+#define __FUNCT__ "ISSortRemoveDups"
+/*@
+  ISSortRemoveDups - Sorts the indices of an index set, removing duplicates.
+
+  Collective on IS
+
+  Input Parameters:
+. is - the index set
+
+  Level: intermediate
+
+  Concepts: index sets^sorting
+  Concepts: sorting^index set
+
+.seealso: ISSort(), ISSorted()
+@*/
+PetscErrorCode ISSortRemoveDups(IS is)
+{
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(is,IS_CLASSID,1);
+  ierr = (*is->ops->sortremovedups)(is);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
 #define __FUNCT__ "ISToGeneral"
 /*@
    ISToGeneral - Converts an IS object of any type to ISGENERAL type
 
    Level: intermediate
 
-.seealso: ISSort()
+.seealso: ISSort(), ISSortRemoveDups()
 @*/
 PetscErrorCode  ISSorted(IS is,PetscBool  *flg)
 {

src/vec/is/utils/isio.c

   hsize_t         rdim, dim;
   hsize_t         dims[3], count[3], offset[3];
   herr_t          status;
-  PetscInt        n, N, bs = 1, bsInd, lenInd, low, timestep;
+  PetscInt        n, N, bs, bsInd, lenInd, low, timestep;
   const PetscInt *ind;
   const char     *isname;
   PetscErrorCode  ierr;
   dim = 0;
   if (timestep >= 0) ++dim;
   ++dim;
-  if (bs >= 1) ++dim;
+  ++dim;
   rdim = H5Sget_simple_extent_dims(filespace, dims, NULL);
   bsInd = rdim-1;
   lenInd = timestep >= 0 ? 1 : 0;
-  if (rdim != dim) {
-    if (rdim == dim+1 && bs == -1) bs = dims[bsInd];
-    else SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Dimension of array in file %d not %d as expected",rdim,dim);
-  } else if (bs >= 1 && bs != (PetscInt) dims[bsInd]) {
-    ierr = ISSetBlockSize(is, dims[bsInd]);CHKERRQ(ierr);
+  if (rdim != dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Dimension of array in file %d not %d as expected",rdim,dim);
+  else if (bs != (PetscInt) dims[bsInd]) {
+    ierr = ISSetBlockSize(is, dims[bsInd]);
     if (ierr) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Block size %d specified for IS does not match blocksize in file %d",bs,dims[bsInd]);
     bs = dims[bsInd];
   }

src/vec/vec/utils/vecio.c

 }
 
 #undef __FUNCT__
+#define __FUNCT__ "PetscViewerHDF5ReadSizes"
+PetscErrorCode PetscViewerHDF5ReadSizes(PetscViewer viewer, const char name[], PetscInt *bs, PetscInt *N)
+{
+  hid_t          file_id, group, dset_id, filespace;
+  hsize_t        rdim, dim;
+  hsize_t        dims[4];
+  herr_t         status;
+  PetscInt       bsInd, lenInd, timestep;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  ierr = PetscViewerHDF5OpenGroup(viewer, &file_id, &group);CHKERRQ(ierr);
+  ierr = PetscViewerHDF5GetTimestep(viewer, &timestep);CHKERRQ(ierr);
+#if (H5_VERS_MAJOR * 10000 + H5_VERS_MINOR * 100 + H5_VERS_RELEASE >= 10800)
+  dset_id = H5Dopen2(group, name, H5P_DEFAULT);