diff --git a/proteus/Archiver.py b/proteus/Archiver.py
index 697256b6d9..74509a4c4a 100644
--- a/proteus/Archiver.py
+++ b/proteus/Archiver.py
@@ -10,25 +10,9 @@
 import numpy
 import os
 import h5py
-from xml.etree.ElementTree import *
 
 memory = Profiling.memory
 
-def indentXML(elem, level=0):
-    i = "\n" + level*"  "
-    if len(elem):
-        if not elem.text or not elem.text.strip():
-            elem.text = i + "  "
-        for e in elem:
-            indentXML(e, level+1)
-            if not e.tail or not e.tail.strip():
-                e.tail = i + "  "
-        if not e.tail or not e.tail.strip():
-            e.tail = i
-    else:
-        if level and (not elem.tail or not elem.tail.strip()):
-            elem.tail = i
-
 #: The dtype of the index arrays the XDMF writers number themselves -- the
 #: point numbering of quadrature, DG and particle meshes -- as opposed to a
 #: mesh's own connectivity, whose dtype comes from the mesh. The writers
@@ -55,6 +39,104 @@ def xdmf_number_type(*arrays):
     item = data_item_for(numpy.empty(0, numpy.result_type(*arrays)), inline=True)
     return {"DataType": item["data_type"], "Precision": str(item["precision"])}
 
+
+class ArchiveCollection(object):
+    """One temporal collection of the running archive.
+
+    A collection is one finite-element space's sequence of grids, one per
+    archived instant, and it holds only the *current* step's grid: the
+    writers build it, :meth:`AR_base.sync` stores it in the HDF5 file and
+    clears it, so the running document is never more than one step deep.
+    The grid is a ymf grid dict (``topology``, ``geometry``, ``attributes``
+    and an optional ``name``), built by :meth:`AR_base.write_grid`,
+    :meth:`AR_base.write_topology`, :meth:`AR_base.write_geometry` and
+    :meth:`AR_base.write_field`, and stored as it is.
+
+    Before the XDMF retirement this was a ``<Grid GridType="Collection"
+    CollectionType="Temporal">`` element whose last child was the current
+    step; the writers compared that child's ``<Time>`` with the time being
+    written to decide whether to start a new grid. :meth:`is_current` is
+    that test.
+    """
+    def __init__(self, name):
+        self.name = name
+        self.grid = None
+        self.time = None
+        self.tCount = None
+
+    def new_grid(self, t, tCount, name=None):
+        """Start the grid of a new step at time ``t``. Returns it."""
+        grid = {}
+        if name is not None:
+            grid["name"] = name
+        grid["topology"] = None
+        grid["geometry"] = None
+        grid["attributes"] = []
+        self.grid, self.time, self.tCount = grid, float(t), tCount
+        return grid
+
+    def is_current(self, t):
+        """Whether the grid held is the one for time ``t``.
+
+        False after :meth:`clear`, so the first writer of a step builds
+        the grid and the others at the same time attach fields to it.
+        """
+        return self.grid is not None and self.time == float(t)
+
+    def clear(self):
+        """Drop the stored step's grid; the next writer starts a new one."""
+        self.grid = None
+
+
+class ArchiveDomain(object):
+    """The running archive's temporal collections, in the order created.
+
+    What the writers used to reach as ``ar.domain``, the ``<Domain>``
+    element of the in-memory XDMF tree. Collections are found by name
+    (:meth:`collection`) or by position (``domain[0]``, as
+    ``MeshTools.Mesh.writeMeshXdmf`` does for the base mesh and its
+    element-boundary mesh), and are never removed: the order they were
+    created in is the order the archive records them.
+    """
+    def __init__(self):
+        self.collections = []
+
+    def __iter__(self):
+        return iter(self.collections)
+
+    def __len__(self):
+        return len(self.collections)
+
+    def __getitem__(self, index):
+        return self.collections[index]
+
+    def names(self):
+        return [c.name for c in self.collections]
+
+    def get(self, name):
+        """The collection called ``name``, or None."""
+        for c in self.collections:
+            if c.name == name:
+                return c
+        return None
+
+    def collection(self, name, create=True):
+        """The collection called ``name``, created at the end if missing.
+
+        With ``create=False`` a missing collection is None rather than
+        made: a writer continuing a run (``init=False``) expects its
+        collection to exist.
+        """
+        c = self.get(name)
+        if c is None and create:
+            c = ArchiveCollection(name)
+            self.collections.append(c)
+        return c
+
+    def clear_step(self):
+        for c in self.collections:
+            c.clear()
+
 class ArchiveFlags(object):
     EVERY_MODEL_STEP     = 0
     EVERY_USER_STEP      = 1
@@ -70,10 +152,14 @@ class AR_base(object):
     <filename>.ymf`` derives an ``.xmf`` that references the same HDF5
     datasets; ``gatherTimes`` rebuilds both from the ``.h5`` alone.
 
-    ``self.tree`` is in memory only. The writers still build XDMF elements
-    under ``self.domain`` (Phase 1 of the YMF campaign converts them to
-    dicts), and :meth:`allGatherIncremental` turns each step's elements into
-    ymf grids; nothing serializes the tree.
+    ``self.domain`` is the running document: an :class:`ArchiveDomain` of
+    named collections, each holding the ymf grid dict of the current step.
+    The writers build those dicts through :meth:`write_grid`,
+    :meth:`write_topology`, :meth:`write_geometry` and :meth:`write_field`,
+    and :meth:`allGatherIncremental` stores them. Nothing in proteus builds
+    or parses XML any more: ``ymf2xmf`` derives XDMF from the ``.ymf``, and
+    the one XDMF proteus still reads -- the ``.xmf`` beside a <= 1.9.x
+    archive, for the order of its collections -- goes through ``ymf.xdmf``.
 
     ``useGlobalXMF`` and ``gatherAtClose`` are accepted and ignored. They
     chose between one ``.xmf`` and one per rank, gathered at close; the
@@ -101,11 +187,7 @@ def __init__(self,dataDir,filename,
         self._liveYmfWarned = False
         self.global_sync = global_sync
         comm_world = self.comm.comm.tompi4py()
-        self.tree=ElementTree(
-            Element("Xdmf",
-                    {"Version":"2.0",
-                     "xmlns:xi":"http://www.w3.org/2001/XInclude"})
-        )
+        self.domain = ArchiveDomain()
         if not useTextArchive:
             if hotStart:
                 mode = "a"
@@ -135,21 +217,18 @@ def __init__(self,dataDir,filename,
         if hotStart and self.hdfFile is not None:
             self._seed_domain()
     def _seed_domain(self):
-        """Give a hot-started run the ``<Domain>`` its writers append to.
+        """Give a hot-started run the collections its writers append to.
 
         The writers find their temporal collection by name, or as the
-        first ``<Grid>`` under the domain (``MeshTools.writeMeshXdmf``), so
-        a continuing run needs the collections the archive already has, in
+        first one in the domain (``MeshTools.writeMeshXdmf``), so a
+        continuing run needs the collections the archive already has, in
         the order they were written, and nothing else: the grids of earlier
-        steps were never kept in the tree, since :meth:`sync` deletes them
-        after storing each step. This used to come from parsing the
-        ``.xmf``; the HDF5 file records the names and their order itself.
+        steps are stored by :meth:`sync` and not kept. This used to come
+        from parsing the ``.xmf``; the HDF5 file records the names and
+        their order itself.
         """
-        domain = SubElement(self.tree.getroot(), "Domain")
         for name in self._collection_names():
-            SubElement(domain, "Grid", {"Name":name,
-                                        "GridType":"Collection",
-                                        "CollectionType":"Temporal"})
+            self.domain.collection(name)
     def _collection_names(self):
         """The archive's temporal-collection names, in the order written.
 
@@ -171,9 +250,8 @@ def _collection_names(self):
                 "%s was written by proteus <= 1.9.x, which did not record "
                 "the order of its time collections; hot starting from it "
                 "needs the %s written beside it" % (self.hdfFilename, xmf))
-        Domain = ElementTree(file=xmf).getroot().find("Domain")
-        return [collection.attrib["Name"]
-                for collection in Domain.findall("Grid")]
+        from ymf.xdmf import collection_names
+        return collection_names(xmf)
     #: Bumped when the in-archive metadata layout changes incompatibly.
     #: Written as an attribute on the HDF5 root by
     #: :meth:`allGatherIncremental` and checked by
@@ -372,19 +450,18 @@ def _v1_grids(self, grid_array):
           fragment that still has one was written globally. Shape alone is
           not enough -- both layouts give shape (1,) at a single rank.
         """
-        from xml.etree.ElementTree import fromstring
-
-        from ymf.xdmf import parse_grid_element
+        from ymf.xdmf import read_grid_fragment
 
-        elements = [fromstring(grid_array[j])
-                    for j in range(grid_array.shape[0])]
-        first_time = elements[0].find("Time")
+        parsed = [read_grid_fragment(grid_array[j])
+                  for j in range(grid_array.shape[0])]
+        grids = [grid for grid, _ in parsed]
+        first_time = parsed[0][1]
         global_sync = first_time is not None
         if first_time is not None:
-            t = float(first_time.attrib["Value"])
+            t = first_time
         else:
             t = float(grid_array.attrs["Time"])
-        return t, [parse_grid_element(e) for e in elements], global_sync
+        return t, grids, global_sync
 
     def load_archived_domain(self):
         """The archive's domain, assembling it from HDF5 on first use.
@@ -571,11 +648,10 @@ def allGatherIncremental(self):
         """
         from mpi4py import MPI
         from ymf.archive import dump_grid
-        from ymf.xdmf import parse_grid_element
 
         logEvent("Gathering Archive Time step")
         self.comm.barrier()
-        Domain = self.tree.getroot().find("Domain")
+        Domain = self.domain
 
         comm_world = self.comm.comm.tompi4py()
         if self.hdfFile is not None:
@@ -588,16 +664,25 @@ def allGatherIncremental(self):
             #the file's keys is guessy where this is not. Written by every
             #rank with the same value, since attribute writes are collective.
             self.hdfFile.attrs[self.COLLECTIONS_ATTR] = "\n".join(
-                c.attrib['Name'] for c in Domain)
+                c.name for c in Domain)
             self.hdfFile.attrs[self.GLOBAL_SYNC_ATTR] = \
                 1 if self.global_sync else 0
 
         live_steps = []
         for i, TemporalGridCollection in enumerate(Domain):
-            GridLocal = TemporalGridCollection[-1]
-            time_elem = GridLocal.find("Time")
-            TimeAttrib = time_elem.attrib['Value'] if time_elem is not None else "0.0"
-            local_grid = parse_grid_element(GridLocal)
+            local_grid = TemporalGridCollection.grid
+            if local_grid is None:
+                raise RuntimeError(
+                    "archive collection %r has no grid for this step: every "
+                    "collection is written at every archived instant"
+                    % (TemporalGridCollection.name,))
+            if local_grid["topology"] is None or local_grid["geometry"] is None:
+                raise RuntimeError(
+                    "archive collection %r: the grid for t=%r has no %s"
+                    % (TemporalGridCollection.name, TemporalGridCollection.time,
+                       "topology" if local_grid["topology"] is None else "geometry"))
+            t = TemporalGridCollection.time
+            TimeAttrib = repr(t)
 
             if not self.global_sync:
                 grid_dicts = comm_world.gather(local_grid)
@@ -616,7 +701,7 @@ def allGatherIncremental(self):
             width = max(int(width), 1)
 
             dataset_name = self._metadata_dataset_name(
-                TemporalGridCollection.attrib['Name'], self.n_datasets)
+                TemporalGridCollection.name, self.n_datasets)
             if self.hdfFile is not None:
                 try:
                     grid_data = self.hdfFile.create_dataset(
@@ -634,8 +719,8 @@ def allGatherIncremental(self):
                 #document costs no extra read and no extra parse
                 live_grids = (grid_dicts if not self.global_sync
                               else [local_grid])
-                live_steps.append((i, TemporalGridCollection.attrib['Name'],
-                                   float(TimeAttrib), live_grids))
+                live_steps.append((i, TemporalGridCollection.name,
+                                   t, live_grids))
         self.n_datasets += 1
         if live_steps and self.writeIncrementalYmf:
             self._appendLiveStep(live_steps)
@@ -644,9 +729,8 @@ def sync(self):
         logEvent("Syncing Archive",level=3)
         memory()
         self.allGatherIncremental()
-        #the step is stored; drop its grids so the tree stays one step deep
-        for TemporalGridCollection in self.tree.getroot().find("Domain"):
-            del TemporalGridCollection[:]
+        #the step is stored; drop its grids so the domain stays one step deep
+        self.domain.clear_step()
         if self.hdfFile is not None:
             self.comm.barrier()
             self.hdfFile.flush()
@@ -677,12 +761,12 @@ def create_dataset_sync(self,name,offsets,data):
     def write_field(self, grid, name, data, tCount,
                     center="Node", rank="Scalar", dimensions=None,
                     sync_offsets=None, sync_data=None, dataset=None,
-                    text_stem=None):
-        """Attach one field to ``grid`` and write its array. Returns the DataItem.
+                    text_stem=None, write_data=True):
+        """Attach one field to ``grid`` and write its array. Returns the data item.
 
-        This is the single place an ``Attribute`` + ``DataItem`` pair gets
-        built. Before it existed, that eight-line block was written out by
-        hand at roughly fifty call sites across Archiver.py, FemTools.py and
+        This is the single place an attribute and its data item get built.
+        Before it existed, that eight-line block was written out by hand at
+        roughly fifty call sites across Archiver.py, FemTools.py and
         MeshTools.py, each repeating the same four-way branch on
         ``global_sync`` and ``hdfFile``. Two consequences of that
         duplication worth knowing about:
@@ -693,18 +777,19 @@ def write_field(self, grid, name, data, tCount,
           garbage. Here they come from the array's own dtype via
           :func:`ymf.archive.data_item_for`, so they cannot disagree with
           the data.
-        * The ``xi:include`` for the text fallback had to be attached to the
-          right ``DataItem`` by hand, and at ``Archiver.py:1475`` it wasn't
-          -- the reference for one field landed on another field's
-          DataItem. Here the include is attached to the DataItem this call
-          just created, so that class of bug is not expressible.
+        * The sidecar reference for the text fallback had to be attached
+          to the right DataItem by hand, and at ``Archiver.py:1475`` it
+          wasn't -- the reference for one field landed on another field's
+          DataItem. Here it is part of the data item this call creates,
+          so that class of bug is not expressible.
 
         Parameters
         ----------
-        grid : Element
-            The XDMF ``Grid`` element to attach the Attribute to.
+        grid : dict
+            The ymf grid dict to attach the attribute to, from
+            :meth:`write_grid`.
         name : str
-            Field name, used for both the Attribute and the dataset.
+            Field name, used for both the attribute and the dataset.
         data : numpy.ndarray
             The values to write. In the ``global_sync`` case this is the
             full local array and ``sync_data`` is the owned slice actually
@@ -727,8 +812,12 @@ def write_field(self, grid, name, data, tCount,
         text_stem : str, optional
             Overrides the sidecar filename stem for the text fallback,
             which defaults to ``<name><tCount>``.
+        write_data : bool
+            False attaches the attribute without writing its array, for a
+            field whose array was written at an earlier step and is
+            referenced again (the mesh writers' local-to-global maps).
         """
-        from ymf.archive import data_item_for
+        from ymf.archive import attribute, data_item_for
 
         if dataset is not None:
             dataset_name = dataset
@@ -764,20 +853,12 @@ def write_field(self, grid, name, data, tCount,
             #disabling keeps that safety net for the common serial run.
             check=not (self.global_sync and self.size > 1),
         )
+        grid["attributes"].append(
+            attribute(name, item, center=center, attribute_type=rank))
 
-        attribute = SubElement(grid, "Attribute",
-                               {"Name": name,
-                                "AttributeType": rank,
-                                "Center": center})
-        values = SubElement(attribute, "DataItem",
-                            {"Format": self.dataItemFormat,
-                             "DataType": item["data_type"],
-                             "Precision": str(item["precision"]),
-                             "Dimensions": " ".join(
-                                 str(d) for d in item["dimensions"])})
-
+        if not write_data:
+            return item
         if self.hdfFile is not None:
-            values.text = item["data"]
             if self.global_sync:
                 self.create_dataset_sync(dataset_name,
                                          offsets=sync_offsets,
@@ -788,30 +869,22 @@ def write_field(self, grid, name, data, tCount,
             assert not self.global_sync, \
                 "global_sync is not supported with text heavy data"
             numpy.savetxt(text_path, data)
-            # Attached to the DataItem this call created -- see the note above.
-            SubElement(values, "xi:include",
-                       {"parse": "text", "href": item["include"]})
-        return values
+        return item
 
     def write_grid(self, collection, grid_name, t, tCount):
-        """Create the ``Grid``/``Time`` pair for one timestep of a mesh.
+        """Start the grid of one timestep of a mesh. Returns the grid dict.
 
-        Returns ``(grid, time)``. The ten ``writeMeshXdmf_*`` methods in
-        :class:`XdmfWriter` each open with this same pair, spelled either
-        ``"%e" % (t,)``/``str(tCount)`` or
-        ``"{0:e}".format(t)``/``"{0:d}".format(tCount)`` depending on the
-        method -- both produce identical strings, so one spelling serves.
+        ``collection`` is the :class:`ArchiveCollection` of the space; the
+        grid becomes its current one, and the collection records the time
+        (what used to be the ``<Time>`` element). The writers then fill
+        it with :meth:`write_topology`, :meth:`write_geometry` and
+        :meth:`write_field`.
         """
-        grid = SubElement(collection, "Grid",
-                          {"Name": grid_name, "GridType": "Uniform"})
-        time = SubElement(grid, "Time",
-                          {"Value": "{0:e}".format(t),
-                           "Name": "{0:d}".format(tCount)})
-        return grid, time
+        return collection.new_grid(t, tCount, grid_name)
 
     def write_topology(self, grid, topology_type, n_elements, dimensions,
                        dataset, text_stem, nodes_per_element=None, *, dtype):
-        """Write a ``Topology`` and its ``DataItem``. Returns the DataItem.
+        """Set the grid's topology, naming its array. Returns the data item.
 
         The array reference is filled in here, but the array itself is
         **not** written: mesh datasets are only created when ``init or
@@ -825,45 +898,43 @@ def write_topology(self, grid, topology_type, n_elements, dimensions,
         :func:`xdmf_number_type`). It has no default, so a caller cannot
         describe an array without saying what it holds.
         """
-        attrs = {"Type": topology_type, "NumberOfElements": str(n_elements)}
-        if nodes_per_element is not None:
-            attrs["NodesPerElement"] = str(nodes_per_element)
-        topology = SubElement(grid, "Topology", attrs)
-        return self._mesh_data_item(topology, dimensions, dataset, text_stem,
-                                    dtype)
+        from ymf.archive import topology
+        item = self._mesh_data_item(dimensions, dataset, text_stem, dtype)
+        grid["topology"] = topology(topology_type, n_elements, item,
+                                    nodes_per_element=nodes_per_element)
+        return item
 
     def write_geometry(self, grid, dimensions, dataset, text_stem,
                        geometry_type="XYZ", *, dtype):
-        """Write a ``Geometry`` and its ``DataItem``. Returns the DataItem.
+        """Set the grid's geometry, naming its array. Returns the data item.
 
         Same division of labour, and the same ``dtype``, as
         :meth:`write_topology`.
         """
-        geometry = SubElement(grid, "Geometry", {"Type": geometry_type})
-        return self._mesh_data_item(geometry, dimensions, dataset, text_stem,
-                                    dtype)
+        from ymf.archive import geometry
+        item = self._mesh_data_item(dimensions, dataset, text_stem, dtype)
+        grid["geometry"] = geometry(item, geometry_type=geometry_type)
+        return item
 
-    def _mesh_data_item(self, parent, dimensions, dataset, text_stem, dtype):
-        """The DataItem shared by :meth:`write_topology`/:meth:`write_geometry`.
+    def _mesh_data_item(self, dimensions, dataset, text_stem, dtype):
+        """The data item shared by :meth:`write_topology`/:meth:`write_geometry`.
 
         Unlike :meth:`write_field` this cannot read the dtype off the array,
         because the array often does not exist yet at this point -- the
         reference is written on every pass and the data only when the mesh
         changed. So the caller names what the array is made from.
         """
+        from ymf.archive import data_item
         dtypes = dtype if isinstance(dtype, tuple) else (dtype,)
-        attrs = {"Format": self.dataItemFormat,
-                 **xdmf_number_type(*dtypes),
-                 "Dimensions": " ".join(str(d) for d in dimensions)}
-        item = SubElement(parent, "DataItem", attrs)
+        kind = xdmf_number_type(*dtypes)
         if self.hdfFile is not None:
-            item.text = "{0:s}:/{1:s}".format(self.hdfFilename, dataset)
+            reference = {"data": "{0:s}:/{1:s}".format(self.hdfFilename, dataset)}
         else:
-            SubElement(item, "xi:include",
-                       {"parse": "text",
-                        "href": "./{0:s}/{1:s}.txt".format(self.textDataDir,
-                                                           text_stem)})
-        return item
+            reference = {"include": "./{0:s}/{1:s}.txt".format(self.textDataDir,
+                                                               text_stem)}
+        return data_item([int(d) for d in dimensions],
+                         data_type=kind["DataType"],
+                         precision=int(kind["Precision"]), **reference)
 
 def readArchiveDomain(filename, dataDir='.'):
     """The ymf domain of an existing archive, read from its HDF5 file.
@@ -897,67 +968,61 @@ def readArchiveDomain(filename, dataDir='.'):
 XdmfArchive=AR_base
 
 ########################################################################
-#for writing out various quantities in Xdmf format
-#import xml.etree.ElementTree as ElementTree
-#from xml.etree.ElementTree import SubElement
-import numpy
+#writers for the function spaces whose nodes are not the mesh's own
 class XdmfWriter(object):
     """
-    collect functionality for writing data to Xdmf format
+    collect functionality for writing data to the archive
 
     Writer is supposed to keep track of grid collection (temporal collection)
 
     as well as current grid under grid collection where data belong,
     since data are associated with a grid of specific type
     (e.g., P1 Lagrange, P2 Lagrange, elementQuadrature dictionary, ...)
+
+    Each ``writeMeshXdmf_*`` method describes one space's mesh as a ymf
+    grid dict through the archive's builders (``write_grid``,
+    ``write_topology``, ``write_geometry``) and writes the arrays behind it
+    when ``init or meshChanged``; the ``write*Xdmf_*`` methods attach fields
+    to the current grid through ``write_field``. The names survive from
+    when these methods built XDMF elements; what they build now is the
+    archive's own representation, which ``ymf2xmf`` turns into XDMF.
     """
     def __init__(self,shareSingleGrid=True,arGridCollection=None,arGrid=None,arTime=None):
         self.arGridCollection = arGridCollection #collection of "grids" (at least one per time level)
         self.arGrid           = arGrid #grid in collection that data should be associated with
-        self.arTime           = arTime #time level for data
+        self.arTime           = arTime #kept for callers that read it; the collection records the time
         self.shareSingleGrid  = shareSingleGrid
 
     def setGridCollectionAndGridElements(self,init,ar,arGrid,t,spaceSuffix):
+        """Find or create this space's collection; returns the grid name.
+
+        On ``init`` the collection ``Mesh<spaceSuffix>`` is created (or
+        found, if another writer of the same space made it); otherwise it
+        is looked up, and a writer that already holds one keeps it. With
+        ``arGrid`` given -- another writer's grid for this space at this
+        time -- that grid becomes the current one, so a second component
+        attaches its fields to the grid the first one built.
         """
-        attempt at boiler plate code to grab current arGridCollection and grid for
-        a given type and time level t
-        returns gridName to use in writing mesh
+        found = ar.domain.collection("Mesh"+spaceSuffix, create=init)
+        if found is not None:
+            self.arGridCollection = found
+        assert self.arGridCollection is not None, \
+            "no collection Mesh%s in the archive; it is created on init" % (spaceSuffix,)
+        if arGrid is not None:
+            self.arGrid = arGrid
+        return "Grid"+spaceSuffix
+
+    def _needsGrid(self, t):
+        """Whether this space still lacks a grid for time ``t``.
+
+        The first writer of a step builds the grid; the others at the same
+        time reuse it. Replaces the comparison of the current grid's
+        ``<Time>`` with ``"{0:e}".format(t)`` that every writer made.
         """
-        if init:
-            #ar should have domain now and mesh should have gridCollection
-            #but want own grid collection to start
-            self.arGridCollection = SubElement(ar.domain,"Grid",{"Name":"Mesh"+spaceSuffix,
-                                                                 "GridType":"Collection",
-                                                                 "CollectionType":"Temporal"})
-        elif self.arGridCollection is None:#try to get existing grid collection
-            for child in ar.domain:
-                if child.tag == "Grid" and child.attrib["Name"] == "Mesh"+spaceSuffix:
-                    self.arGridCollection = child
-        assert self.arGridCollection is not None
-
-        #see if dgp1 grid exists with current time?
-
-        if self.shareSingleGrid:
-            gridName = "Grid"+spaceSuffix
-            if arGrid is not None:
-                self.arGrid = arGrid
-                gt = arGrid.find("Time")
-                self.arTime = gt
-            #brute force search through child grids of arGridCollection
-            #grids = self.arGridCollection.findall("Grid")
-            #for g in grids:
-            #    for gt in g:
-            #        if gt.tag == "Time" and gt.attrib["Value"] == str(t):
-            #            self.arTime = gt
-            #            self.arGrid = g
-            #            break
-            #end brute force search
-
-        else:
-            #`name` was not in scope here; this branch is only reachable
-            #with shareSingleGrid=False, which no caller sets.
-            gridName = "Grid"+spaceSuffix
-        return gridName
+        if self.arGridCollection.is_current(t):
+            self.arGrid = self.arGridCollection.grid
+            return False
+        return True
 
     def writeMeshXdmf_elementQuadrature(self,ar,mesh,spaceDim,x,t=0.0,
                                        init=False,meshChanged=False,arGrid=None,tCount=0):
@@ -1030,106 +1095,49 @@ def writeMeshXdmf_quadraturePoints(self,ar,mesh,spaceDim,x,t=0.0,quadratureType=
         gridName = self.setGridCollectionAndGridElements(init,ar,arGrid,t,spaceSuffix)
 
         assert len(x.shape) == 3 #make sure have right type of dictionary
-        if self.arGrid is None or self.arTime.get('Value') != "{0:e}".format(t):
+        self.mesh = mesh
+        if self._needsGrid(t):
+            Xdmf_NodesPerElement = x.shape[1]
             if ar.global_sync:
-                self.mesh = mesh
-                Xdmf_ElementTopology = "Polyvertex"
+                assert ar.hdfFile is not None, "global_sync not supported with text  heavy data"
                 Xdmf_NumberOfElements= mesh.globalMesh.nElements_global
-                Xdmf_NodesPerElement = x.shape[1]
-                Xdmf_NodesGlobal     = Xdmf_NumberOfElements*Xdmf_NodesPerElement
-
-                self.arGrid = SubElement(self.arGridCollection,"Grid",{"Name":gridName,"GridType":"Uniform"})
-                self.arTime = SubElement(self.arGrid,"Time",{"Value":"%e" % (t,),"Name":"%i" % (tCount)})
-                topology    = SubElement(self.arGrid,"Topology",
-                                         {"Type":Xdmf_ElementTopology,
-                                          "NumberOfElements":str(Xdmf_NumberOfElements),
-                                          "NodesPerElement":str(Xdmf_NodesPerElement)})
-                elements    = SubElement(topology,"DataItem",
-                                         {"Format":ar.dataItemFormat,
-                                          **xdmf_number_type(GENERATED_INDEX_DTYPE),
-                                          "Dimensions":"%i %i" % (Xdmf_NumberOfElements,Xdmf_NodesPerElement)})
-                geometry    = SubElement(self.arGrid,"Geometry",{"Type":"XYZ"})
-                nodes       = SubElement(geometry,"DataItem",
-                                         {"Format":ar.dataItemFormat,
-                                          **xdmf_number_type(x),
-                                          "Dimensions":"%i %i" % (Xdmf_NodesGlobal,3)})
-                if ar.hdfFile is not None:
-                    elements.text = ar.hdfFilename+":/elements"+spaceSuffix+str(tCount)
-                    nodes.text    = ar.hdfFilename+":/nodes"+spaceSuffix+str(tCount)
-                    if init or meshChanged:
-                        #this will fail if elements_dgp1 already exists
-                        #q_l2g = numpy.zeros((Xdmf_NumberOfElements,Xdmf_NodesPerElement),'i')
-                        #brute force to start
-                        #for eN in range(Xdmf_NumberOfElements):
-                        #    for nN in range(Xdmf_NodesPerElement):
-                        #        q_l2g[eN,nN] = eN*Xdmf_NodesPerElement + nN
-                        #
-                        from proteus import Comm
-                        comm = Comm.get()
-                        q_l2g = numpy.arange(mesh.globalMesh.elementOffsets_subdomain_owned[ar.rank]*Xdmf_NodesPerElement,
-                                             mesh.globalMesh.elementOffsets_subdomain_owned[ar.rank+1]*Xdmf_NodesPerElement,
-                                             dtype=GENERATED_INDEX_DTYPE).reshape((mesh.nElements_owned,Xdmf_NodesPerElement))
-                        ar.create_dataset_sync('elements'+spaceSuffix+str(tCount),
-                                               offsets=mesh.globalMesh.elementOffsets_subdomain_owned,
-                                               data = q_l2g)
-                        ar.create_dataset_sync('nodes'+spaceSuffix+str(tCount),
-                                               offsets=mesh.globalMesh.elementOffsets_subdomain_owned*Xdmf_NodesPerElement,
-                                               data = x[:mesh.nElements_owned].reshape((mesh.nElements_owned*Xdmf_NodesPerElement,3)))
-                else:
-                    assert False, "global_sync not supported with text  heavy data"
+                elements_dataset = 'elements'+spaceSuffix+str(tCount)
+                nodes_dataset    = 'nodes'+spaceSuffix+str(tCount)
             else:
-                Xdmf_ElementTopology = "Polyvertex"
                 Xdmf_NumberOfElements= x.shape[0]
-                Xdmf_NodesPerElement = x.shape[1]
-                Xdmf_NodesGlobal     = Xdmf_NumberOfElements*Xdmf_NodesPerElement
-
-                self.arGrid = SubElement(self.arGridCollection,"Grid",{"Name":gridName,"GridType":"Uniform"})
-                self.arTime = SubElement(self.arGrid,"Time",{"Value":"%e" % (t,),"Name":"%i" % (tCount)})
-                topology    = SubElement(self.arGrid,"Topology",
-                                         {"Type":Xdmf_ElementTopology,
-                                          "NumberOfElements":str(Xdmf_NumberOfElements),
-                                          "NodesPerElement":str(Xdmf_NodesPerElement)})
-                elements    = SubElement(topology,"DataItem",
-                                         {"Format":ar.dataItemFormat,
-                                          **xdmf_number_type(GENERATED_INDEX_DTYPE),
-                                          "Dimensions":"%i %i" % (Xdmf_NumberOfElements,Xdmf_NodesPerElement)})
-                geometry    = SubElement(self.arGrid,"Geometry",{"Type":"XYZ"})
-                nodes       = SubElement(geometry,"DataItem",
-                                         {"Format":ar.dataItemFormat,
-                                          **xdmf_number_type(x),
-                                          "Dimensions":"%i %i" % (Xdmf_NodesGlobal,3)})
-                if ar.hdfFile is not None:
-                    elements.text = ar.hdfFilename+":/elements"+str(ar.rank)+spaceSuffix+str(tCount)
-                    nodes.text    = ar.hdfFilename+":/nodes"+str(ar.rank)+spaceSuffix+str(tCount)
-                    if init or meshChanged:
-                        #this will fail if elements_dgp1 already exists
-                        #q_l2g = numpy.zeros((Xdmf_NumberOfElements,Xdmf_NodesPerElement),'i')
-                        #brute force to start
-                        #for eN in range(Xdmf_NumberOfElements):
-                        #    for nN in range(Xdmf_NodesPerElement):
-                        #        q_l2g[eN,nN] = eN*Xdmf_NodesPerElement + nN
-                        #
-                        q_l2g = numpy.arange(Xdmf_NumberOfElements*Xdmf_NodesPerElement,dtype=GENERATED_INDEX_DTYPE).reshape((Xdmf_NumberOfElements,Xdmf_NodesPerElement))
-                        ar.create_dataset_async('elements'+str(ar.rank)+spaceSuffix+str(tCount), data = q_l2g)
-                        ar.create_dataset_async('nodes'+str(ar.rank)+spaceSuffix+str(tCount), data = x.flat[:])
+                elements_dataset = 'elements'+str(ar.rank)+spaceSuffix+str(tCount)
+                nodes_dataset    = 'nodes'+str(ar.rank)+spaceSuffix+str(tCount)
+            Xdmf_NodesGlobal     = Xdmf_NumberOfElements*Xdmf_NodesPerElement
+            elements_stem = 'elements'+spaceSuffix+str(tCount)
+            nodes_stem    = 'nodes'+spaceSuffix+str(tCount)
+
+            self.arGrid = ar.write_grid(self.arGridCollection, gridName, t, tCount)
+            ar.write_topology(self.arGrid, "Polyvertex", Xdmf_NumberOfElements,
+                              [Xdmf_NumberOfElements, Xdmf_NodesPerElement],
+                              elements_dataset, elements_stem,
+                              nodes_per_element=Xdmf_NodesPerElement,
+                              dtype=GENERATED_INDEX_DTYPE)
+            ar.write_geometry(self.arGrid, [Xdmf_NodesGlobal, 3],
+                              nodes_dataset, nodes_stem, dtype=x)
+            if init or meshChanged:
+                if ar.global_sync:
+                    q_l2g = numpy.arange(mesh.globalMesh.elementOffsets_subdomain_owned[ar.rank]*Xdmf_NodesPerElement,
+                                         mesh.globalMesh.elementOffsets_subdomain_owned[ar.rank+1]*Xdmf_NodesPerElement,
+                                         dtype=GENERATED_INDEX_DTYPE).reshape((mesh.nElements_owned,Xdmf_NodesPerElement))
+                    ar.create_dataset_sync(elements_dataset,
+                                           offsets=mesh.globalMesh.elementOffsets_subdomain_owned,
+                                           data = q_l2g)
+                    ar.create_dataset_sync(nodes_dataset,
+                                           offsets=mesh.globalMesh.elementOffsets_subdomain_owned*Xdmf_NodesPerElement,
+                                           data = x[:mesh.nElements_owned].reshape((mesh.nElements_owned*Xdmf_NodesPerElement,3)))
                 else:
-                    SubElement(elements,"xi:include",{"parse":"text","href":"./"+ar.textDataDir+"/elements"+spaceSuffix+str(tCount)+".txt"})
-                    SubElement(nodes,"xi:include",{"parse":"text","href":"./"+ar.textDataDir+"/nodes"+spaceSuffix+str(tCount)+".txt"})
-                    if init or meshChanged:
-                        #this will fail if elements_dgp1 already exists
-                        #q_l2g = numpy.zeros((Xdmf_NumberOfElements,Xdmf_NodesPerElement),'i')
-                        #brute force to start
-                        #for eN in range(Xdmf_NumberOfElements):
-                        #    for nN in range(Xdmf_NodesPerElement):
-                        #        q_l2g[eN,nN] = eN*Xdmf_NodesPerElement + nN
-                        #
-                        q_l2g = numpy.arange(Xdmf_NumberOfElements*Xdmf_NodesPerElement,dtype=GENERATED_INDEX_DTYPE).reshape((Xdmf_NumberOfElements,Xdmf_NodesPerElement))
-                        numpy.savetxt(ar.textDataDir+"/elements"+spaceSuffix+str(tCount)+".txt",q_l2g,fmt='%d')
-                        numpy.savetxt(ar.textDataDir+"/nodes"+spaceSuffix+str(tCount)+".txt",x.flat[:])
-
-                    #
-                #hdfile
-        #need to write a grid
+                    q_l2g = numpy.arange(Xdmf_NumberOfElements*Xdmf_NodesPerElement,dtype=GENERATED_INDEX_DTYPE).reshape((Xdmf_NumberOfElements,Xdmf_NodesPerElement))
+                    if ar.hdfFile is not None:
+                        ar.create_dataset_async(elements_dataset, data = q_l2g)
+                        ar.create_dataset_async(nodes_dataset, data = x.flat[:])
+                    else:
+                        numpy.savetxt(ar.textDataDir+"/"+elements_stem+".txt",q_l2g,fmt='%d')
+                        numpy.savetxt(ar.textDataDir+"/"+nodes_stem+".txt",x.flat[:])
         return self.arGrid
     #def
     def writeScalarXdmf_quadrature(self,ar,u,name,tCount=0,init=True):
@@ -1225,7 +1233,7 @@ def writeMeshXdmf_DGP1Lagrange(self,ar,name,mesh,spaceDim,dofMap,CGDOFMap,t=0.0,
         spaceSuffix = "_dgp1_Lagrange"
         gridName = self.setGridCollectionAndGridElements(init,ar,arGrid,t,spaceSuffix)
 
-        if self.arGrid is None or self.arTime.get('Value') != "{0:e}".format(t):
+        if self._needsGrid(t):
             if spaceDim == 1:
                 Xdmf_ElementTopology = "Polyline"
             elif spaceDim == 2:
@@ -1233,84 +1241,56 @@ def writeMeshXdmf_DGP1Lagrange(self,ar,name,mesh,spaceDim,dofMap,CGDOFMap,t=0.0,
             elif spaceDim == 3:
                 Xdmf_ElementTopology = "Tetrahedron"
             if ar.global_sync:
-                self.arGrid = SubElement(self.arGridCollection,"Grid",{"Name":gridName,"GridType":"Uniform"})
-                self.arTime = SubElement(self.arGrid,"Time",{"Value":"{0:e}".format(t),"Name":"{0:d}".format(tCount)})
-                topology    = SubElement(self.arGrid,"Topology",
-                                         {"Type":Xdmf_ElementTopology,
-                                          "NumberOfElements": "{0:d}".format(mesh.globalMesh.nElements_global)})
-                elements    = SubElement(topology,"DataItem",
-                                         {"Format":ar.dataItemFormat,
-                                          **xdmf_number_type(dofMap.dof_offsets_subdomain_owned,dofMap.l2g),
-                                          "Dimensions":"{0:d} {1:d}".format(mesh.globalMesh.nElements_global,mesh.nNodes_element)})
-                geometry    = SubElement(self.arGrid,"Geometry",{"Type":"XYZ"})
-                nodes       = SubElement(geometry,"DataItem",
-                                         {"Format":ar.dataItemFormat,
-                                          **xdmf_number_type(mesh.nodeArray),
-                                          "Dimensions":"{0:d} {1:d}".format(dofMap.nDOF_all_processes,3)})
-                if ar.hdfFile is not None:
-                    elements.text = ar.hdfFilename+":/elements"+spaceSuffix+str(tCount)
-                    nodes.text    = ar.hdfFilename+":/nodes"+spaceSuffix+str(tCount)
-                    if init or meshChanged:
-                        #this will fail if elements_dgp1 already exists
-                        ar.create_dataset_sync('elements{0}{1:d}'.format(spaceSuffix,tCount),
-                                               offsets = mesh.globalMesh.elementOffsets_subdomain_owned,
-                                               data = dofMap.dof_offsets_subdomain_owned[ar.rank]+dofMap.l2g[:mesh.globalMesh.elementOffsets_subdomain_owned[ar.rank+1]
-                                                                                                             -mesh.globalMesh.elementOffsets_subdomain_owned[ar.rank]])
-                        #bad
-                        dgnodes = numpy.zeros((dofMap.dof_offsets_subdomain_owned[ar.rank+1]-dofMap.dof_offsets_subdomain_owned[ar.rank],3),mesh.nodeArray.dtype)
-                        for eN in range(mesh.globalMesh.elementOffsets_subdomain_owned[ar.rank+1]
-                                        -mesh.globalMesh.elementOffsets_subdomain_owned[ar.rank]):
-                            for nN in range(mesh.nNodes_element):
-                                dgnodes[dofMap.l2g[eN,nN],:]=mesh.nodeArray[mesh.elementNodesArray[eN,nN]]
-                        #make more pythonic loop
-                        ar.create_dataset_sync('nodes{0}{1:d}'.format(spaceSuffix,tCount),
-                                               offsets=dofMap.dof_offsets_subdomain_owned,
-                                               data = dgnodes)
-                else:
-                    assert False, "Global sync  not implemented for text heavy data"
+                assert ar.hdfFile is not None, "Global sync  not implemented for text heavy data"
+                n_elements       = mesh.globalMesh.nElements_global
+                elements_dims    = [n_elements, mesh.nNodes_element]
+                nodes_dims       = [dofMap.nDOF_all_processes, 3]
+                elements_dataset = 'elements'+spaceSuffix+str(tCount)
+                nodes_dataset    = 'nodes'+spaceSuffix+str(tCount)
+                elements_dtype   = (dofMap.dof_offsets_subdomain_owned, dofMap.l2g)
             else:
-                self.arGrid = SubElement(self.arGridCollection,"Grid",{"Name":gridName,"GridType":"Uniform"})
-                self.arTime = SubElement(self.arGrid,"Time",{"Value":"%e" % (t,),"Name":str(tCount)})
-                topology    = SubElement(self.arGrid,"Topology",
-                                         {"Type":Xdmf_ElementTopology,
-                                          "NumberOfElements":str(mesh.nElements_global)})
-                elements    = SubElement(topology,"DataItem",
-                                         {"Format":ar.dataItemFormat,
-                                          **xdmf_number_type(dofMap.l2g),
-                                          "Dimensions":"%i %i" % (mesh.nElements_global,mesh.nNodes_element)})
-                geometry    = SubElement(self.arGrid,"Geometry",{"Type":"XYZ"})
-                nodes       = SubElement(geometry,"DataItem",
-                                         {"Format":ar.dataItemFormat,
-                                          **xdmf_number_type(mesh.nodeArray),
-                                          "Dimensions":"%i %i" % (dofMap.nDOF,3)})
-                if ar.hdfFile is not None:
-                    elements.text = ar.hdfFilename+":/elements"+str(ar.rank)+spaceSuffix+str(tCount)
-                    nodes.text    = ar.hdfFilename+":/nodes"+str(ar.rank)+spaceSuffix+str(tCount)
-                    if init or meshChanged:
-                        #this will fail if elements_dgp1 already exists
-                        ar.create_dataset_async('elements'+str(ar.rank)+spaceSuffix+str(tCount), data = dofMap.l2g)
-                        #bad
-                        dgnodes = numpy.zeros((dofMap.nDOF,3),mesh.nodeArray.dtype)
-                        for eN in range(mesh.nElements_global):
-                            for nN in range(mesh.nNodes_element):
-                                dgnodes[dofMap.l2g[eN,nN],:]=mesh.nodeArray[mesh.elementNodesArray[eN,nN]]
-                        ar.create_dataset_async('nodes'+str(ar.rank)+spaceSuffix+str(tCount), data = dgnodes)
+                n_elements       = mesh.nElements_global
+                elements_dims    = [n_elements, mesh.nNodes_element]
+                nodes_dims       = [dofMap.nDOF, 3]
+                elements_dataset = 'elements'+str(ar.rank)+spaceSuffix+str(tCount)
+                nodes_dataset    = 'nodes'+str(ar.rank)+spaceSuffix+str(tCount)
+                elements_dtype   = dofMap.l2g
+            elements_stem = 'elements'+spaceSuffix+str(tCount)
+            nodes_stem    = 'nodes'+spaceSuffix+str(tCount)
+
+            self.arGrid = ar.write_grid(self.arGridCollection, gridName, t, tCount)
+            ar.write_topology(self.arGrid, Xdmf_ElementTopology, n_elements,
+                              elements_dims, elements_dataset, elements_stem,
+                              dtype=elements_dtype)
+            ar.write_geometry(self.arGrid, nodes_dims, nodes_dataset, nodes_stem,
+                              dtype=mesh.nodeArray)
+            if init or meshChanged:
+                if ar.global_sync:
+                    nElements_owned = (mesh.globalMesh.elementOffsets_subdomain_owned[ar.rank+1]
+                                       -mesh.globalMesh.elementOffsets_subdomain_owned[ar.rank])
+                    #this will fail if elements_dgp1 already exists
+                    ar.create_dataset_sync(elements_dataset,
+                                           offsets = mesh.globalMesh.elementOffsets_subdomain_owned,
+                                           data = dofMap.dof_offsets_subdomain_owned[ar.rank]+dofMap.l2g[:nElements_owned])
+                    dgnodes = numpy.zeros((dofMap.dof_offsets_subdomain_owned[ar.rank+1]-dofMap.dof_offsets_subdomain_owned[ar.rank],3),mesh.nodeArray.dtype)
+                    for eN in range(nElements_owned):
+                        for nN in range(mesh.nNodes_element):
+                            dgnodes[dofMap.l2g[eN,nN],:]=mesh.nodeArray[mesh.elementNodesArray[eN,nN]]
+                    ar.create_dataset_sync(nodes_dataset,
+                                           offsets=dofMap.dof_offsets_subdomain_owned,
+                                           data = dgnodes)
                 else:
-                    SubElement(elements,"xi:include",{"parse":"text","href":"./"+ar.textDataDir+"/elements"+spaceSuffix+str(tCount)+".txt"})
-                    SubElement(nodes,"xi:include",{"parse":"text","href":"./"+ar.textDataDir+"/nodes"+spaceSuffix+str(tCount)+".txt"})
-                    if init or meshChanged:
-                        numpy.savetxt(ar.textDataDir+"/elements"+spaceSuffix+str(tCount)+".txt",dofMap.l2g,fmt='%d')
-                        #bad
-                        dgnodes = numpy.zeros((dofMap.nDOF,3),mesh.nodeArray.dtype)
-                        for eN in range(mesh.nElements_global):
-                            for nN in range(mesh.nNodes_element):
-                                dgnodes[dofMap.l2g[eN,nN],:]=mesh.nodeArray[mesh.elementNodesArray[eN,nN]]
-                        #make more pythonic loop
-                        numpy.savetxt(ar.textDataDir+"/nodes"+spaceSuffix+str(tCount)+".txt",dgnodes)
-
-                    #
-                #hdfile
-            #need to write a grid
+                    dgnodes = numpy.zeros((dofMap.nDOF,3),mesh.nodeArray.dtype)
+                    for eN in range(mesh.nElements_global):
+                        for nN in range(mesh.nNodes_element):
+                            dgnodes[dofMap.l2g[eN,nN],:]=mesh.nodeArray[mesh.elementNodesArray[eN,nN]]
+                    if ar.hdfFile is not None:
+                        #this will fail if elements_dgp1 already exists
+                        ar.create_dataset_async(elements_dataset, data = dofMap.l2g)
+                        ar.create_dataset_async(nodes_dataset, data = dgnodes)
+                    else:
+                        numpy.savetxt(ar.textDataDir+"/"+elements_stem+".txt",dofMap.l2g,fmt='%d')
+                        numpy.savetxt(ar.textDataDir+"/"+nodes_stem+".txt",dgnodes)
         return self.arGrid
     #def
     def writeMeshXdmf_DGP2Lagrange(self,ar,name,mesh,spaceDim,dofMap,CGDOFMap,t=0.0,
@@ -1322,7 +1302,7 @@ def writeMeshXdmf_DGP2Lagrange(self,ar,name,mesh,spaceDim,dofMap,CGDOFMap,t=0.0,
         spaceSuffix = "_dgp2_Lagrange"
         gridName = self.setGridCollectionAndGridElements(init,ar,arGrid,t,spaceSuffix)
 
-        if self.arGrid is None or self.arTime.get('Value') != "{0:e}".format(t):
+        if self._needsGrid(t):
             if spaceDim == 1:
                 Xdmf_ElementTopology = "Edge_3"
             elif spaceDim == 2:
@@ -1330,105 +1310,55 @@ def writeMeshXdmf_DGP2Lagrange(self,ar,name,mesh,spaceDim,dofMap,CGDOFMap,t=0.0,
             elif spaceDim == 3:
                 Xdmf_ElementTopology = "Tet_10"
             if ar.global_sync:
-                self.arGrid = SubElement(self.arGridCollection,"Grid",{"Name":gridName,"GridType":"Uniform"})
-                self.arTime = SubElement(self.arGrid,"Time",{"Value":"{0:e}".format(t),"Name":"{0:d}".format(tCount)})
-                topology    = SubElement(self.arGrid,"Topology",
-                                         {"Type":Xdmf_ElementTopology,
-                                          "NumberOfElements":"{0:d}".format(mesh.globalMesh.nElements_global)})
-                elements    = SubElement(topology,"DataItem",
-                                         {"Format":ar.dataItemFormat,
-                                          **xdmf_number_type(dofMap.dof_offsets_subdomain_owned,dofMap.l2g),
-                                          "Dimensions":"{0:d} {1:d}".format(mesh.globalMesh.nElements_global,dofMap.l2g.shape[-1])})
-                geometry    = SubElement(self.arGrid,"Geometry",{"Type":"XYZ"})
-                #try to use fancy functions later
-                nodes       = SubElement(geometry,"DataItem",
-                                         {"Format":ar.dataItemFormat,
-                                          **xdmf_number_type(CGDOFMap.lagrangeNodesArray),
-                                          "Dimensions":"{0:d} {1:d}".format(dofMap.nDOF_all_processes,3)})
-                if ar.hdfFile is not None:
-                    elements.text = ar.hdfFilename+":/elements"+spaceSuffix+str(tCount)
-                    nodes.text    = ar.hdfFilename+":/nodes"+spaceSuffix+str(tCount)
-                    if init or meshChanged:
-                        #this will fail if elements_dgp1 already exists
-                        ar.create_dataset_sync('elements'+spaceSuffix+str(tCount),
-                                               offsets = mesh.globalMesh.elementOffsets_subdomain_owned,
-                                               data = dofMap.dof_offsets_subdomain_owned[ar.rank]+dofMap.l2g[:mesh.globalMesh.elementOffsets_subdomain_owned[ar.rank+1]-
-                                                                                                             mesh.globalMesh.elementOffsets_subdomain_owned[ar.rank]])
-                        #bad
-                        dgnodes = numpy.zeros((dofMap.dof_offsets_subdomain_owned[ar.rank+1]-dofMap.dof_offsets_subdomain_owned[ar.rank],3),CGDOFMap.lagrangeNodesArray.dtype)
-                        for eN in range(mesh.globalMesh.elementOffsets_subdomain_owned[ar.rank+1]-mesh.globalMesh.elementOffsets_subdomain_owned[ar.rank]):
-                            for i in range(dofMap.l2g.shape[1]):
-                                #for nN in range(mesh.nNodes_element):
-                                #dgnodes[dofMap.l2g[eN,nN],:]=mesh.nodeArray[mesh.elementNodesArray[eN,nN]]
-                                #now changed lagrange nodes to hold all nodes
-                                dgnodes[dofMap.l2g[eN,i],:]= CGDOFMap.lagrangeNodesArray[CGDOFMap.l2g[eN,i],:]
-                            #next loop over extra dofs on element and write out
-                            #for nN in range(mesh.nNodes_element,dofMap.l2g.shape[1]):
-                            #    dgnodes[dofMap.l2g[eN,nN],:]= CGDOFMap.lagrangeNodesArray[CGDOFMap.l2g[eN,nN]-mesh.nNodes_global,:]
-
-                        #make more pythonic loop
-                        ar.create_dataset_sync('nodes'+spaceSuffix+str(tCount),
-                                               offsets = dofMap.dof_offsets_subdomain_owned,
-                                               data = dgnodes)
-                else:
-                    assert False, "global_sync not implemented for  text heavy data"
+                assert ar.hdfFile is not None, "global_sync not implemented for  text heavy data"
+                n_elements       = mesh.globalMesh.nElements_global
+                elements_dims    = [n_elements, dofMap.l2g.shape[-1]]
+                nodes_dims       = [dofMap.nDOF_all_processes, 3]
+                elements_dataset = 'elements'+spaceSuffix+str(tCount)
+                nodes_dataset    = 'nodes'+spaceSuffix+str(tCount)
+                elements_dtype   = (dofMap.dof_offsets_subdomain_owned, dofMap.l2g)
             else:
-                self.arGrid = SubElement(self.arGridCollection,"Grid",{"Name":gridName,"GridType":"Uniform"})
-                self.arTime = SubElement(self.arGrid,"Time",{"Value":"%e" % (t,),"Name":str(tCount)})
-                topology    = SubElement(self.arGrid,"Topology",
-                                         {"Type":Xdmf_ElementTopology,
-                                          "NumberOfElements":str(mesh.nElements_global)})
-                elements    = SubElement(topology,"DataItem",
-                                         {"Format":ar.dataItemFormat,
-                                          **xdmf_number_type(dofMap.l2g),
-                                          "Dimensions":"%i %i" % dofMap.l2g.shape})
-                geometry    = SubElement(self.arGrid,"Geometry",{"Type":"XYZ"})
-                #try to use fancy functions later
-                nodes       = SubElement(geometry,"DataItem",
-                                         {"Format":ar.dataItemFormat,
-                                          **xdmf_number_type(CGDOFMap.lagrangeNodesArray),
-                                          "Dimensions":"%i %i" % (dofMap.nDOF,3)})
-                if ar.hdfFile is not None:
-                    elements.text = ar.hdfFilename+":/elements"+str(ar.rank)+spaceSuffix+str(tCount)
-                    nodes.text    = ar.hdfFilename+":/nodes"+str(ar.rank)+spaceSuffix+str(tCount)
-                    if init or meshChanged:
-                        #this will fail if elements_dgp1 already exists
-                        ar.create_dataset_async('elements'+str(ar.rank)+spaceSuffix+str(tCount), data = dofMap.l2g)
-                        #bad
-                        dgnodes = numpy.zeros((dofMap.nDOF,3),CGDOFMap.lagrangeNodesArray.dtype)
-                        for eN in range(mesh.nElements_global):
-                            for i in range(dofMap.l2g.shape[1]):
-                                #for nN in range(mesh.nNodes_element):
-                                #dgnodes[dofMap.l2g[eN,nN],:]=mesh.nodeArray[mesh.elementNodesArray[eN,nN]]
-                                #now changed lagrange nodes to hold all nodes
-                                dgnodes[dofMap.l2g[eN,i],:]= CGDOFMap.lagrangeNodesArray[CGDOFMap.l2g[eN,i],:]
-                            #next loop over extra dofs on element and write out
-                            #for nN in range(mesh.nNodes_element,dofMap.l2g.shape[1]):
-                            #    dgnodes[dofMap.l2g[eN,nN],:]= CGDOFMap.lagrangeNodesArray[CGDOFMap.l2g[eN,nN]-mesh.nNodes_global,:]
-
-                        ar.create_dataset_async('nodes'+str(ar.rank)+spaceSuffix+str(tCount), data = dgnodes)
+                n_elements       = mesh.nElements_global
+                elements_dims    = list(dofMap.l2g.shape)
+                nodes_dims       = [dofMap.nDOF, 3]
+                elements_dataset = 'elements'+str(ar.rank)+spaceSuffix+str(tCount)
+                nodes_dataset    = 'nodes'+str(ar.rank)+spaceSuffix+str(tCount)
+                elements_dtype   = dofMap.l2g
+            elements_stem = 'elements'+spaceSuffix+str(tCount)
+            nodes_stem    = 'nodes'+spaceSuffix+str(tCount)
+
+            self.arGrid = ar.write_grid(self.arGridCollection, gridName, t, tCount)
+            ar.write_topology(self.arGrid, Xdmf_ElementTopology, n_elements,
+                              elements_dims, elements_dataset, elements_stem,
+                              dtype=elements_dtype)
+            ar.write_geometry(self.arGrid, nodes_dims, nodes_dataset, nodes_stem,
+                              dtype=CGDOFMap.lagrangeNodesArray)
+            if init or meshChanged:
+                if ar.global_sync:
+                    nElements_owned = (mesh.globalMesh.elementOffsets_subdomain_owned[ar.rank+1]
+                                       -mesh.globalMesh.elementOffsets_subdomain_owned[ar.rank])
+                    ar.create_dataset_sync(elements_dataset,
+                                           offsets = mesh.globalMesh.elementOffsets_subdomain_owned,
+                                           data = dofMap.dof_offsets_subdomain_owned[ar.rank]+dofMap.l2g[:nElements_owned])
+                    dgnodes = numpy.zeros((dofMap.dof_offsets_subdomain_owned[ar.rank+1]-dofMap.dof_offsets_subdomain_owned[ar.rank],3),CGDOFMap.lagrangeNodesArray.dtype)
+                    for eN in range(nElements_owned):
+                        for i in range(dofMap.l2g.shape[1]):
+                            #lagrange nodes hold all nodes
+                            dgnodes[dofMap.l2g[eN,i],:]= CGDOFMap.lagrangeNodesArray[CGDOFMap.l2g[eN,i],:]
+                    ar.create_dataset_sync(nodes_dataset,
+                                           offsets = dofMap.dof_offsets_subdomain_owned,
+                                           data = dgnodes)
                 else:
-                    SubElement(elements,"xi:include",{"parse":"text","href":"./"+ar.textDataDir+"/elements"+spaceSuffix+str(tCount)+".txt"})
-                    SubElement(nodes,"xi:include",{"parse":"text","href":"./"+ar.textDataDir+"/nodes"+spaceSuffix+str(tCount)+".txt"})
-                    if init or meshChanged:
-                        numpy.savetxt(ar.textDataDir+"/elements"+spaceSuffix+str(tCount)+".txt",dofMap.l2g,fmt='%d')
-                        #bad
-                        dgnodes = numpy.zeros((dofMap.nDOF,3),CGDOFMap.lagrangeNodesArray.dtype)
-                        for eN in range(mesh.nElements_global):
-                            for i in range(dofMap.l2g.shape[1]):
-                                #for nN in range(mesh.nNodes_element):
-                                #dgnodes[dofMap.l2g[eN,nN],:]=mesh.nodeArray[mesh.elementNodesArray[eN,nN]]
-                                dgnodes[dofMap.l2g[eN,i],:]= CGDOFMap.lagrangeNodesArray[CGDOFMap.l2g[eN,i],:]
-                            #next loop over extra dofs on element and write out
-                            #for nN in range(mesh.nNodes_element,dofMap.l2g.shape[1]):
-                            #    dgnodes[dofMap.l2g[eN,nN],:]= CGDOFMap.lagrangeNodesArray[CGDOFMap.l2g[eN,nN]-mesh.nNodes_global,:]
-
-                        #make more pythonic loop
-                        numpy.savetxt(ar.textDataDir+"/nodes"+spaceSuffix+str(tCount)+".txt",dgnodes)
-
-                    #
-                #hdfile
-            #need to write a grid
+                    dgnodes = numpy.zeros((dofMap.nDOF,3),CGDOFMap.lagrangeNodesArray.dtype)
+                    for eN in range(mesh.nElements_global):
+                        for i in range(dofMap.l2g.shape[1]):
+                            dgnodes[dofMap.l2g[eN,i],:]= CGDOFMap.lagrangeNodesArray[CGDOFMap.l2g[eN,i],:]
+                    if ar.hdfFile is not None:
+                        ar.create_dataset_async(elements_dataset, data = dofMap.l2g)
+                        ar.create_dataset_async(nodes_dataset, data = dgnodes)
+                    else:
+                        numpy.savetxt(ar.textDataDir+"/"+elements_stem+".txt",dofMap.l2g,fmt='%d')
+                        numpy.savetxt(ar.textDataDir+"/"+nodes_stem+".txt",dgnodes)
         return self.arGrid
 
     def writeMeshXdmf_C0P2Lagrange(self,ar,name,mesh,spaceDim,dofMap,t=0.0,
@@ -1440,7 +1370,7 @@ def writeMeshXdmf_C0P2Lagrange(self,ar,name,mesh,spaceDim,dofMap,t=0.0,
         mesh.writeMeshXdmf(ar,"Spatial_Domain",t,init,meshChanged,tCount=tCount)
         spaceSuffix = "_c0p2_Lagrange"
         gridName = self.setGridCollectionAndGridElements(init,ar,arGrid,t,spaceSuffix)
-        if self.arGrid is None or self.arTime.get('Value') != "{0:e}".format(t):
+        if self._needsGrid(t):
             if spaceDim == 1:
                 Xdmf_ElementTopology = "Edge_3"
             elif spaceDim == 2:
@@ -1468,8 +1398,7 @@ def writeMeshXdmf_C0P2Lagrange(self,ar,name,mesh,spaceDim,dofMap,t=0.0,
             elements_stem = 'elements'+spaceSuffix+str(tCount)
             nodes_stem    = 'nodes'+spaceSuffix+str(tCount)
 
-            self.arGrid, self.arTime = ar.write_grid(self.arGridCollection,
-                                                     gridName, t, tCount)
+            self.arGrid = ar.write_grid(self.arGridCollection, gridName, t, tCount)
             ar.write_topology(self.arGrid, Xdmf_ElementTopology, n_elements,
                               elements_dims, elements_dataset, elements_stem,
                               dtype=(dofMap.subdomain2global if ar.global_sync
@@ -1526,7 +1455,7 @@ def writeMeshXdmf_C0Q2Lagrange(self,ar,name,mesh,spaceDim,dofMap,t=0.0,init=Fals
         #mesh.writeMeshXdmf(ar,"Spatial_Domain",t,init,meshChanged,tCount=tCount)
         spaceSuffix = "_c0q2_Lagrange"
         gridName = self.setGridCollectionAndGridElements(init,ar,arGrid,t,spaceSuffix)
-        if self.arGrid is None or self.arTime.get('Value') != "{0:e}".format(t):
+        if self._needsGrid(t):
             if spaceDim == 1:
                 print("No writeMeshXdmf_C0Q2Lagrange for 1D")
                 return 0
@@ -1554,71 +1483,110 @@ def writeMeshXdmf_C0Q2Lagrange(self,ar,name,mesh,spaceDim,dofMap,t=0.0,init=Fals
                     for i in range(8): #loop over subelements
                         for j in range(8): # loop over nodes of subelement
                             l2g[8*eN+i,j] = dofs[e2s[i][j]]
+            lagrangeNodesArray = dofMap.lagrangeNodesArray
             if ar.global_sync:
-                self.arGrid = SubElement(self.arGridCollection,"Grid",{"Name":gridName,"GridType":"Uniform"})
-                self.arTime = SubElement(self.arGrid,"Time",{"Value":"%e" % (t,),"Name":str(tCount)})
-
-                lagrangeNodesArray = dofMap.lagrangeNodesArray
-
-                topology = SubElement(self.arGrid,"Topology",
-                                      {"Type":Xdmf_ElementTopology,
-                                       "NumberOfElements":str(mesh.globalMesh.nElements_global*nsubelements)})
-                elements = SubElement(topology,"DataItem",
-                                      {"Format":ar.dataItemFormat,
-                                       **xdmf_number_type(dofMap.subdomain2global),
-                                       "Dimensions":"%i %i" % (mesh.globalMesh.nElements_global*nsubelements, l2g.shape[-1])})
-                geometry = SubElement(self.arGrid,"Geometry",{"Type":"XYZ"})
-
-                allNodes = SubElement(geometry,"DataItem",
-                                      {"Format":ar.dataItemFormat,
-                                       **xdmf_number_type(lagrangeNodesArray),
-                                       "Dimensions":"%i %i" % (dofMap.nDOF_all_processes, lagrangeNodesArray.shape[-1])})
-                if ar.hdfFile is not None:
-                    elements.text = ar.hdfFilename+":/elements"+spaceSuffix+str(tCount)
-                    allNodes.text = ar.hdfFilename+":/nodes"+spaceSuffix+str(tCount)
-                    import copy
-                    if init or meshChanged:
-                        ar.create_dataset_sync('elements'+spaceSuffix+str(tCount),
-                                               offsets = [x*nsubelements for x in mesh.globalMesh.elementOffsets_subdomain_owned],
-                                               data = dofMap.subdomain2global[l2g[:mesh.nElements_owned*nsubelements]])
-                        ar.create_dataset_sync('nodes'+spaceSuffix+str(tCount),
-                                               offsets = dofMap.dof_offsets_subdomain_owned,
-                                               data = lagrangeNodesArray[:dofMap.dof_offsets_subdomain_owned[ar.rank+1]-dofMap.dof_offsets_subdomain_owned[ar.rank]])
-                else:
-                    assert False, "not implemented  for text heavy data"
+                assert ar.hdfFile is not None, "not implemented  for text heavy data"
+                n_elements       = mesh.globalMesh.nElements_global*nsubelements
+                elements_dims    = [n_elements, l2g.shape[-1]]
+                nodes_dims       = [dofMap.nDOF_all_processes, lagrangeNodesArray.shape[-1]]
+                elements_dataset = 'elements'+spaceSuffix+str(tCount)
+                nodes_dataset    = 'nodes'+spaceSuffix+str(tCount)
+                elements_dtype   = dofMap.subdomain2global
             else:
-                self.arGrid = SubElement(self.arGridCollection,"Grid",{"Name":gridName,"GridType":"Uniform"})
-                self.arTime = SubElement(self.arGrid,"Time",{"Value":"%e" % (t,),"Name":str(tCount)})
-
-                lagrangeNodesArray = dofMap.lagrangeNodesArray
-
-                topology = SubElement(self.arGrid,"Topology",
-                                      {"Type":Xdmf_ElementTopology,
-                                       "NumberOfElements":str(l2g.shape[0])})
-                elements = SubElement(topology,"DataItem",
-                                      {"Format":ar.dataItemFormat,
-                                       **xdmf_number_type(l2g),
-                                       "Dimensions":"%i %i" % l2g.shape})
-                geometry = SubElement(self.arGrid,"Geometry",{"Type":"XYZ"})
-
-                allNodes = SubElement(geometry,"DataItem",
-                                      {"Format":ar.dataItemFormat,
-                                       **xdmf_number_type(lagrangeNodesArray),
-                                       "Dimensions":"%i %i" % lagrangeNodesArray.shape})
+                n_elements       = l2g.shape[0]
+                elements_dims    = list(l2g.shape)
+                nodes_dims       = list(lagrangeNodesArray.shape)
+                elements_dataset = 'elements'+str(ar.rank)+spaceSuffix+str(tCount)
+                nodes_dataset    = 'nodes'+str(ar.rank)+spaceSuffix+str(tCount)
+                elements_dtype   = l2g
+            elements_stem = 'elements'+spaceSuffix+str(tCount)
+            nodes_stem    = 'nodes'+spaceSuffix+str(tCount)
 
-                if ar.hdfFile is not None:
-                    elements.text = ar.hdfFilename+":/elements"+str(ar.rank)+spaceSuffix+str(tCount)
-                    allNodes.text = ar.hdfFilename+":/nodes"+str(ar.rank)+spaceSuffix+str(tCount)
-                    import copy
-                    if init or meshChanged:
-                        ar.create_dataset_async('elements'+str(ar.rank)+spaceSuffix+str(tCount), data = l2g)
-                        ar.create_dataset_async('nodes'+str(ar.rank)+spaceSuffix+str(tCount), data = lagrangeNodesArray)
+            self.arGrid = ar.write_grid(self.arGridCollection, gridName, t, tCount)
+            ar.write_topology(self.arGrid, Xdmf_ElementTopology, n_elements,
+                              elements_dims, elements_dataset, elements_stem,
+                              dtype=elements_dtype)
+            ar.write_geometry(self.arGrid, nodes_dims, nodes_dataset, nodes_stem,
+                              dtype=lagrangeNodesArray)
+            if init or meshChanged:
+                if ar.global_sync:
+                    ar.create_dataset_sync(elements_dataset,
+                                           offsets = [x*nsubelements for x in mesh.globalMesh.elementOffsets_subdomain_owned],
+                                           data = dofMap.subdomain2global[l2g[:mesh.nElements_owned*nsubelements]])
+                    ar.create_dataset_sync(nodes_dataset,
+                                           offsets = dofMap.dof_offsets_subdomain_owned,
+                                           data = lagrangeNodesArray[:dofMap.dof_offsets_subdomain_owned[ar.rank+1]-dofMap.dof_offsets_subdomain_owned[ar.rank]])
+                elif ar.hdfFile is not None:
+                    ar.create_dataset_async(elements_dataset, data = l2g)
+                    ar.create_dataset_async(nodes_dataset, data = lagrangeNodesArray)
                 else:
-                    SubElement(elements,"xi:include",{"parse":"text","href":"./"+ar.textDataDir+"/elements"+spaceSuffix+str(tCount)+".txt"})
-                    SubElement(allNodes,"xi:include",{"parse":"text","href":"./"+ar.textDataDir+"/nodes"+spaceSuffix+str(tCount)+".txt"})
-                    if init or meshChanged:
-                        numpy.savetxt(ar.textDataDir+"/elements"+spaceSuffix+str(tCount)+".txt",l2g,fmt='%d')
-                        numpy.savetxt(ar.textDataDir+"/nodes"+spaceSuffix+str(tCount)+".txt",lagrangeNodesArray)
+                    numpy.savetxt(ar.textDataDir+"/"+elements_stem+".txt",l2g,fmt='%d')
+                    numpy.savetxt(ar.textDataDir+"/"+nodes_stem+".txt",lagrangeNodesArray)
+        return self.arGrid
+
+    def _writeDiscontinuousVertexMesh(self,ar,gridName,mesh,spaceDim,t,tCount,
+                                      spaceSuffix,init,meshChanged,
+                                      nodes_every_pass=False):
+        """The element-by-element vertex mesh several spaces display on.
+
+        Each element gets its own copy of its vertices, numbered
+        consecutively, so a field with one value per element vertex (a
+        non-conforming P1 function written as discontinuous P1, an RT0
+        velocity) can be shown. Shared by the Crouzeix-Raviart and
+        lowest-order mixed writers, which differed only in dataset naming
+        and in whether the nodes were rewritten on every pass.
+        """
+        if spaceDim == 1:
+            Xdmf_ElementTopology = "Polyline"
+        elif spaceDim == 2:
+            Xdmf_ElementTopology = "Triangle"
+        elif spaceDim == 3:
+            Xdmf_ElementTopology = "Tetrahedron"
+        Xdmf_NodesPerElement = spaceDim+1
+        if ar.global_sync:
+            assert ar.hdfFile is not None, "global_sync not implemented for text heavy data"
+            Xdmf_NumberOfElements= mesh.globalMesh.nElements_global
+            elements_dataset = 'elements'+spaceSuffix+str(tCount)
+            nodes_dataset    = 'nodes'+spaceSuffix+str(tCount)
+        else:
+            Xdmf_NumberOfElements= mesh.nElements_global
+            elements_dataset = 'elements'+str(ar.rank)+spaceSuffix+str(tCount)
+            nodes_dataset    = 'nodes'+str(ar.rank)+spaceSuffix+str(tCount)
+        elements_stem = 'elements'+spaceSuffix+str(tCount)
+        nodes_stem    = 'nodes'+spaceSuffix+str(tCount)
+
+        self.arGrid = ar.write_grid(self.arGridCollection, gridName, t, tCount)
+        ar.write_topology(self.arGrid, Xdmf_ElementTopology, Xdmf_NumberOfElements,
+                          [Xdmf_NumberOfElements, Xdmf_NodesPerElement],
+                          elements_dataset, elements_stem,
+                          dtype=GENERATED_INDEX_DTYPE)
+        ar.write_geometry(self.arGrid, [Xdmf_NumberOfElements*Xdmf_NodesPerElement, 3],
+                          nodes_dataset, nodes_stem, dtype=mesh.nodeArray)
+        if ar.global_sync:
+            if init or meshChanged:
+                #simple dg l2g mapping
+                dg_l2g = numpy.arange(mesh.nElements_owned*Xdmf_NodesPerElement,dtype=GENERATED_INDEX_DTYPE).reshape((mesh.nElements_owned,Xdmf_NodesPerElement))
+                ar.create_dataset_sync(elements_dataset,
+                                       offsets = mesh.globalMesh.elementOffsets_subdomain_owned,
+                                       data = dg_l2g)
+            if init or meshChanged or nodes_every_pass:
+                dgnodes = numpy.reshape(mesh.nodeArray[mesh.elementNodesArray[:mesh.nElements_owned]],
+                                        (mesh.nElements_owned*Xdmf_NodesPerElement,3))
+                ar.create_dataset_sync(nodes_dataset,
+                                       offsets = mesh.globalMesh.elementOffsets_subdomain_owned*Xdmf_NodesPerElement,
+                                       data = dgnodes)
+        elif ar.hdfFile is not None:
+            if init or meshChanged:
+                dg_l2g = numpy.arange(Xdmf_NumberOfElements*Xdmf_NodesPerElement,dtype=GENERATED_INDEX_DTYPE).reshape((Xdmf_NumberOfElements,Xdmf_NodesPerElement))
+                ar.create_dataset_async(elements_dataset, data = dg_l2g)
+            if init or meshChanged or nodes_every_pass:
+                dgnodes = numpy.reshape(mesh.nodeArray[mesh.elementNodesArray],(Xdmf_NumberOfElements*Xdmf_NodesPerElement,3))
+                ar.create_dataset_async(nodes_dataset, data = dgnodes)
+        elif init or meshChanged:
+            dg_l2g = numpy.arange(Xdmf_NumberOfElements*Xdmf_NodesPerElement,dtype=GENERATED_INDEX_DTYPE).reshape((Xdmf_NumberOfElements,Xdmf_NodesPerElement))
+            numpy.savetxt(ar.textDataDir+"/"+elements_stem+".txt",dg_l2g,fmt='%d')
+            dgnodes = numpy.reshape(mesh.nodeArray[mesh.elementNodesArray],(Xdmf_NumberOfElements*Xdmf_NodesPerElement,3))
+            numpy.savetxt(ar.textDataDir+"/"+nodes_stem+".txt",dgnodes)
         return self.arGrid
 
     def writeMeshXdmf_CrouzeixRaviartP1(self,ar,mesh,spaceDim,dofMap,t=0.0,
@@ -1635,86 +1603,9 @@ def writeMeshXdmf_CrouzeixRaviartP1(self,ar,mesh,spaceDim,dofMap,t=0.0,
         spaceSuffix = "_ncp1_CrouzeixRaviart"
         gridName = self.setGridCollectionAndGridElements(init,ar,arGrid,t,spaceSuffix)
 
-        if self.arGrid is None or self.arTime.get('Value') != "{0:e}".format(t):
-            if spaceDim == 1:
-                Xdmf_ElementTopology = "Polyline"
-            elif spaceDim == 2:
-                Xdmf_ElementTopology = "Triangle"
-            elif spaceDim == 3:
-                Xdmf_ElementTopology = "Tetrahedron"
-            if ar.global_sync:
-                Xdmf_NodesPerElement = spaceDim+1
-                Xdmf_NumberOfElements= mesh.globalMesh.nElements_global
-                self.arGrid = SubElement(self.arGridCollection,"Grid",{"Name":gridName,"GridType":"Uniform"})
-                self.arTime = SubElement(self.arGrid,"Time",{"Value":"%e" % (t,),"Name":str(tCount)})
-                topology    = SubElement(self.arGrid,"Topology",
-                                         {"Type":Xdmf_ElementTopology,
-                                          "NumberOfElements":str(Xdmf_NumberOfElements)})
-                elements    = SubElement(topology,"DataItem",
-                                         {"Format":ar.dataItemFormat,
-                                          **xdmf_number_type(GENERATED_INDEX_DTYPE),
-                                          "Dimensions":"%i %i" % (Xdmf_NumberOfElements,Xdmf_NodesPerElement)})
-                geometry    = SubElement(self.arGrid,"Geometry",{"Type":"XYZ"})
-                nodes       = SubElement(geometry,"DataItem",
-                                         {"Format":ar.dataItemFormat,
-                                          **xdmf_number_type(mesh.nodeArray),
-                                          "Dimensions":"%i %i" % (Xdmf_NumberOfElements*Xdmf_NodesPerElement,3)})
-                if ar.hdfFile is not None:
-                    elements.text = ar.hdfFilename+":/elements"+spaceSuffix+str(tCount)
-                    nodes.text    = ar.hdfFilename+":/nodes"+spaceSuffix+str(tCount)
-                    if init or meshChanged:
-                        #simple dg l2g mapping
-                        dg_l2g = numpy.arange(mesh.nElements_owned*Xdmf_NodesPerElement,dtype=GENERATED_INDEX_DTYPE).reshape((mesh.nElements_owned,Xdmf_NodesPerElement))
-                        ar.create_dataset_sync('elements'+spaceSuffix+str(tCount),
-                                               offsets = mesh.globalMesh.elementOffsets_subdomain_owned,
-                                               data = dg_l2g)
-                        dgnodes = numpy.reshape(mesh.nodeArray[mesh.elementNodesArray[:mesh.nElements_owned]],
-                                                (mesh.nElements_owned*Xdmf_NodesPerElement,3))
-                        ar.create_dataset_sync('nodes'+spaceSuffix+str(tCount),
-                                               offsets = mesh.globalMesh.elementOffsets_subdomain_owned*Xdmf_NodesPerElement,
-                                               data = dgnodes)
-                else:
-                    assert False, "global_sync not implemented for text heavy data"
-            else:
-                Xdmf_NodesPerElement = spaceDim+1
-                Xdmf_NumberOfElements= mesh.nElements_global
-                self.arGrid = SubElement(self.arGridCollection,"Grid",{"Name":gridName,"GridType":"Uniform"})
-                self.arTime = SubElement(self.arGrid,"Time",{"Value":"%e" % (t,),"Name":str(tCount)})
-                topology    = SubElement(self.arGrid,"Topology",
-                                         {"Type":Xdmf_ElementTopology,
-                                          "NumberOfElements":str(Xdmf_NumberOfElements)})
-                elements    = SubElement(topology,"DataItem",
-                                         {"Format":ar.dataItemFormat,
-                                          **xdmf_number_type(GENERATED_INDEX_DTYPE),
-                                          "Dimensions":"%i %i" % (Xdmf_NumberOfElements,Xdmf_NodesPerElement)})
-                geometry    = SubElement(self.arGrid,"Geometry",{"Type":"XYZ"})
-                nodes       = SubElement(geometry,"DataItem",
-                                         {"Format":ar.dataItemFormat,
-                                          **xdmf_number_type(mesh.nodeArray),
-                                          "Dimensions":"%i %i" % (Xdmf_NumberOfElements*Xdmf_NodesPerElement,3)})
-                if ar.hdfFile is not None:
-                    elements.text = ar.hdfFilename+":/elements"+str(ar.rank)+spaceSuffix+str(tCount)
-                    nodes.text    = ar.hdfFilename+":/nodes"+str(ar.rank)+spaceSuffix+str(tCount)
-                    if init or meshChanged:
-                        #simple dg l2g mapping
-                        dg_l2g = numpy.arange(Xdmf_NumberOfElements*Xdmf_NodesPerElement,dtype=GENERATED_INDEX_DTYPE).reshape((Xdmf_NumberOfElements,Xdmf_NodesPerElement))
-                        ar.create_dataset_async('elements'+str(ar.rank)+spaceSuffix+str(tCount), data = dg_l2g)
-
-                        dgnodes = numpy.reshape(mesh.nodeArray[mesh.elementNodesArray],(Xdmf_NumberOfElements*Xdmf_NodesPerElement,3))
-                        ar.create_dataset_async('nodes'+str(ar.rank)+spaceSuffix+str(tCount), data = dgnodes)
-                else:
-                    SubElement(elements,"xi:include",{"parse":"text","href":"./"+ar.textDataDir+"/elements"+spaceSuffix+str(tCount)+".txt"})
-                    SubElement(nodes,"xi:include",{"parse":"text","href":"./"+ar.textDataDir+"/nodes"+spaceSuffix+str(tCount)+".txt"})
-                    if init or meshChanged:
-                        dg_l2g = numpy.arange(Xdmf_NumberOfElements*Xdmf_NodesPerElement,dtype=GENERATED_INDEX_DTYPE).reshape((Xdmf_NumberOfElements,Xdmf_NodesPerElement))
-                        numpy.savetxt(ar.textDataDir+"/elements"+spaceSuffix+str(tCount)+".txt",dg_l2g,fmt='%d')
-
-                        dgnodes = numpy.reshape(mesh.nodeArray[mesh.elementNodesArray],(Xdmf_NumberOfElements*Xdmf_NodesPerElement,3))
-                        numpy.savetxt(ar.textDataDir+"/nodes"+spaceSuffix+str(tCount)+".txt",dgnodes)
-
-                    #
-                #hdfile
-            #need to write a grid
+        if self._needsGrid(t):
+            self._writeDiscontinuousVertexMesh(ar,gridName,mesh,spaceDim,t,tCount,
+                                               spaceSuffix,init,meshChanged)
         return self.arGrid
     #def
     def writeFunctionXdmf_DGP1Lagrange(self,ar,u,tCount=0,init=True, dofMap=None):
@@ -1752,185 +1643,98 @@ def writeFunctionXdmf_DGP2Lagrange(self,ar,u,tCount=0,init=True, dofMap=None):
                            dataset=dataset)
 
     def writeFunctionXdmf_CrouzeixRaviartP1(self,ar,u,tCount=0,init=True, dofMap=None):
-        if ar.global_sync:
-            Xdmf_NumberOfElements = u.femSpace.mesh.globalMesh.nElements_global
-            Xdmf_NodesPerElement  = u.femSpace.mesh.nNodes_element
-            name = u.name.replace(' ','_')
-            #if writing as dgp1
-            attribute = SubElement(self.arGrid,"Attribute",{"Name":name,
-                                                            "AttributeType":"Scalar",
-                                                            "Center":"Node"})
-            values    = SubElement(attribute,"DataItem",
-                                   {"Format":ar.dataItemFormat,
-                                    **xdmf_number_type(u.dof),
-                                    "Dimensions":"%i" % (Xdmf_NumberOfElements*Xdmf_NodesPerElement,)})
-            nElements_owned = u.femSpace.mesh.globalMesh.elementOffsets_subdomain_owned[ar.rank+1] - u.femSpace.mesh.globalMesh.elementOffsets_subdomain_owned[ar.rank]
-            u_tmp = numpy.zeros((nElements_owned*Xdmf_NodesPerElement,),u.dof.dtype)
-            if u.femSpace.nSpace_global == 1:
-                for eN in range(nElements_owned):
-                    #dof associated with face id, so opposite usual C0P1 ordering here
-                    u_tmp[eN*Xdmf_NodesPerElement + 0] = u.dof[u.femSpace.dofMap.l2g[eN,1]]
-                    u_tmp[eN*Xdmf_NodesPerElement + 1] = u.dof[u.femSpace.dofMap.l2g[eN,0]]
-            elif u.femSpace.nSpace_global == 2:
-                for eN in range(nElements_owned):
-                    #assume vertex associated with face across from it
-                    u_tmp[eN*Xdmf_NodesPerElement + 0] = u.dof[u.femSpace.dofMap.l2g[eN,1]]
-                    u_tmp[eN*Xdmf_NodesPerElement + 0]+= u.dof[u.femSpace.dofMap.l2g[eN,2]]
-                    u_tmp[eN*Xdmf_NodesPerElement + 0]-= u.dof[u.femSpace.dofMap.l2g[eN,0]]
-
-                    u_tmp[eN*Xdmf_NodesPerElement + 1] = u.dof[u.femSpace.dofMap.l2g[eN,0]]
-                    u_tmp[eN*Xdmf_NodesPerElement + 1]+= u.dof[u.femSpace.dofMap.l2g[eN,2]]
-                    u_tmp[eN*Xdmf_NodesPerElement + 1]-= u.dof[u.femSpace.dofMap.l2g[eN,1]]
-
-                    u_tmp[eN*Xdmf_NodesPerElement + 2] = u.dof[u.femSpace.dofMap.l2g[eN,0]]
-                    u_tmp[eN*Xdmf_NodesPerElement + 2]+= u.dof[u.femSpace.dofMap.l2g[eN,1]]
-                    u_tmp[eN*Xdmf_NodesPerElement + 2]-= u.dof[u.femSpace.dofMap.l2g[eN,2]]
-            else:
-                for eN in range(nElements_owned):
-                    for i in range(Xdmf_NodesPerElement):
-                        u_tmp[eN*Xdmf_NodesPerElement + i] = u.dof[u.femSpace.dofMap.l2g[eN,i]]*(1.0-float(u.femSpace.nSpace_global)) + \
-                                                             sum([u.dof[u.femSpace.dofMap.l2g[eN,j]] for j in range(Xdmf_NodesPerElement) if j != i])
-
+        """A non-conforming P1 function, two ways.
 
-            if ar.hdfFile is not None:
-                values.text = ar.hdfFilename+":/"+name+"_t"+str(tCount)
-                ar.create_dataset_sync(name+"_t"+str(tCount),
-                                       offsets = u.femSpace.mesh.globalMesh.elementOffsets_subdomain_owned*Xdmf_NodesPerElement,
-                                       data = u_tmp)
-            else:
-                assert False, "global_sync not implemented for text heavy data"
-            #if writing as face centered (true nc p1)
-            #could be under self.arGrid or mesh.arGrid
-            grid = u.femSpace.mesh.arGrid #self.arGrid
-            attribute_dof = SubElement(grid,"Attribute",{"Name":name+"_dof",
-                                                         "AttributeType":"Scalar",
-                                                         "Center":"Face"})
-            values_dof    = SubElement(attribute_dof,"DataItem",
-                                   {"Format":ar.dataItemFormat,
-                                    **xdmf_number_type(u.dof),
-                                    "Dimensions":"%i" % (dofMap.nDOF_all_processes,)})
-            if ar.hdfFile is not None:
-                values_dof.text = ar.hdfFilename+":/"+name+"_dof"+"_t"+str(tCount)
-                ar.create_dataset_sync(name+"_dof"+"_t"+str(tCount),
-                                       offsets = dofMap.dof_offsets_subdomain_owned,
-                                       data = u.dof[:dofMap.dof_offsets_subdomain_owned[ar.rank+1] - dofMap.dof_offsets_subdomain_owned[ar.rank]])
-            else:
-                assert False, "global_sync not implemented for text heavy data"
+        As a discontinuous P1 function on the element-by-element vertex
+        mesh (the vertex value is the sum over the element's faces minus
+        the opposite face's), which viewers can draw; and as the dofs
+        themselves, face-centred on the base mesh.
+        """
+        mesh = u.femSpace.mesh
+        Xdmf_NodesPerElement  = mesh.nNodes_element
+        name = u.name.replace(' ','_')
+        if ar.global_sync:
+            Xdmf_NumberOfElements = mesh.globalMesh.nElements_global
+            nElements_owned = (mesh.globalMesh.elementOffsets_subdomain_owned[ar.rank+1]
+                               - mesh.globalMesh.elementOffsets_subdomain_owned[ar.rank])
         else:
-            Xdmf_NumberOfElements = u.femSpace.mesh.nElements_global
-            Xdmf_NodesPerElement  = u.femSpace.mesh.nNodes_element
-
-            name = u.name.replace(' ','_')
-            #if writing as dgp1
-            attribute = SubElement(self.arGrid,"Attribute",{"Name":name,
-                                                            "AttributeType":"Scalar",
-                                                            "Center":"Node"})
-            values    = SubElement(attribute,"DataItem",
-                                   {"Format":ar.dataItemFormat,
-                                    **xdmf_number_type(u.dof),
-                                    "Dimensions":"%i" % (Xdmf_NumberOfElements*Xdmf_NodesPerElement,)})
-            u_tmp = numpy.zeros((Xdmf_NumberOfElements*Xdmf_NodesPerElement,),u.dof.dtype)
-            if u.femSpace.nSpace_global == 1:
-                for eN in range(Xdmf_NumberOfElements):
-                    #dof associated with face id, so opposite usual C0P1 ordering here
-                    u_tmp[eN*Xdmf_NodesPerElement + 0] = u.dof[u.femSpace.dofMap.l2g[eN,1]]
-                    u_tmp[eN*Xdmf_NodesPerElement + 1] = u.dof[u.femSpace.dofMap.l2g[eN,0]]
-            elif u.femSpace.nSpace_global == 2:
-                for eN in range(Xdmf_NumberOfElements):
-                    #assume vertex associated with face across from it
-                    u_tmp[eN*Xdmf_NodesPerElement + 0] = u.dof[u.femSpace.dofMap.l2g[eN,1]]
-                    u_tmp[eN*Xdmf_NodesPerElement + 0]+= u.dof[u.femSpace.dofMap.l2g[eN,2]]
-                    u_tmp[eN*Xdmf_NodesPerElement + 0]-= u.dof[u.femSpace.dofMap.l2g[eN,0]]
-
-                    u_tmp[eN*Xdmf_NodesPerElement + 1] = u.dof[u.femSpace.dofMap.l2g[eN,0]]
-                    u_tmp[eN*Xdmf_NodesPerElement + 1]+= u.dof[u.femSpace.dofMap.l2g[eN,2]]
-                    u_tmp[eN*Xdmf_NodesPerElement + 1]-= u.dof[u.femSpace.dofMap.l2g[eN,1]]
-
-                    u_tmp[eN*Xdmf_NodesPerElement + 2] = u.dof[u.femSpace.dofMap.l2g[eN,0]]
-                    u_tmp[eN*Xdmf_NodesPerElement + 2]+= u.dof[u.femSpace.dofMap.l2g[eN,1]]
-                    u_tmp[eN*Xdmf_NodesPerElement + 2]-= u.dof[u.femSpace.dofMap.l2g[eN,2]]
-            else:
-                for eN in range(Xdmf_NumberOfElements):
-                    for i in range(Xdmf_NodesPerElement):
-                        u_tmp[eN*Xdmf_NodesPerElement + i] = u.dof[u.femSpace.dofMap.l2g[eN,i]]*(1.0-float(u.femSpace.nSpace_global)) + \
-                                                             sum([u.dof[u.femSpace.dofMap.l2g[eN,j]] for j in range(Xdmf_NodesPerElement) if j != i])
-
+            Xdmf_NumberOfElements = mesh.nElements_global
+            nElements_owned = Xdmf_NumberOfElements
+        u_tmp = numpy.zeros((nElements_owned*Xdmf_NodesPerElement,),u.dof.dtype)
+        l2g = u.femSpace.dofMap.l2g
+        if u.femSpace.nSpace_global == 1:
+            for eN in range(nElements_owned):
+                #dof associated with face id, so opposite usual C0P1 ordering here
+                u_tmp[eN*Xdmf_NodesPerElement + 0] = u.dof[l2g[eN,1]]
+                u_tmp[eN*Xdmf_NodesPerElement + 1] = u.dof[l2g[eN,0]]
+        elif u.femSpace.nSpace_global == 2:
+            for eN in range(nElements_owned):
+                #assume vertex associated with face across from it
+                u_tmp[eN*Xdmf_NodesPerElement + 0] = u.dof[l2g[eN,1]]
+                u_tmp[eN*Xdmf_NodesPerElement + 0]+= u.dof[l2g[eN,2]]
+                u_tmp[eN*Xdmf_NodesPerElement + 0]-= u.dof[l2g[eN,0]]
+
+                u_tmp[eN*Xdmf_NodesPerElement + 1] = u.dof[l2g[eN,0]]
+                u_tmp[eN*Xdmf_NodesPerElement + 1]+= u.dof[l2g[eN,2]]
+                u_tmp[eN*Xdmf_NodesPerElement + 1]-= u.dof[l2g[eN,1]]
+
+                u_tmp[eN*Xdmf_NodesPerElement + 2] = u.dof[l2g[eN,0]]
+                u_tmp[eN*Xdmf_NodesPerElement + 2]+= u.dof[l2g[eN,1]]
+                u_tmp[eN*Xdmf_NodesPerElement + 2]-= u.dof[l2g[eN,2]]
+        else:
+            for eN in range(nElements_owned):
+                for i in range(Xdmf_NodesPerElement):
+                    u_tmp[eN*Xdmf_NodesPerElement + i] = u.dof[l2g[eN,i]]*(1.0-float(u.femSpace.nSpace_global)) + \
+                                                         sum([u.dof[l2g[eN,j]] for j in range(Xdmf_NodesPerElement) if j != i])
+        #as dgp1, on this writer's grid
+        if ar.global_sync:
+            ar.write_field(self.arGrid, name, u_tmp, tCount,
+                           dimensions=[Xdmf_NumberOfElements*Xdmf_NodesPerElement],
+                           sync_offsets=mesh.globalMesh.elementOffsets_subdomain_owned*Xdmf_NodesPerElement,
+                           sync_data=u_tmp)
+        else:
+            ar.write_field(self.arGrid, name, u_tmp, tCount)
+        #as face centered dofs (true nc p1), on the mesh's grid
+        grid = mesh.arGrid
+        if ar.global_sync:
+            owned = dofMap.dof_offsets_subdomain_owned
+            ar.write_field(grid, name+"_dof", u.dof, tCount, center="Face",
+                           dimensions=[dofMap.nDOF_all_processes],
+                           sync_offsets=owned,
+                           sync_data=u.dof[:owned[ar.rank+1] - owned[ar.rank]])
+        else:
+            ar.write_field(grid, name+"_dof", u.dof, tCount, center="Face",
+                           dimensions=[u.nDOF_global])
 
-            if ar.hdfFile is not None:
-                values.text = ar.hdfFilename+":/"+name+"_p"+str(ar.rank)+"_t"+str(tCount)
-                ar.create_dataset_async(name+"_p"+str(ar.rank)+"_t"+str(tCount), data = u_tmp)
-            else:
-                numpy.savetxt(ar.textDataDir+"/"+name+str(tCount)+".txt",u_tmp)
-                SubElement(values,"xi:include",{"parse":"text","href":"./"+ar.textDataDir+"/"+name+str(tCount)+".txt"})
-
-
-            #if writing as face centered (true nc p1)
-            #could be under self.arGrid or mesh.arGrid
-            grid = u.femSpace.mesh.arGrid #self.arGrid
-            attribute_dof = SubElement(grid,"Attribute",{"Name":name+"_dof",
-                                                         "AttributeType":"Scalar",
-                                                         "Center":"Face"})
-            values_dof    = SubElement(attribute_dof,"DataItem",
-                                   {"Format":ar.dataItemFormat,
-                                    **xdmf_number_type(u.dof),
-                                    "Dimensions":"%i" % (u.nDOF_global,)})
-            if ar.hdfFile is not None:
-                values_dof.text = ar.hdfFilename+":/"+name+"_dof"+"_p"+str(ar.rank)+"_t"+str(tCount)
-                ar.create_dataset_async(name+"_dof"+"_p"+str(ar.rank)+"_t"+str(tCount), data = u.dof)
-            else:
-                numpy.savetxt(ar.textDataDir+"/"+name+"_dof"+str(tCount)+".txt",u.dof)
-                SubElement(values_dof,"xi:include",{"parse":"text","href":"./"+ar.textDataDir+"/"+name+"_dof"+str(tCount)+".txt"})
+    def _nodalVector(self,uList,components,n_rows=None):
+        """A 3-component array from up to three dof vectors; missing components are zero."""
+        u_dof = uList[components[0]].dof
+        rows = slice(None) if n_rows is None else slice(0, n_rows)
+        u_dof = u_dof[rows]
+        if len(components) < 2:
+            v_dof = numpy.zeros(u_dof.shape,dtype='d')
+        else:
+            v_dof = uList[components[1]].dof[rows]
+        if len(components) < 3:
+            w_dof = numpy.zeros(u_dof.shape,dtype='d')
+        else:
+            w_dof = uList[components[2]].dof[rows]
+        return numpy.column_stack((u_dof,v_dof,w_dof))
 
     def writeVectorFunctionXdmf_nodal(self,ar,uList,components,vectorName,spaceSuffix,tCount=0,init=True):
         nDOF_global = uList[components[0]].nDOF_global
+        velocity = self._nodalVector(uList,components)
         if ar.global_sync:
-            attribute = SubElement(self.arGrid,"Attribute",{"Name":vectorName,
-                                                            "AttributeType":"Vector",
-                                                            "Center":"Node"})
-            u_dof = uList[components[0]].dof
-            if len(components) < 2:
-                v_dof = numpy.zeros(u_dof.shape,dtype='d')
-            else:
-                v_dof = uList[components[1]].dof
-            if len(components) < 3:
-                w_dof = numpy.zeros(u_dof.shape,dtype='d')
-            else:
-                w_dof = uList[components[2]].dof
-            velocity = numpy.column_stack((u_dof,v_dof,w_dof))
-            #after velocity, so the DataItem describes what is written; it
-            #declared 4-byte floats over these 8-byte ones
-            values    = SubElement(attribute,"DataItem",
-                                   {"Format":ar.dataItemFormat,
-                                    **xdmf_number_type(velocity),
-                                    "Dimensions":"%i %i" % (uList[components[0]].femSpace.dofMap.nDOF_all_processes,3)})
-            if ar.hdfFile is not None:
-                values.text = ar.hdfFilename+":/"+vectorName+"_t"+str(tCount)
-                nDOF_owned = (uList[components[0]].femSpace.dofMap.dof_offsets_subdomain_owned[ar.rank+1] -
-                              uList[components[0]].femSpace.dofMap.dof_offsets_subdomain_owned[ar.rank] )
-                ar.create_dataset_sync(vectorName+"_t"+str(tCount),
-                                       offsets = uList[components[0]].femSpace.dofMap.dof_offsets_subdomain_owned,
-                                       data = velocity[:nDOF_owned])
+            dofMap = uList[components[0]].femSpace.dofMap
+            nDOF_owned = (dofMap.dof_offsets_subdomain_owned[ar.rank+1] -
+                          dofMap.dof_offsets_subdomain_owned[ar.rank] )
+            ar.write_field(self.arGrid, vectorName, velocity, tCount, rank="Vector",
+                           dimensions=[dofMap.nDOF_all_processes,3],
+                           sync_offsets=dofMap.dof_offsets_subdomain_owned,
+                           sync_data=velocity[:nDOF_owned])
         else:
-            attribute = SubElement(self.arGrid,"Attribute",{"Name":vectorName,
-                                                            "AttributeType":"Vector",
-                                                            "Center":"Node"})
-            u_dof = uList[components[0]].dof
-            if len(components) < 2:
-                v_dof = numpy.zeros(u_dof.shape,dtype='d')
-            else:
-                v_dof = uList[components[1]].dof
-            if len(components) < 3:
-                w_dof = numpy.zeros(u_dof.shape,dtype='d')
-            else:
-                w_dof = uList[components[2]].dof
-            velocity = numpy.column_stack((u_dof,v_dof,w_dof))
-            values    = SubElement(attribute,"DataItem",
-                                   {"Format":ar.dataItemFormat,
-                                    **xdmf_number_type(velocity),
-                                    "Dimensions":"%i %i" % (nDOF_global,3)})
-            if ar.hdfFile is not None:
-                values.text = ar.hdfFilename+":/"+vectorName+"_p"+str(ar.rank)+"_t"+str(tCount)
-                ar.create_dataset_async(vectorName+"_p"+str(ar.rank)+"_t"+str(tCount), data = velocity)
+            ar.write_field(self.arGrid, vectorName, velocity, tCount, rank="Vector",
+                           dimensions=[nDOF_global,3])
 
     def writeVectorFunctionXdmf_CrouzeixRaviartP1(self,ar,uList,components,spaceSuffix,tCount=0,init=True):
         return self.writeVectorFunctionXdmf_nodal(ar,uList,components,"_ncp1_CrouzeixRaviart",
@@ -1951,79 +1755,47 @@ def writeMeshXdmf_MonomialDGPK(self,ar,mesh,spaceDim,interpolationPoints,t=0.0,
         gridName = self.setGridCollectionAndGridElements(init,ar,arGrid,t,spaceSuffix)
 
         assert len(interpolationPoints.shape) == 3 #make sure have right type of dictionary
-        if self.arGrid is None or self.arTime.get('Value') != "{0:e}".format(t):
-            Xdmf_ElementTopology = "Polyvertex"
+        self.mesh = mesh
+        if self._needsGrid(t):
+            Xdmf_NodesPerElement = interpolationPoints.shape[1]
             if ar.global_sync:
+                assert ar.hdfFile is not None, "global_sync not implementedf or text heavy data"
                 Xdmf_NumberOfElements= mesh.globalMesh.nElements_global
-                Xdmf_NodesPerElement = interpolationPoints.shape[1]
-                Xdmf_NodesGlobal     = Xdmf_NumberOfElements*Xdmf_NodesPerElement
-
-                self.arGrid = SubElement(self.arGridCollection,"Grid",{"Name":gridName,"GridType":"Uniform"})
-                self.arTime = SubElement(self.arGrid,"Time",{"Value":"%e" % (t,),"Name":str(tCount)})
-                topology    = SubElement(self.arGrid,"Topology",
-                                         {"Type":Xdmf_ElementTopology,
-                                          "NumberOfElements":str(Xdmf_NumberOfElements),
-                                          "NodesPerElement":str(Xdmf_NodesPerElement)})
-                elements    = SubElement(topology,"DataItem",
-                                         {"Format":ar.dataItemFormat,
-                                          **xdmf_number_type(GENERATED_INDEX_DTYPE),
-                                          "Dimensions":"%i %i" % (Xdmf_NumberOfElements,Xdmf_NodesPerElement)})
-                geometry    = SubElement(self.arGrid,"Geometry",{"Type":"XYZ"})
-                nodes       = SubElement(geometry,"DataItem",
-                                         {"Format":ar.dataItemFormat,
-                                          **xdmf_number_type(interpolationPoints),
-                                          "Dimensions":"%i %i" % (Xdmf_NodesGlobal,3)})
-                if ar.hdfFile is not None:
-                    elements.text = ar.hdfFilename+":/elements"+spaceSuffix+str(tCount)
-                    nodes.text    = ar.hdfFilename+":/nodes"+spaceSuffix+str(tCount)
-                    if init or meshChanged:
-                        q_l2g = numpy.arange(mesh.nElements_owned*Xdmf_NodesPerElement,dtype=GENERATED_INDEX_DTYPE).reshape((mesh.nElements_owned,Xdmf_NodesPerElement))
-                        ar.create_dataset_sync('elements'+spaceSuffix+str(tCount),
-                                               offsets = mesh.globalMesh.elementOffsets_subdomain_owned,
-                                               data = q_l2g)
-                        ar.create_dataset_sync('nodes'+spaceSuffix+str(tCount),
-                                               offsets = mesh.globalMesh.elementOffsets_subdomain_owned,
-                                               data = interpolationPoints[:mesh.nElements_owned])
-                else:
-                    assert False, "global_sync not implementedf or text heavy data"
+                elements_dataset = 'elements'+spaceSuffix+str(tCount)
+                nodes_dataset    = 'nodes'+spaceSuffix+str(tCount)
             else:
                 Xdmf_NumberOfElements= interpolationPoints.shape[0]
-                Xdmf_NodesPerElement = interpolationPoints.shape[1]
-                Xdmf_NodesGlobal     = Xdmf_NumberOfElements*Xdmf_NodesPerElement
-
-                self.arGrid = SubElement(self.arGridCollection,"Grid",{"Name":gridName,"GridType":"Uniform"})
-                self.arTime = SubElement(self.arGrid,"Time",{"Value":"%e" % (t,),"Name":str(tCount)})
-                topology    = SubElement(self.arGrid,"Topology",
-                                         {"Type":Xdmf_ElementTopology,
-                                          "NumberOfElements":str(Xdmf_NumberOfElements),
-                                          "NodesPerElement":str(Xdmf_NodesPerElement)})
-                elements    = SubElement(topology,"DataItem",
-                                         {"Format":ar.dataItemFormat,
-                                          **xdmf_number_type(GENERATED_INDEX_DTYPE),
-                                          "Dimensions":"%i %i" % (Xdmf_NumberOfElements,Xdmf_NodesPerElement)})
-                geometry    = SubElement(self.arGrid,"Geometry",{"Type":"XYZ"})
-                nodes       = SubElement(geometry,"DataItem",
-                                         {"Format":ar.dataItemFormat,
-                                          **xdmf_number_type(interpolationPoints),
-                                          "Dimensions":"%i %i" % (Xdmf_NodesGlobal,3)})
-                if ar.hdfFile is not None:
-                    elements.text = ar.hdfFilename+":/elements"+str(ar.rank)+spaceSuffix+str(tCount)
-                    nodes.text    = ar.hdfFilename+":/nodes"+str(ar.rank)+spaceSuffix+str(tCount)
-                    if init or meshChanged:
-                        q_l2g = numpy.arange(Xdmf_NumberOfElements*Xdmf_NodesPerElement,dtype=GENERATED_INDEX_DTYPE).reshape((Xdmf_NumberOfElements,Xdmf_NodesPerElement))
-                        ar.create_dataset_async('elements'+str(ar.rank)+spaceSuffix+str(tCount), data = q_l2g)
-                        ar.create_dataset_async('nodes'+str(ar.rank)+spaceSuffix+str(tCount), data = interpolationPoints.flat[:])
+                elements_dataset = 'elements'+str(ar.rank)+spaceSuffix+str(tCount)
+                nodes_dataset    = 'nodes'+str(ar.rank)+spaceSuffix+str(tCount)
+            Xdmf_NodesGlobal     = Xdmf_NumberOfElements*Xdmf_NodesPerElement
+            elements_stem = 'elements'+spaceSuffix+str(tCount)
+            nodes_stem    = 'nodes'+spaceSuffix+str(tCount)
+
+            self.arGrid = ar.write_grid(self.arGridCollection, gridName, t, tCount)
+            ar.write_topology(self.arGrid, "Polyvertex", Xdmf_NumberOfElements,
+                              [Xdmf_NumberOfElements, Xdmf_NodesPerElement],
+                              elements_dataset, elements_stem,
+                              nodes_per_element=Xdmf_NodesPerElement,
+                              dtype=GENERATED_INDEX_DTYPE)
+            ar.write_geometry(self.arGrid, [Xdmf_NodesGlobal, 3],
+                              nodes_dataset, nodes_stem, dtype=interpolationPoints)
+            if init or meshChanged:
+                if ar.global_sync:
+                    q_l2g = numpy.arange(mesh.nElements_owned*Xdmf_NodesPerElement,dtype=GENERATED_INDEX_DTYPE).reshape((mesh.nElements_owned,Xdmf_NodesPerElement))
+                    ar.create_dataset_sync(elements_dataset,
+                                           offsets = mesh.globalMesh.elementOffsets_subdomain_owned,
+                                           data = q_l2g)
+                    ar.create_dataset_sync(nodes_dataset,
+                                           offsets = mesh.globalMesh.elementOffsets_subdomain_owned,
+                                           data = interpolationPoints[:mesh.nElements_owned])
                 else:
-                    SubElement(elements,"xi:include",{"parse":"text","href":"./"+ar.textDataDir+"/elements"+spaceSuffix+str(tCount)+".txt"})
-                    SubElement(nodes,"xi:include",{"parse":"text","href":"./"+ar.textDataDir+"/nodes"+spaceSuffix+str(tCount)+".txt"})
-                    if init or meshChanged:
-                        q_l2g = numpy.arange(Xdmf_NumberOfElements*Xdmf_NodesPerElement,dtype=GENERATED_INDEX_DTYPE).reshape((Xdmf_NumberOfElements,Xdmf_NodesPerElement))
-                        numpy.savetxt(ar.textDataDir+"/elements"+spaceSuffix+str(tCount)+".txt",q_l2g,fmt='%d')
-                        numpy.savetxt(ar.textDataDir+"/nodes"+spaceSuffix+str(tCount)+".txt",interpolationPoints.flat[:])
-
-                    #
-                #hdfile
-            #need to write a grid
+                    q_l2g = numpy.arange(Xdmf_NumberOfElements*Xdmf_NodesPerElement,dtype=GENERATED_INDEX_DTYPE).reshape((Xdmf_NumberOfElements,Xdmf_NodesPerElement))
+                    if ar.hdfFile is not None:
+                        ar.create_dataset_async(elements_dataset, data = q_l2g)
+                        ar.create_dataset_async(nodes_dataset, data = interpolationPoints.flat[:])
+                    else:
+                        numpy.savetxt(ar.textDataDir+"/"+elements_stem+".txt",q_l2g,fmt='%d')
+                        numpy.savetxt(ar.textDataDir+"/"+nodes_stem+".txt",interpolationPoints.flat[:])
         return self.arGrid
     #def
     def writeFunctionXdmf_MonomialDGPK(self,ar,interpolationValues,name,tCount=0,init=True, mesh=None):
@@ -2049,48 +1821,27 @@ def writeVectorFunctionXdmf_MonomialDGPK(self,ar,interpolationValues,name,tCount
         need to check on way to save dofs as well
         """
         assert len(interpolationValues.shape) == 3
+        Xdmf_NumberOfComponents = interpolationValues.shape[2]
+        n_points = interpolationValues.shape[1]
         if ar.global_sync:
-            Xdmf_NodesGlobal = self.mesh.globalMesh.nElements_global*interpolationValues.shape[1]
-            Xdmf_NumberOfComponents = interpolationValues.shape[2]
-            attribute = SubElement(self.arGrid,"Attribute",{"Name":name,
-                                                            "AttributeType":"Vector",
-                                                            "Center":"Node"})
-            values    = SubElement(attribute,"DataItem",
-                                   {"Format":ar.dataItemFormat,
-                                    **xdmf_number_type(interpolationValues),
-                                    "Dimensions":"%i %i" % (Xdmf_NodesGlobal,3)})#force 3d vector since points 3d
-            #mwf brute force
-            tmp = numpy.zeros((Xdmf_NodesGlobal,3),interpolationValues.dtype)
-            tmp[:,:Xdmf_NumberOfComponents]=numpy.reshape(interpolationValues.flat,(interpolationValues.shape[0]*interpolationValues.shape[1],
-                                                                                    Xdmf_NumberOfComponents))
-
-            if ar.hdfFile is not None:
-                values.text = ar.hdfFilename+":/"+name+"_t"+str(tCount)
-                ar.create_dataset_sync(name+"_t"+str(tCount),
-                                       offsets = self.mesh.globalMesh.elementOffsets_subdomain_owned*interpolationValues.shape[1],
-                                       data = tmp)
-            else:
-                assert False, "global_sync  not  implemented for text  heavy data"
+            #this rank contributes the elements it owns; the DataItem covers
+            #every element's points globally (the previous version sized the
+            #temporary by the global count and filled it from the local
+            #array, which only agrees on one rank)
+            n_owned = self.mesh.nElements_owned
+            tmp = numpy.zeros((n_owned*n_points,3),interpolationValues.dtype)
+            tmp[:,:Xdmf_NumberOfComponents]=numpy.reshape(interpolationValues[:n_owned].flat,
+                                                          (n_owned*n_points,Xdmf_NumberOfComponents))
+            ar.write_field(self.arGrid, name, tmp, tCount, rank="Vector",
+                           dimensions=[self.mesh.globalMesh.nElements_global*n_points,3],
+                           sync_offsets=self.mesh.globalMesh.elementOffsets_subdomain_owned*n_points,
+                           sync_data=tmp)
         else:
-            Xdmf_NodesGlobal = interpolationValues.shape[0]*interpolationValues.shape[1]
-            Xdmf_NumberOfComponents = interpolationValues.shape[2]
-            attribute = SubElement(self.arGrid,"Attribute",{"Name":name,
-                                                            "AttributeType":"Vector",
-                                                            "Center":"Node"})
-            values    = SubElement(attribute,"DataItem",
-                                   {"Format":ar.dataItemFormat,
-                                    **xdmf_number_type(interpolationValues),
-                                    "Dimensions":"%i %i" % (Xdmf_NodesGlobal,3)})#force 3d vector since points 3d
-            #mwf brute force
+            Xdmf_NodesGlobal = interpolationValues.shape[0]*n_points
             tmp = numpy.zeros((Xdmf_NodesGlobal,3),interpolationValues.dtype)
             tmp[:,:Xdmf_NumberOfComponents]=numpy.reshape(interpolationValues.flat,(Xdmf_NodesGlobal,Xdmf_NumberOfComponents))
-
-            if ar.hdfFile is not None:
-                values.text = ar.hdfFilename+":/"+name+"_p"+str(ar.rank)+"_t"+str(tCount)
-                ar.create_dataset_async(name+"_p"+str(ar.rank)+"_t"+str(tCount), data = tmp)
-            else:
-                numpy.savetxt(ar.textDataDir+"/"+name+str(tCount)+".txt",tmp)
-                SubElement(values,"xi:include",{"parse":"text","href":"./"+ar.textDataDir+"/"+name+str(tCount)+".txt"})
+            ar.write_field(self.arGrid, name, tmp, tCount, rank="Vector",
+                           dimensions=[Xdmf_NodesGlobal,3])
 
 
     def writeMeshXdmf_DGP0(self,ar,mesh,spaceDim,
@@ -2106,93 +1857,41 @@ def writeMeshXdmf_DGP0(self,ar,mesh,spaceDim,
         spaceSuffix = "_dgp0"
         gridName = self.setGridCollectionAndGridElements(init,ar,arGrid,t,spaceSuffix)
 
-        if ar.global_sync:
-            if self.arGrid is None or self.arTime.get('Value') != "{0:e}".format(t):
-                #mwf hack
-                #allow for other types of topologies if the mesh has specified one
-                if 'elementTopologyName' in dir(mesh):
-                    Xdmf_ElementTopology = mesh.elementTopologyName
-                else:
-                    if spaceDim == 1:
-                        Xdmf_ElementTopology = "Polyline"
-                    elif spaceDim == 2:
-                        Xdmf_ElementTopology = "Triangle"
-                    elif spaceDim == 3:
-                        Xdmf_ElementTopology = "Tetrahedron"
-                self.arGrid = SubElement(self.arGridCollection,"Grid",{"Name":gridName,"GridType":"Uniform"})
-                self.arTime = SubElement(self.arGrid,"Time",{"Value":"{0:e}".format(t),"Name":"{0:d}".format(tCount)})
-                topology    = SubElement(self.arGrid,"Topology",
-                                         {"Type":Xdmf_ElementTopology,
-                                          "NumberOfElements":"{0:d}".format(mesh.nElements_global)})
-                #mwf hack, allow for a mixed element mesh
-                if mesh.nNodes_element is None:
-                    assert 'xdmf_topology' in dir(mesh)
-                    elements = SubElement(topology,"DataItem",
-                                          {"Format":ar.dataItemFormat,
-                                           **xdmf_number_type(mesh.globalMesh.nodeNumbering_subdomain2global),
-                                           "Dimensions":"{0:d}".format(len(self.xdmf_topology))})
-                else:
-                    elements    = SubElement(topology,"DataItem",
-                                             {"Format":ar.dataItemFormat,
-                                              **xdmf_number_type(mesh.globalMesh.nodeNumbering_subdomain2global),
-                                              "Dimensions":"{0:d} {1:d}".format(mesh.nElements_global,mesh.nNodes_element)})
-                geometry    = SubElement(self.arGrid,"Geometry",{"Type":"XYZ"})
-                nodes       = SubElement(geometry,"DataItem",
-                                         {"Format":ar.dataItemFormat,
-                                          **xdmf_number_type(mesh.nodeArray),
-                                          "Dimensions":"{0:d} {1:d}".format(mesh.nNodes_global,3)})
-                #just reuse spatial mesh entries
-                if ar.hdfFile is not None:
-                    elements.text = ar.hdfFilename+":/elements"+meshSpaceSuffix+str(tCount)
-                    nodes.text    = ar.hdfFilename+":/nodes"+meshSpaceSuffix+str(tCount)
-                else:
-                    assert False, "global_sync not implemented  for text heavy data"
-                #hdfile
-            #need to write a grid
-        else:
-            if self.arGrid is None or self.arTime.get('Value') != "{0:e}".format(t):
-                #mwf hack
-                #allow for other types of topologies if the mesh has specified one
-                if 'elementTopologyName' in dir(mesh):
-                    Xdmf_ElementTopology = mesh.elementTopologyName
-                else:
-                    if spaceDim == 1:
-                        Xdmf_ElementTopology = "Polyline"
-                    elif spaceDim == 2:
-                        Xdmf_ElementTopology = "Triangle"
-                    elif spaceDim == 3:
-                        Xdmf_ElementTopology = "Tetrahedron"
-                self.arGrid = SubElement(self.arGridCollection,"Grid",{"Name":gridName,"GridType":"Uniform"})
-                self.arTime = SubElement(self.arGrid,"Time",{"Value":"%e" % (t,),"Name":str(tCount)})
-                topology    = SubElement(self.arGrid,"Topology",
-                                         {"Type":Xdmf_ElementTopology,
-                                          "NumberOfElements":str(mesh.nElements_global)})
-                #mwf hack, allow for a mixed element mesh
-                if mesh.nNodes_element is None:
-                    assert 'xdmf_topology' in dir(mesh)
-                    elements = SubElement(topology,"DataItem",
-                                          {"Format":ar.dataItemFormat,
-                                           **xdmf_number_type(mesh.elementNodesArray),
-                                           "Dimensions":"%i" % len(self.xdmf_topology)})
-                else:
-                    elements    = SubElement(topology,"DataItem",
-                                             {"Format":ar.dataItemFormat,
-                                              **xdmf_number_type(mesh.elementNodesArray),
-                                              "Dimensions":"%i %i" % (mesh.nElements_global,mesh.nNodes_element)})
-                geometry    = SubElement(self.arGrid,"Geometry",{"Type":"XYZ"})
-                nodes       = SubElement(geometry,"DataItem",
-                                         {"Format":ar.dataItemFormat,
-                                          **xdmf_number_type(mesh.nodeArray),
-                                          "Dimensions":"%i %i" % (mesh.nNodes_global,3)})
-                #just reuse spatial mesh entries
-                if ar.hdfFile is not None:
-                    elements.text = ar.hdfFilename+":/elements"+str(ar.rank)+meshSpaceSuffix+str(tCount)
-                    nodes.text    = ar.hdfFilename+":/nodes"+str(ar.rank)+meshSpaceSuffix+str(tCount)
-                else:
-                    SubElement(elements,"xi:include",{"parse":"text","href":"./"+ar.textDataDir+"/elements"+meshSpaceSuffix+".txt"})
-                    SubElement(nodes,"xi:include",{"parse":"text","href":"./"+ar.textDataDir+"/nodes"+meshSpaceSuffix+".txt"})
-                #hdfile
-            #need to write a grid
+        if self._needsGrid(t):
+            #mwf hack
+            #allow for other types of topologies if the mesh has specified one
+            if 'elementTopologyName' in dir(mesh):
+                Xdmf_ElementTopology = mesh.elementTopologyName
+            else:
+                if spaceDim == 1:
+                    Xdmf_ElementTopology = "Polyline"
+                elif spaceDim == 2:
+                    Xdmf_ElementTopology = "Triangle"
+                elif spaceDim == 3:
+                    Xdmf_ElementTopology = "Tetrahedron"
+            #the spatial mesh's own datasets are reused; nothing is written here
+            if ar.global_sync:
+                assert ar.hdfFile is not None, "global_sync not implemented  for text heavy data"
+                elements_dataset = 'elements'+meshSpaceSuffix+str(tCount)
+                nodes_dataset    = 'nodes'+meshSpaceSuffix+str(tCount)
+                elements_dtype   = mesh.globalMesh.nodeNumbering_subdomain2global
+            else:
+                elements_dataset = 'elements'+str(ar.rank)+meshSpaceSuffix+str(tCount)
+                nodes_dataset    = 'nodes'+str(ar.rank)+meshSpaceSuffix+str(tCount)
+                elements_dtype   = mesh.elementNodesArray
+            #mwf hack, allow for a mixed element mesh
+            if mesh.nNodes_element is None:
+                assert 'xdmf_topology' in dir(mesh)
+                elements_dims = [len(mesh.xdmf_topology)]
+            else:
+                elements_dims = [mesh.nElements_global, mesh.nNodes_element]
+            self.arGrid = ar.write_grid(self.arGridCollection, gridName, t, tCount)
+            ar.write_topology(self.arGrid, Xdmf_ElementTopology, mesh.nElements_global,
+                              elements_dims, elements_dataset,
+                              'elements'+meshSpaceSuffix, dtype=elements_dtype)
+            ar.write_geometry(self.arGrid, [mesh.nNodes_global, 3],
+                              nodes_dataset, 'nodes'+meshSpaceSuffix,
+                              dtype=mesh.nodeArray)
         return self.arGrid
     #def
     def writeFunctionXdmf_DGP0(self,ar,u,tCount=0,init=True):
@@ -2213,16 +1912,7 @@ def writeVectorFunctionXdmf_DGP0(self,ar,uList,components,vectorName,tCount=0,in
         #default global_sync configuration. The mesh comes from the first
         #component, matching u_dof below.
         mesh = uList[components[0]].femSpace.elementMaps.mesh
-        u_dof = uList[components[0]].dof
-        if len(components) < 2:
-            v_dof = numpy.zeros(u_dof.shape,dtype='d')
-        else:
-            v_dof = uList[components[1]].dof
-        if len(components) < 3:
-            w_dof = numpy.zeros(u_dof.shape,dtype='d')
-        else:
-            w_dof = uList[components[2]].dof
-        velocity = numpy.column_stack((u_dof,v_dof,w_dof))
+        velocity = self._nodalVector(uList,components)
 
         if ar.global_sync:
             ar.write_field(self.arGrid, vectorName, velocity, tCount,
@@ -2247,7 +1937,7 @@ def writeMeshXdmf_P1Bubble(self,ar,mesh,spaceDim,dofMap,t=0.0,
         spaceSuffix = "_c0p1_Bubble%s" % tCount
         gridName = self.setGridCollectionAndGridElements(init,ar,arGrid,t,spaceSuffix)
 
-        if self.arGrid is None or self.arTime.get('Value') != "{0:e}".format(t):
+        if self._needsGrid(t):
             if spaceDim == 1:
                 Xdmf_ElementTopology = "Polyline"
             elif spaceDim == 2:
@@ -2257,113 +1947,55 @@ def writeMeshXdmf_P1Bubble(self,ar,mesh,spaceDim,dofMap,t=0.0,
             Xdmf_NodesPerElement = spaceDim+1 #just handle vertices for now
             Xdmf_NumberOfNodes   = mesh.nNodes_global
             if  ar.global_sync:
+                assert ar.hdfFile is not None, "global_sync not implemented for text heavy data"
                 Xdmf_NumberOfElements= mesh.globalMesh.nElements_global
-                self.arGrid = SubElement(self.arGridCollection,"Grid",{"Name":gridName,"GridType":"Uniform"})
-                self.arTime = SubElement(self.arGrid,"Time",{"Value":"%e" % (t,),"Name":str(tCount)})
-                topology    = SubElement(self.arGrid,"Topology",
-                                         {"Type":Xdmf_ElementTopology,
-                                          "NumberOfElements":str(Xdmf_NumberOfElements)})
-                elements    = SubElement(topology,"DataItem",
-                                         {"Format":ar.dataItemFormat,
-                                          **xdmf_number_type(mesh.elementNodesArray),
-                                          "Dimensions":"%i %i" % (Xdmf_NumberOfElements,Xdmf_NodesPerElement)})
-                geometry    = SubElement(self.arGrid,"Geometry",{"Type":"XYZ"})
-                nodes       = SubElement(geometry,"DataItem",
-                                         {"Format":ar.dataItemFormat,
-                                          **xdmf_number_type(mesh.nodeArray),
-                                          "Dimensions":"%i %i" % (Xdmf_NumberOfNodes,3)})
-                if ar.hdfFile is not None:
-                    elements.text = ar.hdfFilename+":/elements"+spaceSuffix+str(tCount)
-                    nodes.text    = ar.hdfFilename+":/nodes"+spaceSuffix+str(tCount)
-                    if init or meshChanged:
-                        #c0p1 mapping for now
-                        ar.create_dataset_sync('elements'+spaceSuffix+str(tCount),
-                                               offsets = mesh.globalMesh.elementOffsets_subdomain_owned,
-                                               data = mesh.elementNodesArray)
-                        ar.create_dataset_sync('nodes'+spaceSuffix+str(tCount),
-                                               offsets = mesh.globalMesh.nodeOffsets_subdomain_owned,
-                                               data = mesh.nodeArray)
-                else:
-                    assert False, "global_sync not implemented for text heavy data"
+                elements_dataset = 'elements'+spaceSuffix+str(tCount)
+                nodes_dataset    = 'nodes'+spaceSuffix+str(tCount)
             else:
                 Xdmf_NumberOfElements= mesh.nElements_global
-                self.arGrid = SubElement(self.arGridCollection,"Grid",{"Name":gridName,"GridType":"Uniform"})
-                self.arTime = SubElement(self.arGrid,"Time",{"Value":"%e" % (t,),"Name":str(tCount)})
-                topology    = SubElement(self.arGrid,"Topology",
-                                         {"Type":Xdmf_ElementTopology,
-                                          "NumberOfElements":str(Xdmf_NumberOfElements)})
-                elements    = SubElement(topology,"DataItem",
-                                         {"Format":ar.dataItemFormat,
-                                          **xdmf_number_type(mesh.elementNodesArray),
-                                          "Dimensions":"%i %i" % (Xdmf_NumberOfElements,Xdmf_NodesPerElement)})
-                geometry    = SubElement(self.arGrid,"Geometry",{"Type":"XYZ"})
-                nodes       = SubElement(geometry,"DataItem",
-                                         {"Format":ar.dataItemFormat,
-                                          **xdmf_number_type(mesh.nodeArray),
-                                          "Dimensions":"%i %i" % (Xdmf_NumberOfNodes,3)})
-                if ar.hdfFile is not None:
-                    elements.text = ar.hdfFilename+":/elements"+str(ar.rank)+spaceSuffix+str(tCount)
-                    nodes.text    = ar.hdfFilename+":/nodes"+str(ar.rank)+spaceSuffix+str(tCount)
-                    if init or meshChanged:
-                        #c0p1 mapping for now
-                        ar.create_dataset_async('elements'+str(ar.rank)+spaceSuffix+str(tCount), data = mesh.elementNodesArray)
-                        ar.create_dataset_async('nodes'+str(ar.rank)+spaceSuffix+str(tCount), data = mesh.nodeArray)
+                elements_dataset = 'elements'+str(ar.rank)+spaceSuffix+str(tCount)
+                nodes_dataset    = 'nodes'+str(ar.rank)+spaceSuffix+str(tCount)
+            elements_stem = 'elements'+spaceSuffix+str(tCount)
+            nodes_stem    = 'nodes'+spaceSuffix+str(tCount)
+
+            self.arGrid = ar.write_grid(self.arGridCollection, gridName, t, tCount)
+            ar.write_topology(self.arGrid, Xdmf_ElementTopology, Xdmf_NumberOfElements,
+                              [Xdmf_NumberOfElements, Xdmf_NodesPerElement],
+                              elements_dataset, elements_stem,
+                              dtype=mesh.elementNodesArray)
+            ar.write_geometry(self.arGrid, [Xdmf_NumberOfNodes, 3],
+                              nodes_dataset, nodes_stem, dtype=mesh.nodeArray)
+            if init or meshChanged:
+                if ar.global_sync:
+                    #c0p1 mapping for now
+                    ar.create_dataset_sync(elements_dataset,
+                                           offsets = mesh.globalMesh.elementOffsets_subdomain_owned,
+                                           data = mesh.elementNodesArray)
+                    ar.create_dataset_sync(nodes_dataset,
+                                           offsets = mesh.globalMesh.nodeOffsets_subdomain_owned,
+                                           data = mesh.nodeArray)
+                elif ar.hdfFile is not None:
+                    ar.create_dataset_async(elements_dataset, data = mesh.elementNodesArray)
+                    ar.create_dataset_async(nodes_dataset, data = mesh.nodeArray)
                 else:
-                    SubElement(elements,"xi:include",{"parse":"text","href":"./"+ar.textDataDir+"/elements"+spaceSuffix+str(tCount)+".txt"})
-                    SubElement(nodes,"xi:include",{"parse":"text","href":"./"+ar.textDataDir+"/nodes"+spaceSuffix+str(tCount)+".txt"})
-                    if init or meshChanged:
-                        numpy.savetxt(ar.textDataDir+"/elements"+spaceSuffix+str(tCount)+".txt",mesh.elementNodesArray,fmt='%d')
-                        numpy.savetxt(ar.textDataDir+"/nodes"+spaceSuffix+str(tCount)+".txt",mesh.nodeArray)
-
-                    #
-                #hdfile
-            #need to write a grid
+                    numpy.savetxt(ar.textDataDir+"/"+elements_stem+".txt",mesh.elementNodesArray,fmt='%d')
+                    numpy.savetxt(ar.textDataDir+"/"+nodes_stem+".txt",mesh.nodeArray)
         return self.arGrid
     #def
     def writeFunctionXdmf_P1Bubble(self,ar,u,tCount=0,init=True):
         #this read ar.sync_global, which does not exist -- the attribute
         #is global_sync -- so every call raised AttributeError
+        name = u.name.replace(' ','_')
+        mesh = u.femSpace.mesh
+        #just write out nodal part right now
         if ar.global_sync:
-            #just write out nodal part right now
-            Xdmf_NumberOfElements = u.femSpace.mesh.globalMesh.nElements_global
-            Xdmf_NumberOfNodes    = u.femSpace.mesh.globalMesh.nNodes_global
-            name = u.name.replace(' ','_')
-
-            #if writing as dgp1
-            attribute = SubElement(self.arGrid,"Attribute",{"Name":name,
-                                                            "AttributeType":"Scalar",
-                                                            "Center":"Node"})
-            values    = SubElement(attribute,"DataItem",
-                                   {"Format":ar.dataItemFormat,
-                                    **xdmf_number_type(u.dof),
-                                    "Dimensions":"%i" % (Xdmf_NumberOfNodes,)})
-            if ar.hdfFile is not None:
-                values.text = ar.hdfFilename+":/"+name+"_t"+str(tCount)
-                ar.create_dataset_sync(name+"_t"+str(tCount),
-                                       offsets = u.femSpace.mesh.nodeOffsets_subdomain_owned,
-                                       data = u.dof[0:u.femSpace.mesh.nNodes_owned])
-            else:
-                assert False, "global_sync not implemented for text heavy data"
+            ar.write_field(self.arGrid, name, u.dof, tCount,
+                           dimensions=[mesh.globalMesh.nNodes_global],
+                           sync_offsets=mesh.globalMesh.nodeOffsets_subdomain_owned,
+                           sync_data=u.dof[0:mesh.nNodes_owned])
         else:
-            #just write out nodal part right now
-            Xdmf_NumberOfElements = u.femSpace.mesh.nElements_global
-            Xdmf_NumberOfNodes    = u.femSpace.mesh.nNodes_global
-            name = u.name.replace(' ','_')
-
-            #if writing as dgp1
-            attribute = SubElement(self.arGrid,"Attribute",{"Name":name,
-                                                            "AttributeType":"Scalar",
-                                                            "Center":"Node"})
-            values    = SubElement(attribute,"DataItem",
-                                   {"Format":ar.dataItemFormat,
-                                    **xdmf_number_type(u.dof),
-                                    "Dimensions":"%i" % (Xdmf_NumberOfNodes,)})
-            if ar.hdfFile is not None:
-                values.text = ar.hdfFilename+":/"+name+"_p"+str(ar.rank)+"_t"+str(tCount)
-                ar.create_dataset_async(name+"_p"+str(ar.rank)+"_t"+str(tCount), data = u.dof[0:Xdmf_NumberOfNodes])
-            else:
-                numpy.savetxt(ar.textDataDir+"/"+name+str(tCount)+".txt",u.dof[0:Xdmf_NumberOfNodes])
-                SubElement(values,"xi:include",{"parse":"text","href":"./"+ar.textDataDir+"/"+name+str(tCount)+".txt"})
+            ar.write_field(self.arGrid, name, u.dof[0:mesh.nNodes_global], tCount,
+                           dimensions=[mesh.nNodes_global])
 
     def writeFunctionXdmf_C0P2Lagrange(self,ar,u,tCount=0,init=True):
         #the text branch here read `assert Fasle` -- a typo for False, so it
@@ -2383,54 +2015,21 @@ def writeFunctionXdmf_C0P2Lagrange(self,ar,u,tCount=0,init=True):
             ar.write_field(self.arGrid, u.name, u.dof, tCount)
 
     def writeVectorFunctionXdmf_P1Bubble(self,ar,uList,components,vectorName,spaceSuffix,tCount=0,init=True):
+        mesh = uList[components[0]].femSpace.mesh
         if ar.global_sync:
-            nDOF_global = uList[components[0]].femSpace.mesh.globalMesh.nNodes_global
-            nDOF_local = (uList[components[0]].femSpace.mesh.globalMesh.nodeOffsets_subdomain_owned[ar.rank+1] -
-                          uList[components[0]].femSpace.mesh.globalMesh.nodeOffsets_subdomain_owned[ar.rank])
-            attribute = SubElement(self.arGrid,"Attribute",{"Name":vectorName,
-                                                            "AttributeType":"Vector",
-                                                            "Center":"Node"})
-            u_dof = uList[components[0]].dof[0:nDOF_local]
-            if len(components) < 2:
-                v_dof = numpy.zeros((nDOF_global,),dtype='d')
-            else:
-                v_dof = uList[components[1]].dof[0:nDOF_local]
-            if len(components) < 3:
-                w_dof = numpy.zeros((nDOF_global,),dtype='d')
-            else:
-                w_dof = uList[components[2]].dof[0:nDOF_local]
-            velocity = numpy.column_stack((u_dof,v_dof,w_dof))
-            values    = SubElement(attribute,"DataItem",
-                                   {"Format":ar.dataItemFormat,
-                                    **xdmf_number_type(velocity),
-                                    "Dimensions":"%i %i" % (nDOF_global,3)})
-            if ar.hdfFile is not None:
-                values.text = ar.hdfFilename+":/"+vectorName+"_t"+str(tCount)
-                ar.create_dataset_sync(vectorName+"_t"+str(tCount),
-                                       offsets = uList[components[0]].femSpace.mesh.globalMesh.nodeOffsets_subdomain_owned,
-                                       data = velocity)
+            nDOF_global = mesh.globalMesh.nNodes_global
+            nDOF_local = (mesh.globalMesh.nodeOffsets_subdomain_owned[ar.rank+1] -
+                          mesh.globalMesh.nodeOffsets_subdomain_owned[ar.rank])
+            velocity = self._nodalVector(uList,components,n_rows=nDOF_local)
+            ar.write_field(self.arGrid, vectorName, velocity, tCount, rank="Vector",
+                           dimensions=[nDOF_global,3],
+                           sync_offsets=mesh.globalMesh.nodeOffsets_subdomain_owned,
+                           sync_data=velocity)
         else:
-            nDOF_global = uList[components[0]].femSpace.mesh.nNodes_global
-            attribute = SubElement(self.arGrid,"Attribute",{"Name":vectorName,
-                                                            "AttributeType":"Vector",
-                                                            "Center":"Node"})
-            u_dof = uList[components[0]].dof[0:nDOF_global]
-            if len(components) < 2:
-                v_dof = numpy.zeros((nDOF_global,),dtype='d')
-            else:
-                v_dof = uList[components[1]].dof[0:nDOF_global]
-            if len(components) < 3:
-                w_dof = numpy.zeros((nDOF_global,),dtype='d')
-            else:
-                w_dof = uList[components[2]].dof[0:nDOF_global]
-            velocity = numpy.column_stack((u_dof,v_dof,w_dof))
-            values    = SubElement(attribute,"DataItem",
-                                   {"Format":ar.dataItemFormat,
-                                    **xdmf_number_type(velocity),
-                                    "Dimensions":"%i %i" % (nDOF_global,3)})
-            if ar.hdfFile is not None:
-                values.text = ar.hdfFilename+":/"+vectorName+"_p"+str(ar.rank)+"_t"+str(tCount)
-                ar.create_dataset_async(vectorName+"_p"+str(ar.rank)+"_t"+str(tCount), data = velocity)
+            nDOF_global = mesh.nNodes_global
+            velocity = self._nodalVector(uList,components,n_rows=nDOF_global)
+            ar.write_field(self.arGrid, vectorName, velocity, tCount, rank="Vector",
+                           dimensions=[nDOF_global,3])
 
     def writeMeshXdmf_particles(self,ar,mesh,spaceDim,x,t=0.0,
                                 init=False,meshChanged=False,arGrid=None,tCount=0,
@@ -2444,7 +2043,7 @@ def writeMeshXdmf_particles(self,ar,mesh,spaceDim,x,t=0.0,
 
         gridName = self.setGridCollectionAndGridElements(init,ar,arGrid,t,spaceSuffix)
         nPoints = numpy.cumprod(x.shape)[-2]
-        if self.arGrid is None or self.arTime.get('Value') != "{0:e}".format(t):
+        if self._needsGrid(t):
             Xdmf_NodesPerElement = 1
             #the dataset names carry the rank before the space suffix while the
             #sidecar names omit it; both predate write_field's convention
@@ -2453,8 +2052,7 @@ def writeMeshXdmf_particles(self,ar,mesh,spaceDim,x,t=0.0,
             elements_stem    = 'elements'+spaceSuffix+str(tCount)
             nodes_stem       = 'nodes'+spaceSuffix+str(tCount)
 
-            self.arGrid, self.arTime = ar.write_grid(self.arGridCollection,
-                                                     gridName, t, tCount)
+            self.arGrid = ar.write_grid(self.arGridCollection, gridName, t, tCount)
             ar.write_topology(self.arGrid, "Polyvertex", nPoints,
                               [nPoints, Xdmf_NodesPerElement],
                               elements_dataset, elements_stem,
@@ -2497,130 +2095,32 @@ def writeMeshXdmf_LowestOrderMixed(self,ar,mesh,spaceDim,t=0.0,init=False,meshCh
         mesh.writeMeshXdmf(ar,"Spatial_Domain",t,init,meshChanged,tCount=tCount)
         #now try to write out a mesh that matches RT0 velocity as dgp1 lagrange
         gridName = self.setGridCollectionAndGridElements(init,ar,arGrid,t,spaceSuffix)
-
-        if self.arGrid is None or self.arTime.get('Value') != "{0:e}".format(t):
-            if spaceDim == 1:
-                Xdmf_ElementTopology = "Polyline"
-            elif spaceDim == 2:
-                Xdmf_ElementTopology = "Triangle"
-            elif spaceDim == 3:
-                Xdmf_ElementTopology = "Tetrahedron"
-            Xdmf_NodesPerElement = spaceDim+1
-            if ar.global_sync:
-                Xdmf_NumberOfElements= mesh.globalMesh.nElements_global
-                self.arGrid = SubElement(self.arGridCollection,"Grid",{"Name":gridName,"GridType":"Uniform"})
-                self.arTime = SubElement(self.arGrid,"Time",{"Value":"%e" % (t,),"Name":str(tCount)})
-                topology    = SubElement(self.arGrid,"Topology",
-                                         {"Type":Xdmf_ElementTopology,
-                                          "NumberOfElements":str(Xdmf_NumberOfElements)})
-                elements    = SubElement(topology,"DataItem",
-                                         {"Format":ar.dataItemFormat,
-                                          **xdmf_number_type(GENERATED_INDEX_DTYPE),
-                                          "Dimensions":"%i %i" % (Xdmf_NumberOfElements,Xdmf_NodesPerElement)})
-                geometry    = SubElement(self.arGrid,"Geometry",{"Type":"XYZ"})
-                nodes       = SubElement(geometry,"DataItem",
-                                         {"Format":ar.dataItemFormat,
-                                          **xdmf_number_type(mesh.nodeArray),
-                                          "Dimensions":"%i %i" % (Xdmf_NumberOfElements*Xdmf_NodesPerElement,3)})
-                if ar.hdfFile is not None:
-                    elements.text = ar.hdfFilename+":/elements"+spaceSuffix+str(tCount)
-                    nodes.text    = ar.hdfFilename+":/nodes"+spaceSuffix+str(tCount)
-                    if init or meshChanged:
-                        #simple dg l2g mapping
-                        dg_l2g = numpy.arange(Xdmf_NumberOfElements*Xdmf_NodesPerElement,dtype=GENERATED_INDEX_DTYPE).reshape((mesh.nElements_global,
-                                                                                                             Xdmf_NodesPerElement))
-                        ar.create_dataset_sync('elements'+spaceSuffix+str(tCount),
-                                               offsets = mesh.globalMesh.elementOffsets_subdomain_owned,
-                                               data = dg_l2g)
-                        
-                    dgnodes = numpy.reshape(mesh.nodeArray[mesh.elementNodesArray[:mesh.nElements_owned]],(mesh.nElements_owned*Xdmf_NodesPerElement,3))
-                    ar.create_dataset_sync('nodes'+spaceSuffix+str(tCount),
-                                           offsets = mesh.globalMesh.elementOffsets_subdomain_owned*Xdmf_NodesPerElement,
-                                           data = dgnodes)
-                else:
-                    assert False, "global_sync not implemented for text heavy data"
-            else:
-                Xdmf_NumberOfElements= mesh.nElements_global
-                self.arGrid = SubElement(self.arGridCollection,"Grid",{"Name":gridName,"GridType":"Uniform"})
-                self.arTime = SubElement(self.arGrid,"Time",{"Value":"%e" % (t,),"Name":str(tCount)})
-                topology    = SubElement(self.arGrid,"Topology",
-                                         {"Type":Xdmf_ElementTopology,
-                                          "NumberOfElements":str(Xdmf_NumberOfElements)})
-                elements    = SubElement(topology,"DataItem",
-                                         {"Format":ar.dataItemFormat,
-                                          **xdmf_number_type(GENERATED_INDEX_DTYPE),
-                                          "Dimensions":"%i %i" % (Xdmf_NumberOfElements,Xdmf_NodesPerElement)})
-                geometry    = SubElement(self.arGrid,"Geometry",{"Type":"XYZ"})
-                nodes       = SubElement(geometry,"DataItem",
-                                         {"Format":ar.dataItemFormat,
-                                          **xdmf_number_type(mesh.nodeArray),
-                                          "Dimensions":"%i %i" % (Xdmf_NumberOfElements*Xdmf_NodesPerElement,3)})
-                if ar.hdfFile is not None:
-                    elements.text = ar.hdfFilename+":/elements"+str(ar.rank)+spaceSuffix+str(tCount)
-                    nodes.text    = ar.hdfFilename+":/nodes"+str(ar.rank)+spaceSuffix+str(tCount)
-                    if init or meshChanged:
-                        #simple dg l2g mapping
-                        dg_l2g = numpy.arange(Xdmf_NumberOfElements*Xdmf_NodesPerElement,dtype=GENERATED_INDEX_DTYPE).reshape((Xdmf_NumberOfElements,Xdmf_NodesPerElement))
-                        ar.create_dataset_async('elements'+str(ar.rank)+spaceSuffix+str(tCount), data = dg_l2g)
-                    
-                    dgnodes = numpy.reshape(mesh.nodeArray[mesh.elementNodesArray],(Xdmf_NumberOfElements*Xdmf_NodesPerElement,3))
-                    ar.create_dataset_async('nodes'+str(ar.rank)+spaceSuffix+str(tCount), data = dgnodes)
-                else:
-                    SubElement(elements,"xi:include",{"parse":"text","href":"./"+ar.textDataDir+"/elements"+spaceSuffix+str(tCount)+".txt"})
-                    SubElement(nodes,"xi:include",{"parse":"text","href":"./"+ar.textDataDir+"/nodes"+spaceSuffix+str(tCount)+".txt"})
-                    if init or meshChanged:
-                        dg_l2g = numpy.arange(Xdmf_NumberOfElements*Xdmf_NodesPerElement,dtype=GENERATED_INDEX_DTYPE).reshape((Xdmf_NumberOfElements,Xdmf_NodesPerElement))
-                        numpy.savetxt(ar.textDataDir+"/elements"+spaceSuffix+str(tCount)+".txt",dg_l2g,fmt='%d')
-
-                        dgnodes = numpy.reshape(mesh.nodeArray[mesh.elementNodesArray],(Xdmf_NumberOfElements*Xdmf_NodesPerElement,3))
-                        numpy.savetxt(ar.textDataDir+"/nodes"+spaceSuffix+str(tCount)+".txt",dgnodes)
-
-                    #
-                #hdfile
-            #need to write a grid
+        self.mesh = mesh
+        if self._needsGrid(t):
+            #this writer rewrote the nodes on every pass of the HDF5 path
+            self._writeDiscontinuousVertexMesh(ar,gridName,mesh,spaceDim,t,tCount,
+                                               spaceSuffix,init,meshChanged,
+                                               nodes_every_pass=True)
         return self.arGrid
     #def
     def writeVectorFunctionXdmf_LowestOrderMixed(self,ar,u,tCount=0,init=True,spaceSuffix="_RT0",name="velocity"):
+        Xdmf_NumberOfComponents = u.shape[2]
+        Xdmf_StorageDim = 3
+        n_points = u.shape[1]
         if ar.global_sync:
-            Xdmf_NodesGlobal = self.mesh.globalMesh.nElements_global*u.shape[1]
-            Xdmf_NumberOfComponents = u.shape[2]
-            Xdmf_StorageDim = 3
-
-            #if writing as dgp1
-            attribute = SubElement(self.arGrid,"Attribute",{"Name":name,
-                                                            "AttributeType":"Vector",
-                                                            "Center":"Node"})
-            values    = SubElement(attribute,"DataItem",
-                                   {"Format":ar.dataItemFormat,
-                                    **xdmf_number_type(u),
-                                    "Dimensions":"%i %i" % (Xdmf_NodesGlobal,Xdmf_StorageDim)})#force 3d vector since points 3d
-            tmp = numpy.zeros((u.shape[0]*u.shape[1],Xdmf_StorageDim),u.dtype)
-            tmp[:,:Xdmf_NumberOfComponents]=numpy.reshape(u.flat,(u.shape[0]*u.shape[1],Xdmf_NumberOfComponents))
-
-            if ar.hdfFile is not None:
-                values.text = ar.hdfFilename+":/"+name+"_t"+str(tCount)
-                ar.create_dataset_async(name+"_t"+str(tCount), data = tmp)
-            else:
-                assert False, "global_sync not implemented for text heavy data"
+            #this rank contributes the elements it owns; the previous version
+            #wrote the whole local array under the global name, which only
+            #agrees on one rank
+            n_owned = self.mesh.nElements_owned
+            tmp = numpy.zeros((n_owned*n_points,Xdmf_StorageDim),u.dtype)
+            tmp[:,:Xdmf_NumberOfComponents]=numpy.reshape(u[:n_owned].flat,(n_owned*n_points,Xdmf_NumberOfComponents))
+            ar.write_field(self.arGrid, name, tmp, tCount, rank="Vector",
+                           dimensions=[self.mesh.globalMesh.nElements_global*n_points,Xdmf_StorageDim],
+                           sync_offsets=self.mesh.globalMesh.elementOffsets_subdomain_owned*n_points,
+                           sync_data=tmp)
         else:
-            Xdmf_NodesGlobal = u.shape[0]*u.shape[1]
-            Xdmf_NumberOfComponents = u.shape[2]
-            Xdmf_StorageDim = 3
-
-            #if writing as dgp1
-            attribute = SubElement(self.arGrid,"Attribute",{"Name":name,
-                                                            "AttributeType":"Vector",
-                                                            "Center":"Node"})
-            values    = SubElement(attribute,"DataItem",
-                                   {"Format":ar.dataItemFormat,
-                                    **xdmf_number_type(u),
-                                    "Dimensions":"%i %i" % (Xdmf_NodesGlobal,Xdmf_StorageDim)})#force 3d vector since points 3d
+            Xdmf_NodesGlobal = u.shape[0]*n_points
             tmp = numpy.zeros((Xdmf_NodesGlobal,Xdmf_StorageDim),u.dtype)
             tmp[:,:Xdmf_NumberOfComponents]=numpy.reshape(u.flat,(Xdmf_NodesGlobal,Xdmf_NumberOfComponents))
-
-            if ar.hdfFile is not None:
-                values.text = ar.hdfFilename+":/"+name+"_p"+str(ar.rank)+"_t"+str(tCount)
-                ar.create_dataset_async(name+"_p"+str(ar.rank)+"_t"+str(tCount), data = tmp)
-            else:
-                numpy.savetxt(ar.textDataDir+"/"+name+str(tCount)+".txt",tmp)
-                SubElement(values,"xi:include",{"parse":"text","href":"./"+ar.textDataDir+"/"+name+str(tCount)+".txt"})
\ No newline at end of file
+            ar.write_field(self.arGrid, name, tmp, tCount, rank="Vector",
+                           dimensions=[Xdmf_NodesGlobal,Xdmf_StorageDim])
diff --git a/proteus/AuxiliaryVariables.py b/proteus/AuxiliaryVariables.py
index 8c2a2b7da0..0f03e3dcf1 100644
--- a/proteus/AuxiliaryVariables.py
+++ b/proteus/AuxiliaryVariables.py
@@ -7,7 +7,6 @@
 import numpy
 from . import Viewers
 from . import Archiver
-from xml.etree.ElementTree import *
 
 from . import Profiling
 from .Profiling import logEvent
@@ -31,21 +30,19 @@ def attachModel(self,model,ar):
     def attachAuxiliaryVariables(self,avDict):
         pass
     def writeScalarXdmf(self,u,name):
-        #get the current spatial grid
-        for child in self.ar.domain:
-            if child.tag == "Grid" and child.attrib["Name"] == "Mesh Spatial_Domain":
-                self.arGridCollection = child
-        assert self.arGridCollection is not None
-        self.arGrid = self.arGridCollection[-1]
-        #write the attribute
-        attribute = SubElement(self.arGrid,"Attribute",{"Name":name,
-                                                        "AttributeType":"Scalar",
-                                                        "Center":"Grid"})
-        values    = SubElement(attribute,"DataItem",
-                               {"Format":"XML",
-                                **Archiver.xdmf_number_type(numpy.asarray(u)),
-                                "Dimensions":"1"})
-        values.text = str(u)
+        """Record one number for the step, inline in the archive document."""
+        from ymf.archive import attribute, data_item_for
+        #the current spatial grid
+        self.arGridCollection = self.ar.domain.get("Mesh Spatial_Domain")
+        assert self.arGridCollection is not None, \
+            "no Mesh Spatial_Domain collection; the mesh is archived first"
+        self.arGrid = self.arGridCollection.grid
+        assert self.arGrid is not None, \
+            "Mesh Spatial_Domain has no grid for this step yet"
+        #a grid-centred attribute holding the value itself
+        item = data_item_for(numpy.asarray(u).reshape((1,)), inline=True)
+        self.arGrid["attributes"].append(
+            attribute(name, item, center="Grid", attribute_type="Scalar"))
     def calculate_init(self):
         self.calculate()
     def calculate(self):
diff --git a/proteus/FemTools.py b/proteus/FemTools.py
index 819659df8e..1d67788917 100644
--- a/proteus/FemTools.py
+++ b/proteus/FemTools.py
@@ -4176,8 +4176,6 @@ def readFunctionXdmf(self,ar,u,tCount=0):
             u.dof[permute] = dataset[numbering[permute].tolist()]
         else:
             u.dof[:] = dataset[:]
-            #numpy.savetxt(ar.textDataDir+"/"+"{0:s}{0:d}".format(u.name, tCount)+".txt",u.dof)
-            #SubElement(values,"xi:include",{"parse":"text","href":"./"+ar.textDataDir+"/"+"{0:s}{0:d}".format(u.name, tCount)+".txt"})
     def writeVectorFunctionXdmf(self,ar,uList,components,vectorName,tCount=0,init=True):
         #XDMF vectors are always 3-component, so a 1D or 2D problem pads the
         #components it doesn't have with zeros. This was duplicated verbatim
diff --git a/proteus/MeshTools.py b/proteus/MeshTools.py
index b48a10c548..5ef8a88b87 100644
--- a/proteus/MeshTools.py
+++ b/proteus/MeshTools.py
@@ -10,7 +10,6 @@
 import inspect
 import h5py
 import os
-from xml.etree import ElementTree as ET
 from .Archiver import *
 from .LinearAlgebraTools import ParVec_petsc4py
 from .Profiling import logEvent,memory, memHardLimit
@@ -999,231 +998,135 @@ def partitionMeshFromFiles(self,filebase,base,nLayersOfOverlap=1,parallelPartiti
         # logEvent(memory("Without global mesh","Mesh"),level=1)
         # comm.endSequential()
     def writeMeshXdmf(self,ar,name='',t=0.0,init=False,meshChanged=False,Xdmf_ElementTopology="Triangle",tCount=0, EB=False):
+        """Describe this mesh in the archive at time ``t``; with ``EB`` its element boundaries too.
+
+        The grid -- a ymf grid dict in the collection ``Mesh <name>`` (and
+        ``EBMesh <name>``) -- is built on every pass through the archive's
+        builders; the arrays behind it are written only when ``init or
+        meshChanged``. The material types are written on every pass, as
+        they were; the local-to-global maps, which ``collect.py`` and a
+        hot start from a global archive read, only with the arrays.
+        """
         if self.arGridCollection is not None:
             init = False
         elif not init:
-            grids = ar.domain.findall("Grid")
-            self.arGridCollection = grids[0]
+            self.arGridCollection = ar.domain[0]
             if EB:
-                assert(len(grids) > 1)
-                self.arEBGridCollection = grids[1]
+                assert(len(ar.domain) > 1)
+                self.arEBGridCollection = ar.domain[1]
         if init:
-            self.arGridCollection = SubElement(ar.domain,"Grid",{"Name":"Mesh "+name,
-                                                               "GridType":"Collection",
-                                                               "CollectionType":"Temporal"})
+            self.arGridCollection = ar.domain.collection("Mesh "+name)
             if EB:
-                self.arEBGridCollection = SubElement(ar.domain,"Grid",{"Name":"EBMesh "+name,
-                                                                       "GridType":"Collection",
-                                                                       "CollectionType":"Temporal"})
-        if self.arGrid is None or self.arTime.get('Value') != "{0:e}".format(t):
+                self.arEBGridCollection = ar.domain.collection("EBMesh "+name)
+        if not self.arGridCollection.is_current(t):
+            rank = ar.comm.rank()
             #
             #topology and geometry
             #
             if ar.global_sync:
-                self.arGrid = SubElement(self.arGridCollection,"Grid",{"GridType":"Uniform"})
-                self.arTime = SubElement(self.arGrid,"Time",{"Value":"%e" % (t,),"Name":"%i" % (tCount,)})
-                topology = SubElement(self.arGrid,"Topology",
-                                      {"Type":Xdmf_ElementTopology,
-                                       "NumberOfElements":"%i" % (self.globalMesh.nElements_global,)})
-                elements = SubElement(topology,"DataItem",
-                                      {"Format":ar.dataItemFormat,
-                                       **xdmf_number_type(self.globalMesh.nodeNumbering_subdomain2global),
-                                       "Dimensions":"%i %i" % (self.globalMesh.nElements_global,self.nNodes_element)})
-                geometry = SubElement(self.arGrid,"Geometry",{"Type":"XYZ"})
-                nodes    = SubElement(geometry,"DataItem",
-                                      {"Format":ar.dataItemFormat,
-                                       **xdmf_number_type(self.nodeArray),
-                                       "Dimensions":"%i %i" % (self.globalMesh.nNodes_global,3)})
-                if ar.hdfFile is not None:
-                    elements.text = ar.hdfFilename+":/elements"+name+str(tCount)
-                    nodes.text = ar.hdfFilename+":/nodes"+name+str(tCount)
-                    if init or meshChanged:
-                        ar.create_dataset_sync('elements'+name+str(tCount),
-                                               offsets=self.globalMesh.elementOffsets_subdomain_owned,
-                                               data=self.globalMesh.nodeNumbering_subdomain2global[self.elementNodesArray[:self.nElements_owned]])
-                        ar.create_dataset_sync('nodes'+name+str(tCount),
-                                               offsets=self.globalMesh.nodeOffsets_subdomain_owned,
-                                               data=self.nodeArray[:self.nNodes_owned])
-                else:
-                    assert False, "global_sync not  supported with text heavy data"
+                assert ar.hdfFile is not None, "global_sync not  supported with text heavy data"
+                n_elements       = self.globalMesh.nElements_global
+                nodes_dims       = [self.globalMesh.nNodes_global, 3]
+                elements_dataset = 'elements'+name+str(tCount)
+                nodes_dataset    = 'nodes'+name+str(tCount)
+                elements_dtype   = self.globalMesh.nodeNumbering_subdomain2global
             else:
-                self.arGrid = SubElement(self.arGridCollection,"Grid",{"GridType":"Uniform"})
-                self.arTime = SubElement(self.arGrid,"Time",{"Value":"%e" % (t,),"Name":"%i" % (tCount,)})
-                topology = SubElement(self.arGrid,"Topology",
-                                      {"Type":Xdmf_ElementTopology,
-                                       "NumberOfElements":"%i" % (self.nElements_owned,)})
-                elements = SubElement(topology,"DataItem",
-                                      {"Format":ar.dataItemFormat,
-                                       **xdmf_number_type(self.elementNodesArray),
-                                       "Dimensions":"%i %i" % (self.nElements_owned,self.nNodes_element)})
-                geometry = SubElement(self.arGrid,"Geometry",{"Type":"XYZ"})
-                nodes    = SubElement(geometry,"DataItem",
-                                      {"Format":ar.dataItemFormat,
-                                       **xdmf_number_type(self.nodeArray),
-                                       "Dimensions":"%i %i" % (self.nNodes_global,3)})
-                if ar.hdfFile is not None:
-                    elements.text = ar.hdfFilename+":/elements"+str(ar.comm.rank())+name+str(tCount)
-                    nodes.text = ar.hdfFilename+":/nodes"+str(ar.comm.rank())+name+str(tCount)
-                    if init or meshChanged:
-                        ar.create_dataset_async('elements'+str(ar.comm.rank())+name+str(tCount),data=self.elementNodesArray[:self.nElements_owned])
-                        ar.create_dataset_async('nodes'+str(ar.comm.rank())+name+str(tCount),data=self.nodeArray)
+                n_elements       = self.nElements_owned
+                nodes_dims       = [self.nNodes_global, 3]
+                elements_dataset = 'elements'+str(rank)+name+str(tCount)
+                nodes_dataset    = 'nodes'+str(rank)+name+str(tCount)
+                elements_dtype   = self.elementNodesArray
+            self.arGrid = ar.write_grid(self.arGridCollection, None, t, tCount)
+            ar.write_topology(self.arGrid, Xdmf_ElementTopology, n_elements,
+                              [n_elements, self.nNodes_element],
+                              elements_dataset, 'elements'+name, dtype=elements_dtype)
+            ar.write_geometry(self.arGrid, nodes_dims, nodes_dataset, 'nodes'+name,
+                              dtype=self.nodeArray)
+            if init or meshChanged:
+                if ar.global_sync:
+                    ar.create_dataset_sync(elements_dataset,
+                                           offsets=self.globalMesh.elementOffsets_subdomain_owned,
+                                           data=self.globalMesh.nodeNumbering_subdomain2global[self.elementNodesArray[:self.nElements_owned]])
+                    ar.create_dataset_sync(nodes_dataset,
+                                           offsets=self.globalMesh.nodeOffsets_subdomain_owned,
+                                           data=self.nodeArray[:self.nNodes_owned])
+                elif ar.hdfFile is not None:
+                    ar.create_dataset_async(elements_dataset,data=self.elementNodesArray[:self.nElements_owned])
+                    ar.create_dataset_async(nodes_dataset,data=self.nodeArray)
                 else:
-                    SubElement(elements,"xi:include",{"parse":"text","href":"./"+ar.textDataDir+"/elements"+name+".txt"})
-                    SubElement(nodes,"xi:include",{"parse":"text","href":"./"+ar.textDataDir+"/nodes"+name+".txt"})
-                    if init or meshChanged:
-                        numpy.savetxt(ar.textDataDir+"/elements"+name+".txt",self.elementNodesArray[:self.nElements_owned],fmt='%d')
-                        numpy.savetxt(ar.textDataDir+"/nodes"+name+".txt",self.nodeArray)
+                    numpy.savetxt(ar.textDataDir+"/elements"+name+".txt",self.elementNodesArray[:self.nElements_owned],fmt='%d')
+                    numpy.savetxt(ar.textDataDir+"/nodes"+name+".txt",self.nodeArray)
             #
             #element boundary topology and geometry
             #
             if EB:
-                self.arEBGrid = SubElement(self.arEBGridCollection,"Grid",{"GridType":"Uniform"})
-                self.arEBTime = SubElement(self.arEBGrid,"Time",{"Value":"%e" % (t,),"Name":"%i" % (tCount,)})
                 Xdmf_ElementEBTopology = "Triangle" #cek hack
-                ebtopology = SubElement(self.arEBGrid,"Topology",
-                                    {"Type":Xdmf_ElementEBTopology,
-                                     "NumberOfElements":"%i" % (self.nElementBoundaries_global,)})
-                ebelements = SubElement(ebtopology,"DataItem",
-                                    {"Format":ar.dataItemFormat,
-                                     **xdmf_number_type(self.elementBoundaryNodesArray),
-                                     "Dimensions":"%i %i" % (self.nElementBoundaries_global,self.nNodes_elementBoundary)})
-                ebgeometry = SubElement(self.arEBGrid,"Geometry",{"Type":"XYZ"})
-                ebnodes    = SubElement(ebgeometry,"DataItem",
-                                    {"Format":ar.dataItemFormat,
-                                     **xdmf_number_type(self.nodeArray),
-                                     "Dimensions":"%i %i" % (self.nNodes_global,3)})
-                if ar.hdfFile is not None:
-                    ebelements.text = ar.hdfFilename+":/elementBoundaries"+str(ar.comm.rank())+name+str(tCount)
-                    ebnodes.text = ar.hdfFilename+":/nodes"+str(ar.comm.rank())+name+str(tCount)
-                    if init or meshChanged:
-                        ar.create_dataset_async('elementBoundaries'+str(ar.comm.rank())+name+str(tCount), data = self.elementBoundaryNodesArray)
-                        #ar.create_dataset_async('nodes'+`ar.comm.rank()`+name+`tCount`, data = self.nodeArray)
-                else:
-                    SubElement(ebelements,"xi:include",{"parse":"text","href":"./"+ar.textDataDir+"/elementBoundaries"+name+".txt"})
-                    SubElement(ebnodes,"xi:include",{"parse":"text","href":"./"+ar.textDataDir+"/nodes"+name+".txt"})
-                    if init or meshChanged:
+                self.arEBGrid = ar.write_grid(self.arEBGridCollection, None, t, tCount)
+                #the boundary grid shares the mesh's nodes dataset
+                ar.write_topology(self.arEBGrid, Xdmf_ElementEBTopology, self.nElementBoundaries_global,
+                                  [self.nElementBoundaries_global, self.nNodes_elementBoundary],
+                                  'elementBoundaries'+str(rank)+name+str(tCount),
+                                  'elementBoundaries'+name, dtype=self.elementBoundaryNodesArray)
+                ar.write_geometry(self.arEBGrid, [self.nNodes_global, 3],
+                                  'nodes'+str(rank)+name+str(tCount), 'nodes'+name,
+                                  dtype=self.nodeArray)
+                if init or meshChanged:
+                    if ar.hdfFile is not None:
+                        ar.create_dataset_async('elementBoundaries'+str(rank)+name+str(tCount), data = self.elementBoundaryNodesArray)
+                    else:
                         np.savetxt(ar.textDataDir+"/elementBoundaries"+name+".txt",self.elementBoundaryNodesArray,fmt='%d')
 
             # Add the local->global index maps for collect.py and for
-            # reverse mapping in hotstarts from a global XDMF file.
+            # reverse mapping in hotstarts from a global archive.
             if self.globalMesh is not None and not ar.global_sync:
-                nodeMapAtt = SubElement(self.arGrid,"Attribute",
-                                        {"Name":"NodeMapL2G",
-                                         "AttributeType":"Scalar",
-                                         "Center":"Node"})
-                nodeMap = SubElement(nodeMapAtt,"DataItem",
-                                     {"Format":ar.dataItemFormat,
-                                      **xdmf_number_type(self.globalMesh.nodeNumbering_subdomain2global),
-                                      "Dimensions":"%i" % (self.nNodes_global,)})
-                elemMapAtt = SubElement(self.arGrid,"Attribute",
-                                        {"Name":"CellMapL2G",
-                                         "AttributeType":"Scalar",
-                                         "Center":"Cell"})
-                elemMap = SubElement(elemMapAtt,"DataItem",
-                                     {"Format":ar.dataItemFormat,
-                                      **xdmf_number_type(self.globalMesh.elementNumbering_subdomain2global),
-                                      "Dimensions":"%i" % (self.nElements_owned,)})
-
-                if ar.hdfFile is not None:
-                    nodeMap.text = ar.hdfFilename+":/nodeMapL2G"+str(ar.comm.rank())+name+str(tCount)
-                    elemMap.text = ar.hdfFilename+":/cellMapL2G"+str(ar.comm.rank())+name+str(tCount)
-                    if init or meshChanged:
-                        ar.create_dataset_async('nodeMapL2G'+str(ar.comm.rank())+name+str(tCount), data=self.globalMesh.nodeNumbering_subdomain2global)
-                        ar.create_dataset_async('cellMapL2G'+str(ar.comm.rank())+name+str(tCount), data=self.globalMesh.elementNumbering_subdomain2global[:self.nElements_owned])
-                else:
-                    SubElement(nodeMap,"xi:include",{"parse":"text","href":"./"+ar.textDataDir+"/nodeMapL2G"+name+".txt"})
-                    SubElement(nodeMap,"xi:include",{"parse":"text","href":"./"+ar.textDataDir+"/cellMapL2G"+name+".txt"})
-                    if init or meshChanged:
-                        np.savetxt(ar.textDataDir+"/nodeMapL2G"+name+".txt",self.globalMesh.nodeNumbering_subdomain2global)
-                        np.savetxt(ar.textDataDir+"/cellMapL2G"+name+".txt",self.globalMesh.elementNumbering_subdomain2global[:self.nElements_owned])
+                ar.write_field(self.arGrid, "NodeMapL2G",
+                               self.globalMesh.nodeNumbering_subdomain2global, tCount,
+                               dimensions=[self.nNodes_global],
+                               dataset='nodeMapL2G'+str(rank)+name+str(tCount),
+                               text_stem='nodeMapL2G'+name,
+                               write_data=init or meshChanged)
+                ar.write_field(self.arGrid, "CellMapL2G",
+                               self.globalMesh.elementNumbering_subdomain2global[:self.nElements_owned], tCount,
+                               center="Cell", dimensions=[self.nElements_owned],
+                               dataset='cellMapL2G'+str(rank)+name+str(tCount),
+                               text_stem='cellMapL2G'+name,
+                               write_data=init or meshChanged)
             #
             #material types
             #
             if ar.global_sync:
-                nodeMaterialTypes = SubElement(self.arGrid,"Attribute",{"Name":"nodeMaterialTypes",
-                                                                        "AttributeType":"Scalar",
-                                                                        "Center":"Node"})
-                nodeMaterialTypesValues = SubElement(nodeMaterialTypes,"DataItem",
-                                                     {"Format":ar.dataItemFormat,
-                                                      **xdmf_number_type(self.nodeMaterialTypes),
-                                                      "Dimensions":"%i" % (self.globalMesh.nNodes_global,)})
-                elementMaterialTypes = SubElement(self.arGrid,"Attribute",{"Name":"elementMaterialTypes",
-                                                                           "AttributeType":"Scalar",
-                                                                           "Center":"Cell"})
-                elementMaterialTypesValues = SubElement(elementMaterialTypes,"DataItem",
-                                                        {"Format":ar.dataItemFormat,
-                                                         **xdmf_number_type(self.elementMaterialTypes),
-                                                         "Dimensions":"%i" % (self.globalMesh.nElements_global,)})
+                ar.write_field(self.arGrid, "nodeMaterialTypes", self.nodeMaterialTypes, tCount,
+                               dimensions=[self.globalMesh.nNodes_global],
+                               sync_offsets=self.globalMesh.nodeOffsets_subdomain_owned,
+                               sync_data=self.nodeMaterialTypes[:self.nNodes_owned])
+                ar.write_field(self.arGrid, "elementMaterialTypes", self.elementMaterialTypes, tCount,
+                               center="Cell", dimensions=[self.globalMesh.nElements_global],
+                               sync_offsets=self.globalMesh.elementOffsets_subdomain_owned,
+                               sync_data=self.elementMaterialTypes[:self.nElements_owned])
                 if EB:
-                    ebnodeMaterialTypes = SubElement(self.arEBGrid,"Attribute",{"Name":"ebnodeMaterialTypes",
-                                                                          "AttributeType":"Scalar",
-                                                                          "Center":"Node"})
-                    ebnodeMaterialTypesValues = SubElement(ebnodeMaterialTypes,"DataItem",
-                                                       {"Format":ar.dataItemFormat,
-                                                        **xdmf_number_type(self.nodeMaterialTypes),
-                                                        "Dimensions":"%i" % (self.globalMesh.nNodes_global,)})
-                    elementBoundaryMaterialTypes = SubElement(self.arEBGrid,"Attribute",{"Name":"elementBoundaryMaterialTypes",
-                                                                                             "AttributeType":"Scalar",
-                                                                                             "Center":"Cell"})
-                    elementBoundaryMaterialTypesValues = SubElement(elementBoundaryMaterialTypes,"DataItem",
-                                                          {"Format":ar.dataItemFormat,
-                                                           **xdmf_number_type(self.elementBoundaryMaterialTypes),
-                                                           "Dimensions":"%i" % (self.globalMesh.nElementBoundaries_global,)})
-                if ar.hdfFile is not None:
-                    nodeMaterialTypesValues.text = ar.hdfFilename+":/"+"nodeMaterialTypes"+"_t"+str(tCount)
-                    ar.create_dataset_sync("nodeMaterialTypes"+"_t"+str(tCount), offsets=self.globalMesh.nodeOffsets_subdomain_owned, data=self.nodeMaterialTypes[:self.nNodes_owned])
-                    elementMaterialTypesValues.text = ar.hdfFilename+":/"+"elementMaterialTypes"+"_t"+str(tCount)
-                    ar.create_dataset_sync("elementMaterialTypes"+"_t"+str(tCount), offsets=self.globalMesh.elementOffsets_subdomain_owned, data=self.elementMaterialTypes[:self.nElements_owned])
-                    if EB:
-                        ebnodeMaterialTypesValues.text = ar.hdfFilename+":/"+"nodeMaterialTypes"+"_t"+str(tCount)
-                        elementBoundaryMaterialTypesValues.text = ar.hdfFilename+":/"+"elementBoundaryMaterialTypes"+"_t"+str(tCount)
-                        ar.create_dataset_sync("elementBoundaryMaterialTypes"+"_t"+str(tCount), offsets = self.globalMesh.elementBoundaryOffsets_subdomain_owned, data=self.elementBoundaryMaterialTypes[:self.nElementBoundaries_owned])
-                else:
-                    assert False, "global_sync  not  supported with text heavy data"
+                    #the boundary grid's node types are the mesh's, just written
+                    ar.write_field(self.arEBGrid, "ebnodeMaterialTypes", self.nodeMaterialTypes, tCount,
+                                   dimensions=[self.globalMesh.nNodes_global],
+                                   dataset="nodeMaterialTypes_t"+str(tCount), write_data=False)
+                    ar.write_field(self.arEBGrid, "elementBoundaryMaterialTypes",
+                                   self.elementBoundaryMaterialTypes, tCount,
+                                   center="Cell", dimensions=[self.globalMesh.nElementBoundaries_global],
+                                   sync_offsets=self.globalMesh.elementBoundaryOffsets_subdomain_owned,
+                                   sync_data=self.elementBoundaryMaterialTypes[:self.nElementBoundaries_owned])
             else:
-                nodeMaterialTypes = SubElement(self.arGrid,"Attribute",{"Name":"nodeMaterialTypes",
-                                                                        "AttributeType":"Scalar",
-                                                                        "Center":"Node"})
-                nodeMaterialTypesValues = SubElement(nodeMaterialTypes,"DataItem",
-                                                     {"Format":ar.dataItemFormat,
-                                                      **xdmf_number_type(self.nodeMaterialTypes),
-                                                      "Dimensions":"%i" % (self.nNodes_global,)})
-                elementMaterialTypes = SubElement(self.arGrid,"Attribute",{"Name":"elementMaterialTypes",
-                                                                           "AttributeType":"Scalar",
-                                                                           "Center":"Cell"})
-                elementMaterialTypesValues = SubElement(elementMaterialTypes,"DataItem",
-                                                        {"Format":ar.dataItemFormat,
-                                                         **xdmf_number_type(self.elementMaterialTypes),
-                                                         "Dimensions":"%i" % (self.nElements_owned,)})
+                ar.write_field(self.arGrid, "nodeMaterialTypes", self.nodeMaterialTypes, tCount,
+                               dimensions=[self.nNodes_global])
+                ar.write_field(self.arGrid, "elementMaterialTypes",
+                               self.elementMaterialTypes[:self.nElements_owned], tCount,
+                               center="Cell", dimensions=[self.nElements_owned])
                 if EB:
-                    ebnodeMaterialTypes = SubElement(self.arEBGrid,"Attribute",{"Name":"ebnodeMaterialTypes",
-                                                                          "AttributeType":"Scalar",
-                                                                          "Center":"Node"})
-                    ebnodeMaterialTypesValues = SubElement(ebnodeMaterialTypes,"DataItem",
-                                                       {"Format":ar.dataItemFormat,
-                                                        **xdmf_number_type(self.nodeMaterialTypes),
-                                                        "Dimensions":"%i" % (self.nNodes_global,)})
-                    elementBoundaryMaterialTypes = SubElement(self.arEBGrid,"Attribute",{"Name":"elementBoundaryMaterialTypes",
-                                                                                             "AttributeType":"Scalar",
-                                                                                             "Center":"Cell"})
-                    elementBoundaryMaterialTypesValues = SubElement(elementBoundaryMaterialTypes,"DataItem",
-                                                          {"Format":ar.dataItemFormat,
-                                                           **xdmf_number_type(self.elementBoundaryMaterialTypes),
-                                                           "Dimensions":"%i" % (self.nElementBoundaries_global,)})
-                if ar.hdfFile is not None:
-                    nodeMaterialTypesValues.text = ar.hdfFilename+":/"+"nodeMaterialTypes"+"_p"+str(ar.comm.rank())+"_t"+str(tCount)
-                    ar.create_dataset_async("nodeMaterialTypes"+"_p"+str(ar.comm.rank())+"_t"+str(tCount), data=self.nodeMaterialTypes)
-                    elementMaterialTypesValues.text = ar.hdfFilename+":/"+"elementMaterialTypes"+"_p"+str(ar.comm.rank())+"_t"+str(tCount)
-                    ar.create_dataset_async("elementMaterialTypes"+"_p"+str(ar.comm.rank())+"_t"+str(tCount), data=self.elementMaterialTypes[:self.nElements_owned])
-                    if EB:
-                        ebnodeMaterialTypesValues.text = ar.hdfFilename+":/"+"nodeMaterialTypes"+"_p"+str(ar.comm.rank())+"_t"+str(tCount)
-                        elementBoundaryMaterialTypesValues.text = ar.hdfFilename+":/"+"elementBoundaryMaterialTypes"+"_p"+str(ar.comm.rank())+"_t"+str(tCount)
-                        ar.create_dataset_async("elementBoundaryMaterialTypes"+"_p"+str(ar.comm.rank())+"_t"+str(tCount), data=self.elementBoundaryMaterialTypes)
-                else:
-                    numpy.savetxt(ar.textDataDir+"/"+"nodeMaterialTypes"+str(tCount)+".txt",self.nodeMaterialTypes)
-                    SubElement(nodeMaterialTypesValues,"xi:include",{"parse":"text","href":"./"+ar.textDataDir+"/"+"nodeMaterialTypes"+str(tCount)+".txt"})
-                    numpy.savetxt(ar.textDataDir+"/"+"elementMaterialTypes"+str(tCount)+".txt",self.elementMaterialTypes[:self.nElements_owned])
-                    SubElement(elementMaterialTypesValues,"xi:include",{"parse":"text","href":"./"+ar.textDataDir+"/"+"elementMaterialTypes"+str(tCount)+".txt"})
+                    ar.write_field(self.arEBGrid, "ebnodeMaterialTypes", self.nodeMaterialTypes, tCount,
+                                   dimensions=[self.nNodes_global],
+                                   dataset="nodeMaterialTypes_p"+str(rank)+"_t"+str(tCount),
+                                   write_data=False)
+                    ar.write_field(self.arEBGrid, "elementBoundaryMaterialTypes",
+                                   self.elementBoundaryMaterialTypes, tCount,
+                                   center="Cell", dimensions=[self.nElementBoundaries_global])
             #done with material types
     def buildFromC(self,cmesh):
         from . import cmeshTools
@@ -4496,94 +4399,59 @@ def writeMeshXdmf(self,ar,name='',t=0.0,init=False,meshChanged=False,Xdmf_Elemen
         if self.arGridCollection is not None:
             init = False
         elif not init:
-            self.arGridCollection = ar.domain.find("Grid")
+            self.arGridCollection = ar.domain[0]
         if init:
-            self.arGridCollection = SubElement(ar.domain,"Grid",{"Name":"Mesh "+name,
-                                                                "GridType":"Collection",
-                                                                "CollectionType":"Temporal"})
-        if self.arGrid is None or self.arTime.get('Value') != "{0:e}".format(t):
+            self.arGridCollection = ar.domain.collection("Mesh "+name)
+        if not self.arGridCollection.is_current(t):
+            rank = ar.comm.rank()
             #
             #topology and geometry
             #
             if ar.global_sync:
-                self.arGrid = SubElement(self.arGridCollection,"Grid",{"GridType":"Uniform"})
-                self.arTime = SubElement(self.arGrid,"Time",{"Value":str(t),"Name":str(tCount)})
-                topology = SubElement(self.arGrid,"Topology",
-                                      {"Type":Xdmf_ElementTopology,
-                                       "NumberOfElements":"%i" % (self.globalMesh.nElements_global,)})
-                elements = SubElement(topology,"DataItem",
-                                      {"Format":ar.dataItemFormat,
-                                       **xdmf_number_type(self.globalMesh.nodeNumbering_subdomain2global),
-                                       "Dimensions":"%i %i" % (self.globalMesh.nElements_global,
-                                                               self.nNodes_element)})
-                geometry = SubElement(self.arGrid,"Geometry",{"Type":"XYZ"})
-                nodes    = SubElement(geometry,"DataItem",
-                                      {"Format":ar.dataItemFormat,
-                                       **xdmf_number_type(self.nodeArray),
-                                       "Dimensions":"%i %i" % (self.globalMesh.nNodes_global,3)})
-                #material types
-                elementMaterialTypes = SubElement(self.arGrid,"Attribute",{"Name":"elementMaterialTypes",
-                                                                           "AttributeType":"Scalar",
-                                                                           "Center":"Cell"})
-                elementMaterialTypesValues = SubElement(elementMaterialTypes,"DataItem",
-                                                        {"Format":ar.dataItemFormat,
-                                                         **xdmf_number_type(self.elementMaterialTypes),
-                                                         "Dimensions":"%i" % (self.globalMesh.nElements_global,)})
-                if ar.hdfFile is not None:
-                    elements.text = ar.hdfFilename+":/elements"+name+str(tCount)
-                    nodes.text = ar.hdfFilename+":/nodes"+name+str(tCount)
-                    elementMaterialTypesValues.text = ar.hdfFilename+":/"+"elementMaterialTypes"+"_t"+str(tCount)
-                    if init or meshChanged:
-                        ar.create_dataset_sync('elements'+name+str(tCount),
-                                               offsets = self.globalMesh.elementOffsets_subdomain_owned,
-                                               data = self.globalMesh.nodeNumbering_subdomain2global[self.elementNodesArray[:self.nElements_owned]])
-                        ar.create_dataset_sync('nodes'+name+str(tCount),
-                                               offsets = self.globalMesh.nodeOffsets_subdomain_owned,
-                                               data = self.nodeArray[:self.nNodes_owned])
-                        ar.create_dataset_sync("elementMaterialTypes"+"_t"+str(tCount),
-                                               offsets = self.globalMesh.elementOffsets_subdomain_owned,
-                                               data = self.elementMaterialTypes[:self.nElements_owned])
-                else:
-                    assert False, "global_sync with text heavy data not supported"
+                assert ar.hdfFile is not None, "global_sync with text heavy data not supported"
+                n_elements       = self.globalMesh.nElements_global
+                nodes_dims       = [self.globalMesh.nNodes_global, 3]
+                elements_dataset = 'elements'+name+str(tCount)
+                nodes_dataset    = 'nodes'+name+str(tCount)
+                elements_dtype   = self.globalMesh.nodeNumbering_subdomain2global
             else:
-                self.arGrid = SubElement(self.arGridCollection,"Grid",{"GridType":"Uniform"})
-                self.arTime = SubElement(self.arGrid,"Time",{"Value":str(t),"Name":str(tCount)})
-                topology = SubElement(self.arGrid,"Topology",
-                                      {"Type":Xdmf_ElementTopology,
-                                       "NumberOfElements":"%i" % (self.nElements_owned,)})
-                elements = SubElement(topology,"DataItem",
-                                      {"Format":ar.dataItemFormat,
-                                       **xdmf_number_type(self.elementNodesArray),
-                                       "Dimensions":"%i %i" % (self.nElements_owned,self.nNodes_element)})
-                geometry = SubElement(self.arGrid,"Geometry",{"Type":"XYZ"})
-                nodes    = SubElement(geometry,"DataItem",
-                                      {"Format":ar.dataItemFormat,
-                                       **xdmf_number_type(self.nodeArray),
-                                       "Dimensions":"%i %i" % (self.nNodes_global,3)})
-                #material types
-                elementMaterialTypes = SubElement(self.arGrid,"Attribute",{"Name":"elementMaterialTypes",
-                                                                           "AttributeType":"Scalar",
-                                                                           "Center":"Cell"})
-                elementMaterialTypesValues = SubElement(elementMaterialTypes,"DataItem",
-                                                        {"Format":ar.dataItemFormat,
-                                                         **xdmf_number_type(self.elementMaterialTypes),
-                                                         "Dimensions":"%i" % (self.nElements_owned,)})
-                if ar.hdfFile is not None:
-                    elements.text = ar.hdfFilename+":/elements"+str(ar.comm.rank())+name+str(tCount)
-                    nodes.text = ar.hdfFilename+":/nodes"+str(ar.comm.rank())+name+str(tCount)
-                    elementMaterialTypesValues.text = ar.hdfFilename+":/"+"elementMaterialTypes"+"_p"+str(ar.comm.rank())+"_t"+str(tCount)
-                    if init or meshChanged:
-                        ar.create_dataset_async('elements'+str(ar.comm.rank())+name+str(tCount), data = self.elementNodesArray[:self.nElements_owned])
-                        ar.create_dataset_async('nodes'+str(ar.comm.rank())+name+str(tCount), data = self.nodeArray)
-                        ar.create_dataset_async("elementMaterialTypes"+"_p"+str(ar.comm.rank())+"_t"+str(tCount), data = self.elementMaterialTypes[:self.nElements_owned])
+                n_elements       = self.nElements_owned
+                nodes_dims       = [self.nNodes_global, 3]
+                elements_dataset = 'elements'+str(rank)+name+str(tCount)
+                nodes_dataset    = 'nodes'+str(rank)+name+str(tCount)
+                elements_dtype   = self.elementNodesArray
+            self.arGrid = ar.write_grid(self.arGridCollection, None, t, tCount)
+            ar.write_topology(self.arGrid, Xdmf_ElementTopology, n_elements,
+                              [n_elements, self.nNodes_element],
+                              elements_dataset, 'elements'+name, dtype=elements_dtype)
+            ar.write_geometry(self.arGrid, nodes_dims, nodes_dataset, 'nodes'+name,
+                              dtype=self.nodeArray)
+            #material types
+            if ar.global_sync:
+                ar.write_field(self.arGrid, "elementMaterialTypes", self.elementMaterialTypes, tCount,
+                               center="Cell", dimensions=[self.globalMesh.nElements_global],
+                               sync_offsets=self.globalMesh.elementOffsets_subdomain_owned,
+                               sync_data=self.elementMaterialTypes[:self.nElements_owned],
+                               write_data=init or meshChanged)
+            else:
+                ar.write_field(self.arGrid, "elementMaterialTypes",
+                               self.elementMaterialTypes[:self.nElements_owned], tCount,
+                               center="Cell", dimensions=[self.nElements_owned],
+                               write_data=init or meshChanged)
+            if init or meshChanged:
+                if ar.global_sync:
+                    ar.create_dataset_sync(elements_dataset,
+                                           offsets = self.globalMesh.elementOffsets_subdomain_owned,
+                                           data = self.globalMesh.nodeNumbering_subdomain2global[self.elementNodesArray[:self.nElements_owned]])
+                    ar.create_dataset_sync(nodes_dataset,
+                                           offsets = self.globalMesh.nodeOffsets_subdomain_owned,
+                                           data = self.nodeArray[:self.nNodes_owned])
+                elif ar.hdfFile is not None:
+                    ar.create_dataset_async(elements_dataset, data = self.elementNodesArray[:self.nElements_owned])
+                    ar.create_dataset_async(nodes_dataset, data = self.nodeArray)
                 else:
-                    SubElement(elements,"xi:include",{"parse":"text","href":"./"+ar.textDataDir+"/elements"+name+".txt"})
-                    SubElement(nodes,"xi:include",{"parse":"text","href":"./"+ar.textDataDir+"/nodes"+name+".txt"})
-                    SubElement(elementMaterialTypesValues,"xi:include",{"parse":"text","href":"./"+ar.textDataDir+"/"+"elementMaterialTypes"+str(tCount)+".txt"})
-                    if init or meshChanged:
-                        np.savetxt(ar.textDataDir+"/elements"+name+".txt",self.elementNodesArray[:self.nElements_owned],fmt='%d')
-                        np.savetxt(ar.textDataDir+"/nodes"+name+".txt",self.nodeArray)
-                        np.savetxt(ar.textDataDir+"/"+"elementMaterialTypes"+str(tCount)+".txt",self.elementMaterialTypes[:self.nElements_owned])
+                    np.savetxt(ar.textDataDir+"/elements"+name+".txt",self.elementNodesArray[:self.nElements_owned],fmt='%d')
+                    np.savetxt(ar.textDataDir+"/nodes"+name+".txt",self.nodeArray)
 
 
 class Mesh3DM(Mesh):
@@ -4849,97 +4717,60 @@ def writeMeshXdmf(self,ar,name='',t=0.0,init=False,meshChanged=False,Xdmf_Elemen
         if self.arGridCollection is not None:
             init = False
         elif not init:
-            self.arGridCollection = ar.domain.find("Grid")
+            self.arGridCollection = ar.domain[0]
         if init:
-            self.arGridCollection = SubElement(ar.domain,"Grid",{"Name":"Mesh "+name,
-                                                               "GridType":"Collection",
-                                                               "CollectionType":"Temporal"})
-        if self.arGrid is None or self.arTime.get('Value') != "{0:e}".format(t):
+            self.arGridCollection = ar.domain.collection("Mesh "+name)
+        if not self.arGridCollection.is_current(t):
+            rank = ar.comm.rank()
+            #
+            #topology and geometry
+            #
             if ar.global_sync:
-                #
-                #topology and geometry
-                #
-                self.arGrid = SubElement(self.arGridCollection,"Grid",{"GridType":"Uniform"})
-                self.arTime = SubElement(self.arGrid,"Time",{"Value":str(t),"Name":str(tCount)})
-                topology = SubElement(self.arGrid,"Topology",
-                                      {"Type":Xdmf_ElementTopology,
-                                       "NumberOfElements":"%i" % (self.globalMesh.nElements_global,)})
-                elements = SubElement(topology,"DataItem",
-                                      {"Format":ar.dataItemFormat,
-                                       **xdmf_number_type(self.globalMesh.nodeNumbering_subdomain2global),
-                                       "Dimensions":"%i %i" % (self.globalMesh.nElements_global,
-                                                               self.nNodes_element)})
-                geometry = SubElement(self.arGrid,"Geometry",{"Type":"XYZ"})
-                nodes    = SubElement(geometry,"DataItem",
-                                      {"Format":ar.dataItemFormat,
-                                       **xdmf_number_type(self.nodeArray),
-                                       "Dimensions":"%i %i" % (self.globalMesh.nNodes_global,3)})
-                #material types
-                elementMaterialTypes = SubElement(self.arGrid,"Attribute",{"Name":"elementMaterialTypes",
-                                                                           "AttributeType":"Scalar",
-                                                                           "Center":"Cell"})
-                elementMaterialTypesValues = SubElement(elementMaterialTypes,"DataItem",
-                                                        {"Format":ar.dataItemFormat,
-                                                         **xdmf_number_type(self.elementMaterialTypes),
-                                                         "Dimensions":"%i" % (self.globalMesh.nElements_global,)})
-                if ar.hdfFile is not None:
-                    elements.text = ar.hdfFilename+":/elements"+name+str(tCount)
-                    nodes.text = ar.hdfFilename+":/nodes"+name+str(tCount)
-                    elementMaterialTypesValues.text = ar.hdfFilename+":/"+"elementMaterialTypes"+"_t"+str(tCount)
-                    if init or meshChanged:
-                        ar.create_dataset_sync('elements'+name+str(tCount),
-                                               offsets = self.globalMesh.elementOffsets_subdomain_owned,
-                                               data = self.globalMesh.nodeNumbering_subdomain2global[self.elementNodesArray[:self.nElements_owned]])
-                        ar.create_dataset_sync('nodes'+name+str(tCount),
-                                               offsets = self.globalMesh.nodeOffsets_subdomain_owned,
-                                               data = self.nodeArray[:self.nNodes_owned])
-                        ar.create_dataset_sync("elementMaterialTypes"+"_t"+str(tCount),
-                                               offsets = self.globalMesh.elementOffsets_subdomain_owned,
-                                               data = self.elementMaterialTypes[:self.nElements_owned])
-                else:
-                    assert False, "global_sync not supported  with text heavy data"
+                assert ar.hdfFile is not None, "global_sync with text heavy data not supported"
+                n_elements       = self.globalMesh.nElements_global
+                nodes_dims       = [self.globalMesh.nNodes_global, 3]
+                elements_dataset = 'elements'+name+str(tCount)
+                nodes_dataset    = 'nodes'+name+str(tCount)
+                elements_dtype   = self.globalMesh.nodeNumbering_subdomain2global
             else:
-                #
-                #topology and geometry
-                #
-                self.arGrid = SubElement(self.arGridCollection,"Grid",{"GridType":"Uniform"})
-                self.arTime = SubElement(self.arGrid,"Time",{"Value":str(t),"Name":str(tCount)})
-                topology = SubElement(self.arGrid,"Topology",
-                                      {"Type":Xdmf_ElementTopology,
-                                       "NumberOfElements":"%i" % (self.nElements_owned,)})
-                elements = SubElement(topology,"DataItem",
-                                      {"Format":ar.dataItemFormat,
-                                       **xdmf_number_type(self.elementNodesArray),
-                                       "Dimensions":"%i %i" % (self.nElements_owned,self.nNodes_element)})
-                geometry = SubElement(self.arGrid,"Geometry",{"Type":"XYZ"})
-                nodes    = SubElement(geometry,"DataItem",
-                                      {"Format":ar.dataItemFormat,
-                                       **xdmf_number_type(self.nodeArray),
-                                       "Dimensions":"%i %i" % (self.nNodes_global,3)})
-                #material types
-                elementMaterialTypes = SubElement(self.arGrid,"Attribute",{"Name":"elementMaterialTypes",
-                                                                           "AttributeType":"Scalar",
-                                                                           "Center":"Cell"})
-                elementMaterialTypesValues = SubElement(elementMaterialTypes,"DataItem",
-                                                        {"Format":ar.dataItemFormat,
-                                                         **xdmf_number_type(self.elementMaterialTypes),
-                                                         "Dimensions":"%i" % (self.nElements_owned,)})
-                if ar.hdfFile is not None:
-                    elements.text = ar.hdfFilename+":/elements"+str(ar.comm.rank())+name+str(tCount)
-                    nodes.text = ar.hdfFilename+":/nodes"+str(ar.comm.rank())+name+str(tCount)
-                    elementMaterialTypesValues.text = ar.hdfFilename+":/"+"elementMaterialTypes"+"_p"+str(ar.comm.rank())+"_t"+str(tCount)
-                    if init or meshChanged:
-                        ar.create_dataset_async('elements'+str(ar.comm.rank())+name+str(tCount), data = self.elementNodesArray[:self.nElements_owned])
-                        ar.create_dataset_async('nodes'+str(ar.comm.rank())+name+str(tCount), data = self.nodeArray)
-                        ar.create_dataset_async("elementMaterialTypes"+"_p"+str(ar.comm.rank())+"_t"+str(tCount), data = self.elementMaterialTypes[:self.nElements_owned])
+                n_elements       = self.nElements_owned
+                nodes_dims       = [self.nNodes_global, 3]
+                elements_dataset = 'elements'+str(rank)+name+str(tCount)
+                nodes_dataset    = 'nodes'+str(rank)+name+str(tCount)
+                elements_dtype   = self.elementNodesArray
+            self.arGrid = ar.write_grid(self.arGridCollection, None, t, tCount)
+            ar.write_topology(self.arGrid, Xdmf_ElementTopology, n_elements,
+                              [n_elements, self.nNodes_element],
+                              elements_dataset, 'elements'+name, dtype=elements_dtype)
+            ar.write_geometry(self.arGrid, nodes_dims, nodes_dataset, 'nodes'+name,
+                              dtype=self.nodeArray)
+            #material types
+            if ar.global_sync:
+                ar.write_field(self.arGrid, "elementMaterialTypes", self.elementMaterialTypes, tCount,
+                               center="Cell", dimensions=[self.globalMesh.nElements_global],
+                               sync_offsets=self.globalMesh.elementOffsets_subdomain_owned,
+                               sync_data=self.elementMaterialTypes[:self.nElements_owned],
+                               write_data=init or meshChanged)
+            else:
+                ar.write_field(self.arGrid, "elementMaterialTypes",
+                               self.elementMaterialTypes[:self.nElements_owned], tCount,
+                               center="Cell", dimensions=[self.nElements_owned],
+                               write_data=init or meshChanged)
+            if init or meshChanged:
+                if ar.global_sync:
+                    ar.create_dataset_sync(elements_dataset,
+                                           offsets = self.globalMesh.elementOffsets_subdomain_owned,
+                                           data = self.globalMesh.nodeNumbering_subdomain2global[self.elementNodesArray[:self.nElements_owned]])
+                    ar.create_dataset_sync(nodes_dataset,
+                                           offsets = self.globalMesh.nodeOffsets_subdomain_owned,
+                                           data = self.nodeArray[:self.nNodes_owned])
+                elif ar.hdfFile is not None:
+                    ar.create_dataset_async(elements_dataset, data = self.elementNodesArray[:self.nElements_owned])
+                    ar.create_dataset_async(nodes_dataset, data = self.nodeArray)
                 else:
-                    SubElement(elements,"xi:include",{"parse":"text","href":"./"+ar.textDataDir+"/elements"+name+".txt"})
-                    SubElement(nodes,"xi:include",{"parse":"text","href":"./"+ar.textDataDir+"/nodes"+name+".txt"})
-                    SubElement(elementMaterialTypesValues,"xi:include",{"parse":"text","href":"./"+ar.textDataDir+"/"+"elementMaterialTypes"+str(tCount)+".txt"})
-                    if init or meshChanged:
-                        np.savetxt(ar.textDataDir+"/elements"+name+".txt",self.elementNodesArray[:self.nElements_owned],fmt='%d')
-                        np.savetxt(ar.textDataDir+"/nodes"+name+".txt",self.nodeArray)
-                        np.savetxt(ar.textDataDir+"/"+"elementMaterialTypes"+str(tCount)+".txt",self.elementMaterialTypes[:self.nElements_owned])
+                    np.savetxt(ar.textDataDir+"/elements"+name+".txt",self.elementNodesArray[:self.nElements_owned],fmt='%d')
+                    np.savetxt(ar.textDataDir+"/nodes"+name+".txt",self.nodeArray)
+
 
 class MultilevelTetrahedralMesh(MultilevelMesh):
     """A hierarchical multilevel mesh with tetrahedral cells"""
@@ -7144,80 +6975,60 @@ def refine(self):
 
 ## @}
 
-###utility functions for reading meshes from Xdmf
+###utility functions for reading meshes back from an archive
 
 def _dataitem_dataset(node, label):
-    """The HDF5 dataset name a node's DataItem refers to.
+    """The HDF5 dataset name a topology, geometry or attribute's data item refers to.
 
-    Replaces ``node[0].text.split(':')[-1]``, which reached into the first
-    child by position and split the reference on any colon. A DataItem
-    reference is ``<file>:/<dataset>``, so splitting on ``:/`` keeps a
-    dataset path containing a colon intact.
+    A data item reference is ``<file>:/<dataset>``, so splitting on ``:/``
+    keeps a dataset path that itself contains a colon intact.
     """
-    item = node.find('DataItem')
-    assert item is not None, "no <DataItem> under <%s>" % (label,)
-    text = (item.text or "").strip()
-    assert text, "<DataItem> under <%s> has no reference" % (label,)
+    item = node.get("data_item") if node is not None else None
+    assert item is not None, "no data_item under %s" % (label,)
+    text = (item.get("data") or "").strip()
+    assert text, "the data_item under %s has no HDF5 reference" % (label,)
     return text.split(':/')[-1]
 
 
-def findXMLgridElement(xmf,MeshTag='Spatial_Domain',id_in_collection=-1,verbose=0):
-    """Try to find the element of the xml tree xmf that holds a uniform
-    grid with the name given in MeshTag by searching through Temporal
-    Grid Collections and Grid Collections.
+def findArchiveGrid(domain, MeshTag='Spatial_Domain', id_in_collection=-1, verbose=0):
+    """One grid of an archive's domain: a step of the collection named by ``MeshTag``.
 
-    If MeshTag isn't found, uses the first entry in the Domain
+    ``domain`` is the ymf domain dict (:func:`proteus.Archiver.readArchiveDomain`
+    reads one from an archive's ``.h5``). The collection is the first whose
+    name contains ``MeshTag``, or the first of all when none does;
+    ``id_in_collection`` selects the step, -1 being the last written. A
+    per-subdomain step holds one grid per rank, of which the first is
+    returned, which is what a serial reader can use.
     """
-    #Find the Domain by tag rather than as the root's last child: an
-    #<Information> element is legal ahead of it, and a positional read
-    #silently picks up whatever happens to sit last.
-    Domain = xmf.getroot().find('Domain')
-    assert Domain is not None, "no <Domain> in the XDMF document"
-
-    GridCollection = None
-    for collection in Domain.findall('Grid'):
-        if MeshTag in collection.attrib.get('Name', ''):
-            GridCollection = collection
+    collections = domain.get("time_collections", [])
+    assert collections, "the archive has no time collections"
+    collection = None
+    for c in collections:
+        if MeshTag in c.get("name", ""):
+            collection = c
             break
-    if GridCollection is None:
-        candidates = Domain.findall('Grid')
-        assert candidates, "no <Grid> under <Domain>"
-        GridCollection = candidates[0]
+    if collection is None:
+        collection = collections[0]
         logEvent("No grid collection named %r; using %r"
-                 % (MeshTag, GridCollection.attrib.get('Name', '<unnamed>')), 3)
-    logEvent("Trying GridCollection.tag= %s" % (GridCollection.tag),4)
-
-    if GridCollection.attrib['GridType'] == 'Collection':
-        #id_in_collection selects the step; -1 is the last one written
-        steps = GridCollection.findall('Grid')
-        assert steps, ("grid collection %r holds no grids"
-                       % (GridCollection.attrib.get('Name', '<unnamed>'),))
-        Grid = steps[id_in_collection]
-        if Grid.attrib.get('GridType') == 'Collection':
-            #a spatial collection of subdomain grids: take the first, which
-            #is what a serial reader can use
-            subdomains = Grid.findall('Grid')
-            assert subdomains, "spatial collection holds no subdomain grids"
-            Grid = subdomains[0]
-    elif GridCollection.attrib['GridType'] == 'Uniform':
-        Grid = GridCollection
-    assert Grid.tag == 'Grid'
-    assert Grid.attrib['GridType'] == 'Uniform'
-
-    return Grid
-
-def extractPropertiesFromXdmfGridNode(Grid):
-    """unpack the Topology, Geometry, NodeMaterials, and ElementMaterials
-    nodes from xdmf node for a uniform grid
-    """
-    #Look these up by tag and name. The previous version iterated with
-    #enumerate and then re-indexed Grid[i] to get the very leaf it was
-    #already holding.
-    Topology = Grid.find('Topology')
-    Geometry = Grid.find('Geometry')
+                 % (MeshTag, collection.get("name", "<unnamed>")), 3)
+    steps = collection["data"]
+    assert steps, ("grid collection %r holds no grids"
+                   % (collection.get("name", "<unnamed>"),))
+    step = steps[id_in_collection]
+    if "spatial_collection" in step:
+        subdomains = step["spatial_collection"]
+        assert subdomains, "spatial collection holds no subdomain grids"
+        return subdomains[0]
+    return step
+
+
+def extractPropertiesFromGrid(grid):
+    """Unpack the topology, geometry, nodeMaterialTypes and elementMaterialTypes of a grid dict."""
+    Topology = grid.get("topology")
+    Geometry = grid.get("geometry")
     NodeMaterials = ElementMaterials = None
-    for attribute in Grid.findall('Attribute'):
-        name = attribute.attrib.get('Name')
+    for attribute in grid.get("attributes", []):
+        name = attribute.get("name")
         if name == 'nodeMaterialTypes':
             NodeMaterials = attribute
         elif name == 'elementMaterialTypes':
@@ -7231,7 +7042,7 @@ def extractPropertiesFromXdmfGridNode(Grid):
 
 def readUniformElementTopologyFromXdmf(elementTopologyName,Topology,hdf5,topologyid2name,topology2nodes):
     """
-    Read xmdf element topology information when there are uniform elements in the mesh
+    Read the element topology when there are uniform elements in the mesh
     Type of element given by elementTopologyName
     Heavy data stored in hdf5
     topologyid2name -- lookup for number of nodes in a given element type
@@ -7241,7 +7052,7 @@ def readUniformElementTopologyFromXdmf(elementTopologyName,Topology,hdf5,topolog
     nElements_global  -- the number of elements in the mesh
     nNodes_element    -- number of nodes per element
     elementNodesArray -- element --> node connectivity stored as flattened array accessed using elementNodes_offset
-    elementNodes_offset -- offsets into the elementNodesArray storage for element connectivity, 
+    elementNodes_offset -- offsets into the elementNodesArray storage for element connectivity,
     element eN nodes are in elementNodesArray[elementNodes_offset[eN]:elementNodes_offset[eN+1]]
 
     """
@@ -7261,19 +7072,19 @@ def readUniformElementTopologyFromXdmf(elementTopologyName,Topology,hdf5,topolog
 
 def readMixedElementTopologyFromXdmf(elementTopologyName,Topology,hdf5,topologyid2name,topology2nodes):
     """
-    Read xmdf element topology information when there are mixed elements in the mesh
+    Read the element topology when there are mixed elements in the mesh
     Heavy data stored in hdf5
     topologyid2name -- lookup for number of nodes in a given element type
 
     returns
 
     nElements_global  -- the number of elements in the mesh
-    elementNodesArray -- element --> node connectivity stored as flattened 
+    elementNodesArray -- element --> node connectivity stored as flattened
     array accessed using elementNodes_offset
-    elementNodes_offset -- offsets into the elementNodesArray storage for element 
-    connectivity, element eN nodes are 
+    elementNodes_offset -- offsets into the elementNodesArray storage for element
+    connectivity, element eN nodes are
     inelementNodesArray[elementNodes_offset[eN]:elementNodes_offset[eN+1]]
- 
+
     """
     assert elementTopologyName == 'Mixed'
 
@@ -7299,41 +7110,57 @@ def readMixedElementTopologyFromXdmf(elementTopologyName,Topology,hdf5,topologyi
         eN += 1; i += nNodes_local+1
     elementNodesArray = np.zeros((elementNodes_offset[nElements_global],),'i')
     i = 0; eN = 0
-    while i < len(self.xdmf_topology):
+    #this read self.xdmf_topology, which does not exist in a module function
+    while i < len(xdmf_topology):
         nNodes_local = topology2nodes[topologyid2name[xdmf_topology[i]]]
         elementNodesArray[elementNodes_offset[eN]:elementNodes_offset[eN+1]][:] = xdmf_topology[i+1:i+1+nNodes_local][:]
         eN += 1; i += nNodes_local+1
 
     return nElements_global, elementNodesArray, elementNodes_offset
 
+def readMeshFromArchive(archive_base,MeshTag="Spatial_Domain",verbose=0):
+    """Read a mesh back from a proteus archive, from its hdf5 file alone
+
+    ``archive_base`` is the archive's path without the ``.h5`` extension.
+    Its metadata (version 1, the XDMF fragments proteus <= 1.9.x stored, or
+    version 2) describes every step; the last step of the collection named
+    by ``MeshTag`` is read, so no ``.xmf`` is needed. Serial only: a
+    per-subdomain archive yields its first subdomain.
+
+    :return: a BasicMeshInfo object with the minimal information read
+
+    """
+    dataDir, filename = os.path.split(archive_base)
+    dataDir = dataDir or '.'
+    h5_path = os.path.join(dataDir, filename+'.h5')
+    assert os.path.isfile(h5_path)
+    domain = readArchiveDomain(filename, dataDir)
+    assert domain is not None and domain.get("time_collections"), \
+        ("%s holds no proteus archive metadata; a mesh in an XDMF file written "
+         "by another tool is read with readMeshXdmf" % (h5_path,))
+    return _readMeshFromDomain(domain, h5_path, MeshTag=MeshTag, verbose=verbose)
+
 def readMeshXdmf(xmf_archive_base,heavy_file_base,MeshTag="Spatial_Domain",hasHDF5=True,verbose=0):
-    """Read in a mesh from XDMF, assuming heavy data is in hdf5
+    """Read a mesh from an XDMF file, assuming heavy data is in hdf5
+
+    For XDMF written by other tools -- an OpenFOAM mesh converted to XDMF,
+    say -- whose ``.h5`` holds the arrays and nothing else. The XDMF is
+    parsed by ``ymf.xdmf``, so proteus itself imports no XML parser. A
+    proteus archive needs no ``.xmf``: :func:`readMeshFromArchive` reads it
+    from its ``.h5``.
 
     :return: a BasicMeshInfo object with the minimal information read
-    
+
     """
-    # start trying to read an xdmf archive with name xmf_archive_base.xmf
-    # assumes heavy_file_base.h5 has heavy data
-    # root Element is Xdmf
-    #   last child of Xdmf which should be a Domain Element
-    #      find child of Domain that is a Temporal Grid Collection with a name containing MeshTag, if None use first collection
-    #         last child of Temporal Grid Collection should be a Uniform Grid at final time
-    #            Attribute (usually 1) of child is  Topology
-    #               set elementTopologyName to Type
-    #               if Type != Mixed
-    #                 get text attribute and read this entry from  hdf5 file
-    #                 set nNodes_element based on Type, nElements_global from leading dimension of elementNodesArray
-    #                 create elementNodes_offset from Type and flatten elementNodesArray
-    #               else
-    #                 get text attribute and read this entry from  hdf5 file to place in into xdmf_topology
-    #                 generate elementNodesArray from xdmf_topology, calculating the number of elements using
-    #                   walk through xdmf_topology
-    #            Attribute (usually 2) of child is Geometry  --> load data into nodeArray
-    #                set nNodes_global from nodeArray
-    #            If has Attribute nodeMaterials read this from hdf file, else set to default of all zeros
-    #            If has Attribute elementMaterialTypes, read this from hdf file, else set to default of all zeros
+    from ymf.xdmf import read_xdmf
+    assert hasHDF5, "a mesh is read with its heavy data in an hdf5 file"
     assert os.path.isfile(xmf_archive_base+'.xmf')
     assert os.path.isfile(heavy_file_base+'.h5')
+    domain, _ = read_xdmf(xmf_archive_base+'.xmf')
+    return _readMeshFromDomain(domain, heavy_file_base+'.h5', MeshTag=MeshTag, verbose=verbose)
+
+def _readMeshFromDomain(domain, h5_path, MeshTag="Spatial_Domain", verbose=0):
+    """The mesh described by one grid of ``domain``, its arrays read from ``h5_path``."""
 
     ###information about allowed Xdmf topologies
     #Xdmf cell type id to Name
@@ -7360,13 +7187,11 @@ def __init__(self):
     #
     MeshInfo = BasicMeshInfo()
 
-    xmf = ET.parse(xmf_archive_base+'.xmf')
-    hdf5= h5py.File(heavy_file_base+'.h5',"r")
-    assert hasHDF5
+    hdf5= h5py.File(h5_path,"r")
 
-    Grid = findXMLgridElement(xmf,MeshTag,id_in_collection=-1,verbose=verbose)
+    Grid = findArchiveGrid(domain,MeshTag,id_in_collection=-1,verbose=verbose)
 
-    Topology,Geometry,NodeMaterials,ElementMaterials = extractPropertiesFromXdmfGridNode(Grid)
+    Topology,Geometry,NodeMaterials,ElementMaterials = extractPropertiesFromGrid(Grid)
 
     assert Geometry is not None
     entry = _dataitem_dataset(Geometry, "Geometry")
@@ -7383,10 +7208,7 @@ def __init__(self):
         MeshInfo.nodeMaterialTypes = np.zeros((MeshInfo.nNodes_global,),'i')
 
     assert Topology is not None
-    if 'Type' in Topology.attrib:
-        MeshInfo.elementTopologyName = Topology.attrib['Type']
-    elif 'TopologyType' in Topology.attrib:
-        MeshInfo.elementTopologyName = Topology.attrib['TopologyType']
+    MeshInfo.elementTopologyName = Topology.get("type")
     assert MeshInfo.elementTopologyName is not None
 
     logEvent("elementTopologyName= %s " % MeshInfo.elementTopologyName,3)
@@ -7417,6 +7239,7 @@ def __init__(self):
 
     return MeshInfo
 #
+
 def writeHexMesh(mesh_info,hexfile_base,index_base=0):
     """
     Write a hex mesh in Ido's format with base numbering index_base
diff --git a/proteus/NumericalSolution.py b/proteus/NumericalSolution.py
index 38ce98aa85..93309b9ad9 100644
--- a/proteus/NumericalSolution.py
+++ b/proteus/NumericalSolution.py
@@ -608,8 +608,8 @@ def calculateSolution(self,runName):
             if not self.opts.hotStart:
                 logEvent("Archiving initial conditions")
                 self.archiveInitialSolution(m,index)
-            else:
-                self.ar[index].domain = self.ar[index].tree.find("Domain")
+            #a hot-started archive seeded its domain with the collections
+            #the archive already has when it was opened
             #if(not hasattr(self.pList[0].domain,'PUMIMesh') and not self.opts.hotStart):
             self.initializeViewSolution(m)
             logEvent("Estimating initial time derivative and initializing time history for model "+p.name)
@@ -1124,13 +1124,12 @@ def restrictFromFineMesh(self,model):
     ##save model's initial solution values to archive
     def archiveInitialSolution(self,model,index):
         if True if self.fastArchive == False else 'clsvof' in model.name:
-            import xml.etree.ElementTree as ElementTree
             if self.archiveFlag == ArchiveFlags.UNDEFINED:
                 return
             logEvent("Writing initial mesh for  model = "+model.name,level=3)
             logEvent("Writing initial conditions for  model = "+model.name,level=3)
-            if not self.so.useOneArchive or index==0:
-                self.ar[index].domain = ElementTree.SubElement(self.ar[index].tree.getroot(),"Domain")
+            #the archive owns its domain (Archiver.ArchiveDomain) from
+            #construction; the writers add their collections to it
             if self.so.useOneArchive:
                 model.levelModelList[-1].archiveFiniteElementSolutions(self.ar[index],self.tnList[0],self.tCount,initialPhase=True,
                                                                        writeVectors=True,meshChanged=True,femSpaceWritten=self.femSpaceWritten,
diff --git a/proteus/mbd/CouplingFSI.pyx b/proteus/mbd/CouplingFSI.pyx
index 580a73ecbc..37055315f0 100644
--- a/proteus/mbd/CouplingFSI.pyx
+++ b/proteus/mbd/CouplingFSI.pyx
@@ -50,10 +50,33 @@ import h5py
 cimport ChronoHeaders as ch
 # chrono Python headers
 from proteus.mprans import BodyDynamics as bd
-from proteus.Archiver import indentXML
 import pychrono as chrono
 from pychrono import fea as chrono_fea
 
+
+def indentXML(elem, level=0):
+    """Pretty-print an ElementTree in place.
+
+    This module still writes its own XDMF for the bodies and cables it
+    moves (one ``.xmf`` per object, with the per-step fragments stored in
+    the HDF5 file the way proteus <= 1.9.x did). proteus's own archive no
+    longer builds XML, so the helper that used to live in Archiver lives
+    here until these writers go through the archive too.
+    """
+    i = "\n" + level*"  "
+    if len(elem):
+        if not elem.text or not elem.text.strip():
+            elem.text = i + "  "
+        for e in elem:
+            indentXML(e, level+1)
+            if not e.tail or not e.tail.strip():
+                e.tail = i + "  "
+        if not e.tail or not e.tail.strip():
+            e.tail = i
+    else:
+        if level and (not elem.tail or not elem.tail.strip()):
+            elem.tail = i
+
 # needed for sphinx docs
 __all__ = ['ProtChSystem',
            'ProtChBody',
diff --git a/test/mesh_tests/test_interpolatedbathy.py b/test/mesh_tests/test_interpolatedbathy.py
index 1a07d0489d..6922a5b66a 100644
--- a/test/mesh_tests/test_interpolatedbathy.py
+++ b/test/mesh_tests/test_interpolatedbathy.py
@@ -4,7 +4,6 @@
 from proteus.MeshTools import InterpolatedBathymetryMesh
 from proteus.Archiver import XdmfArchive
 
-import xml.etree.ElementTree as ElementTree
 import pytest
 
 @pytest.mark.MeshTools
@@ -99,7 +98,6 @@ def test_L1(self):
                                           bathyAssignmentScheme="localAveraging",errorNormType="L1")
         outfile = "interpolatedBathySimpleTest_L1_"
         archive = XdmfArchive(dataDir='.',filename = outfile, global_sync=False)
-        archive.domain = ElementTree.SubElement(archive.tree.getroot(),"Domain")
         mesh.meshList[-1].writeMeshXdmf(ar=archive,init=True)
         archive.sync(); archive.close()
         self.aux_names.append(outfile)
@@ -117,7 +115,6 @@ def test_L2(self):
                                           errorNormType="L2")
         outfile = "interpolatedBathySimpleTest_L2_"
         archive = XdmfArchive(dataDir='.',filename = outfile, global_sync=False)
-        archive.domain = ElementTree.SubElement(archive.tree.getroot(),"Domain")
         mesh.meshList[-1].writeMeshXdmf(ar=archive,init=True)
         archive.sync(); archive.close()
         self.aux_names.append(outfile)
@@ -134,7 +131,6 @@ def test_Linfty(self):
                                           errorNormType="Linfty")
         outfile = "interpolatedBathySimpleTest_Linfty_"
         archive = XdmfArchive(dataDir='.',filename=outfile, global_sync=False)
-        archive.domain = ElementTree.SubElement(archive.tree.getroot(),"Domain")
         mesh.meshList[-1].writeMeshXdmf(ar=archive,init=True)
         archive.sync(); archive.close()
         self.aux_names.append(outfile)
@@ -152,7 +148,6 @@ def test_L1_interp(self):
                                           errorNormType="L1")
         outfile = "interpolatedBathySimpleTest_L1_interp_"
         archive = XdmfArchive(dataDir='.',filename=outfile, global_sync=False)
-        archive.domain = ElementTree.SubElement(archive.tree.getroot(),"Domain")
         mesh.meshList[-1].writeMeshXdmf(ar=archive,init=True)
         archive.sync(); archive.close()
         self.aux_names.append(outfile)
@@ -170,7 +165,6 @@ def test_L2_interp(self):
                                           errorNormType="L2")
         outfile = "interpolatedBathySimpleTest_L2_interp_"
         archive = XdmfArchive(dataDir='.',filename=outfile, global_sync=False)
-        archive.domain = ElementTree.SubElement(archive.tree.getroot(),"Domain")
         mesh.meshList[-1].writeMeshXdmf(ar=archive,init=True)
         archive.sync(); archive.close()
         self.aux_names.append(outfile)
@@ -188,7 +182,6 @@ def test_Linfty_interp(self):
                                           errorNormType="Linfty")
         outfile = "interpolatedBathySimpleTest_Linfty_interp_"
         archive = XdmfArchive(dataDir='.',filename=outfile, global_sync=False)
-        archive.domain = ElementTree.SubElement(archive.tree.getroot(),"Domain")
         mesh.meshList[-1].writeMeshXdmf(ar=archive,init=True)
         archive.sync(); archive.close()
         self.aux_names.append(outfile)
@@ -206,7 +199,6 @@ def test_L1_grid(self):
                                           errorNormType="L1")
         outfile = "interpolatedBathySimpleTest_grid_L1_"
         archive = XdmfArchive(dataDir='.',filename=outfile, global_sync=False)
-        archive.domain = ElementTree.SubElement(archive.tree.getroot(),"Domain")
         mesh.meshList[-1].writeMeshXdmf(ar=archive,init=True)
         archive.sync(); archive.close()
         self.aux_names.append(outfile)
@@ -224,7 +216,6 @@ def test_L2_grid(self):
                                           errorNormType="L2")
         outfile = "interpolatedBathySimpleTest_grid_L2_"
         archive = XdmfArchive(dataDir='.',filename=outfile, global_sync=False)
-        archive.domain = ElementTree.SubElement(archive.tree.getroot(),"Domain")
         mesh.meshList[-1].writeMeshXdmf(ar=archive,init=True)
         archive.sync(); archive.close()
         self.aux_names.append(outfile)
@@ -242,7 +233,6 @@ def test_Linfty_grid(self):
                                           errorNormType="Linfty")
         outfile = "interpolatedBathySimpleTest_grid_Linfty_"
         archive = XdmfArchive(dataDir='.',filename=outfile, global_sync=False)
-        archive.domain = ElementTree.SubElement(archive.tree.getroot(),"Domain")
         mesh.meshList[-1].writeMeshXdmf(ar=archive,init=True)
         archive.sync(); archive.close()
         self.aux_names.append(outfile)
@@ -260,7 +250,6 @@ def test_L1_interp_grid(self):
                                           errorNormType="L1")
         outfile = "interpolatedBathySimpleTest_grid_L1_interp_"
         archive = XdmfArchive(dataDir='.',filename=outfile, global_sync=False)
-        archive.domain = ElementTree.SubElement(archive.tree.getroot(),"Domain")
         mesh.meshList[-1].writeMeshXdmf(ar=archive,init=True)
         archive.sync(); archive.close()
         self.aux_names.append(outfile)
@@ -278,7 +267,6 @@ def test_L2_interp_grid(self):
                                           errorNormType="L2")
         outfile = "interpolatedBathySimpleTest_grid_L2_interp_"
         archive = XdmfArchive(dataDir='.',filename=outfile, global_sync=False)
-        archive.domain = ElementTree.SubElement(archive.tree.getroot(),"Domain")
         mesh.meshList[-1].writeMeshXdmf(ar=archive,init=True)
         archive.sync(); archive.close()
         self.aux_names.append(outfile)
@@ -296,7 +284,6 @@ def test_Linfty_interp_grid(self):
                                           errorNormType="Linfty")
         outfile = "interpolatedBathySimpleTest_grid_Linfty_interp_"
         archive = XdmfArchive(dataDir='.',filename=outfile, global_sync=False)
-        archive.domain = ElementTree.SubElement(archive.tree.getroot(),"Domain")
         mesh.meshList[-1].writeMeshXdmf(ar=archive,init=True)
         archive.sync(); archive.close()
         mesh.meshList[-1].writeMeshADH("interpolatedBathySimpleTest_grid_Linfty_interp_")
diff --git a/test/test_archive_dtypes.py b/test/test_archive_dtypes.py
index f75cbe6d41..6f19735fed 100644
--- a/test/test_archive_dtypes.py
+++ b/test/test_archive_dtypes.py
@@ -13,7 +13,6 @@
 finding out afterwards with archives already written.
 """
 import os
-from xml.etree.ElementTree import Element, SubElement
 
 import numpy as np
 import pytest
@@ -58,13 +57,14 @@ def test_a_dtype_with_no_xdmf_type_is_refused_not_guessed():
 def test_topology_and_geometry_cannot_be_written_without_a_type():
     ar = Archiver.AR_base.__new__(Archiver.AR_base)
     ar.dataItemFormat, ar.hdfFile, ar.hdfFilename = "HDF", object(), "out.h5"
-    grid = Element("Grid")
+    grid = Archiver.ArchiveCollection("Mesh").new_grid(0.0, 0)
     with pytest.raises(TypeError):
         ar.write_topology(grid, "Triangle", 2, [2, 3], "elements", "elements")
     with pytest.raises(TypeError):
         ar.write_geometry(grid, [4, 3], "nodes", "nodes")
     item = ar.write_geometry(grid, [4, 3], "nodes", "nodes", dtype="float32")
-    assert (item.get("DataType"), item.get("Precision")) == ("Float", "4")
+    assert (item["data_type"], item["precision"]) == ("Float", 4)
+    assert grid["geometry"]["data_item"] is item
 
 
 # --------------------------------------------------------------------------
@@ -100,9 +100,7 @@ def __init__(self, mesh):
 @pytest.fixture
 def archive(tmp_path):
     Comm.init()
-    ar = Archiver.XdmfArchive(str(tmp_path), "dtypes", global_sync=False)
-    ar.domain = SubElement(ar.tree.getroot(), "Domain")
-    return ar
+    return Archiver.XdmfArchive(str(tmp_path), "dtypes", global_sync=False)
 
 
 def written_arrays(path):
diff --git a/test/test_archive_no_xml.py b/test/test_archive_no_xml.py
new file mode 100644
index 0000000000..ea61a21fc7
--- /dev/null
+++ b/test/test_archive_no_xml.py
@@ -0,0 +1,88 @@
+"""proteus builds its archive as ymf dicts and writes no XML.
+
+C2 of the campaign roster: only ``ymf2xmf`` should need XDMF. The archive
+(``Archiver.py``), the mesh writers (``MeshTools.py``), the finite-element
+writers (``FemTools.py``), the models that drive them and the auxiliary
+variables build ymf grid dicts through the archive's builders; the ``.xmf``
+a viewer wants is derived from the ``.ymf`` on demand. These tests pin that
+no module on the archive path imports an XML library, and that the running
+archive is one step of dicts, not a tree.
+
+The one writer that still builds XDMF, ``mbd/CouplingFSI.pyx``, keeps its
+own files for the bodies it moves and is outside this check.
+"""
+import ast
+import os
+
+import numpy as np
+import pytest
+
+PROTEUS = os.path.join(os.path.dirname(os.path.dirname(os.path.abspath(__file__))), "proteus")
+
+ARCHIVE_PATH = ["Archiver.py", "MeshTools.py", "FemTools.py", "AuxiliaryVariables.py",
+                "NumericalSolution.py", "Transport.py", "SplitOperator.py"]
+
+
+def imports_of(path):
+    tree = ast.parse(open(path).read(), path)
+    found = set()
+    for node in ast.walk(tree):
+        if isinstance(node, ast.Import):
+            found.update(alias.name for alias in node.names)
+        elif isinstance(node, ast.ImportFrom) and node.module:
+            found.add(node.module)
+    return found
+
+
+@pytest.mark.parametrize("module", ARCHIVE_PATH)
+def test_the_archive_path_imports_no_xml_library(module):
+    names = imports_of(os.path.join(PROTEUS, module))
+    assert not [n for n in names if n == "xml" or n.startswith("xml.")], module
+
+
+def test_no_module_of_the_package_builds_xdmf_elements():
+    """``SubElement`` was the signature of the hand-built XDMF; none is left."""
+    offenders = []
+    for root, _, files in os.walk(PROTEUS):
+        for name in files:
+            if name.endswith(".py"):
+                path = os.path.join(root, name)
+                with open(path) as f:
+                    if "SubElement(" in f.read():
+                        offenders.append(os.path.relpath(path, PROTEUS))
+    assert offenders == []
+
+
+def test_the_running_archive_is_one_step_of_dicts():
+    from proteus.Archiver import ArchiveCollection, ArchiveDomain
+
+    domain = ArchiveDomain()
+    mesh = domain.collection("Mesh Spatial_Domain")
+    space = domain.collection("Mesh_c0p2_Lagrange")
+    assert domain.collection("Mesh Spatial_Domain") is mesh      # found, not duplicated
+    assert domain.collection("Mesh_dgp1", create=False) is None   # a continuing run expects it
+    assert domain.names() == ["Mesh Spatial_Domain", "Mesh_c0p2_Lagrange"]
+    assert domain[0] is mesh and len(domain) == 2
+    g = mesh.new_grid(0.25, 1)
+    assert mesh.is_current(0.25) and not space.is_current(0.25)
+    domain.clear_step()
+    assert mesh.grid is None and not mesh.is_current(0.25)
+    assert isinstance(g, dict) and g["attributes"] == []
+
+
+def test_a_grid_scalar_is_inlined_in_the_document():
+    """AuxiliaryVariables.writeScalarXdmf used to put the value in a DataItem's text."""
+    from proteus.Archiver import ArchiveDomain
+    from proteus.AuxiliaryVariables import AV_base
+
+    class Archive(object):
+        domain = ArchiveDomain()
+
+    av = AV_base()
+    av.ar = Archive()
+    av.ar.domain.collection("Mesh Spatial_Domain").new_grid(0.0, 0)
+    av.writeScalarXdmf(3.5, "force")
+    (attr,) = av.ar.domain[0].grid["attributes"]
+    assert attr["name"] == "force" and attr["center"] == "Grid"
+    assert attr["data_item"]["values"] == [3.5]
+    assert attr["data_item"]["format"] == "XML"
diff --git a/test/test_archive_parallel_metadata.py b/test/test_archive_parallel_metadata.py
index 08a34ae4f2..6c23b08a49 100644
--- a/test/test_archive_parallel_metadata.py
+++ b/test/test_archive_parallel_metadata.py
@@ -54,7 +54,6 @@ def _find_mpiexec():
 RANK_SCRIPT = textwrap.dedent(
     '''
     import os, shutil, sys
-    from xml.etree.ElementTree import SubElement
     import numpy as np
 
     from proteus import Comm
@@ -79,23 +78,18 @@ def _find_mpiexec():
     if hot:
         # what NumericalSolution does: continue after the archive's last step
         ar.n_datasets = len(ar.archived_times)
-        ar.domain = ar.tree.getroot().find("Domain")
-    else:
-        ar.domain = SubElement(ar.tree.getroot(), "Domain")
 
     # Sizes differ per rank on purpose: the metadata dataset is fixed-width
     # and created collectively, so if the width agreement is wrong the
     # longest rank's YAML is silently truncated.
     n_elements, n_nodes = 4 + rank, 6 + rank
     if hot:
-        # found by name, as XdmfWriter.setGridCollectionAndGridElements does
-        (collection,) = [c for c in ar.domain
-                         if c.attrib["Name"] == "Mesh Spatial_Domain"]
+        # found by name, as XdmfWriter.setGridCollectionAndGridElements does;
+        # the archive seeded it from the .h5 when it was opened
+        collection = ar.domain.collection("Mesh Spatial_Domain", create=False)
+        assert collection is not None, ar.domain.names()
     else:
-        collection = SubElement(ar.domain, "Grid",
-                                {"Name": "Mesh Spatial_Domain",
-                                 "GridType": "Collection",
-                                 "CollectionType": "Temporal"})
+        collection = ar.domain.collection("Mesh Spatial_Domain")
 
     # global_sync describes one assembled global array, so the metadata
     # carries global counts and each rank contributes only the slice it owns.
@@ -107,7 +101,7 @@ def _find_mpiexec():
 
     times = [1.0] if hot else [0.0, 0.5]
     for tCount, t in enumerate(times, start=ar.n_datasets):
-        grid, _ = ar.write_grid(collection, "Grid_p%d" % rank, t, tCount)
+        grid = ar.write_grid(collection, "Grid_p%d" % rank, t, tCount)
         if global_sync:
             ar.write_topology(grid, "Triangle", n_elements_global,
                               [n_elements_global, 3],
diff --git a/test/test_archive_read_paths.py b/test/test_archive_read_paths.py
index a9f0848ce6..af52715f0c 100644
--- a/test/test_archive_read_paths.py
+++ b/test/test_archive_read_paths.py
@@ -5,9 +5,14 @@
 first child as its Time. Those positions are not guaranteed by XDMF:
 ``<Information>`` is legal ahead of ``<Domain>``, a real proteus archive
 commonly carries two grid collections (the linear mesh and a quadratic
-space), and nothing fixes the order of a grid's children. Each of these
-tests puts the thing being looked for somewhere other than the position the
-old code assumed, so a positional reader fails them.
+space), and nothing fixes the order of a grid's children.
+
+The readers now take the archive's own domain dict (what
+``Archiver.readArchiveDomain`` assembles from the ``.h5``), so the
+positional traps are gone with the tree; what these tests pin is the
+lookup by name, the step selection and the dataset references. The
+``gatherArchives`` script still composes ``.xmf`` files, and keeps its
+XML-built fixtures.
 """
 
 import os
@@ -25,9 +30,42 @@
 TWO_COLLECTIONS = ("Mesh Spatial_Domain", "Mesh_c0p2_Lagrange")
 
 
+def build_domain(n_steps=2, collections=TWO_COLLECTIONS, rank=0, spatial=False):
+    """An archive domain shaped like proteus's: two spaces, a few steps."""
+    def grid(tn, name=None):
+        g = {}
+        if name is not None:
+            g["name"] = name
+        g["topology"] = {"type": "Triangle", "number_of_elements": 2,
+                         "data_item": {"format": "HDF", "data_type": "Int",
+                                       "precision": 4, "dimensions": [2, 3],
+                                       "data": "f.h5:/e%d_p%d" % (tn, rank)}}
+        g["geometry"] = {"type": "XYZ",
+                         "data_item": {"format": "HDF", "data_type": "Float",
+                                       "precision": 8, "dimensions": [4, 3],
+                                       "data": "f.h5:/n%d_p%d" % (tn, rank)}}
+        g["attributes"] = [{"name": "nodeMaterialTypes", "attribute_type": "Scalar",
+                            "center": "Node",
+                            "data_item": {"format": "HDF", "data_type": "Int",
+                                          "precision": 4, "dimensions": [4],
+                                          "data": "f.h5:/nodeMat%d" % tn}}]
+        return g
+    domain = {"time_collections": []}
+    for name in collections:
+        steps = []
+        for tn in range(n_steps):
+            if spatial:
+                steps.append({"time": float(tn),
+                              "spatial_collection": [grid(tn, "p%d" % r) for r in range(2)]})
+            else:
+                steps.append(dict(time=float(tn), **grid(tn)))
+        domain["time_collections"].append({"name": name, "data": steps})
+    return domain
+
+
 def build_xmf(path, n_steps=2, collections=TWO_COLLECTIONS, rank=0,
               lead_with_information=False, time_last=False):
-    """A per-rank .xmf shaped like proteus's, with optional awkwardness.
+    """A per-rank .xmf shaped like proteus's, for the gatherArchives script.
 
     ``lead_with_information`` puts an <Information> element ahead of
     <Domain>, which XDMF allows and which breaks ``getroot()[-1]``.
@@ -69,107 +107,76 @@ def build_xmf(path, n_steps=2, collections=TWO_COLLECTIONS, rank=0,
 
 
 # --------------------------------------------------------------------------
-# MeshTools.findXMLgridElement
+# MeshTools.findArchiveGrid
 # --------------------------------------------------------------------------
 
 
-def test_the_named_collection_is_found_not_the_last_one(tmp_path):
+def test_the_named_collection_is_found_not_the_last_one():
     """`Domain[-1]` would return the quadratic space, not the base mesh."""
-    from proteus.MeshTools import findXMLgridElement
+    from proteus.MeshTools import findArchiveGrid
 
-    path = build_xmf(str(tmp_path / "a.xmf"))
-    grid = findXMLgridElement(parse(path), MeshTag="Spatial_Domain")
+    grid = findArchiveGrid(build_domain(), MeshTag="Spatial_Domain")
     # the base mesh's connectivity, not the c0p2 space's
-    assert grid.find("Topology").find("DataItem").text.endswith("e1_p0")
-    assert grid.attrib["GridType"] == "Uniform"
-
-
-def test_an_information_element_before_domain_is_tolerated(tmp_path):
-    """XDMF permits <Information*> ahead of <Grid+> inside <Xdmf>.
+    assert grid["topology"]["data_item"]["data"].endswith("e1_p0")
 
-    `getroot()[-1]` happened to work only because nothing was written
-    after <Domain>.
-    """
-    from proteus.MeshTools import findXMLgridElement
 
-    path = build_xmf(str(tmp_path / "a.xmf"), lead_with_information=True)
-    grid = findXMLgridElement(parse(path), MeshTag="Spatial_Domain")
-    assert grid.attrib["GridType"] == "Uniform"
+def test_an_unknown_mesh_tag_falls_back_to_the_first_collection():
+    from proteus.MeshTools import findArchiveGrid
 
+    grid = findArchiveGrid(build_domain(), MeshTag="NoSuchMesh")
+    assert grid["topology"]["data_item"]["data"].endswith("e1_p0")
 
-def test_an_unknown_mesh_tag_falls_back_to_the_first_collection(tmp_path):
-    from proteus.MeshTools import findXMLgridElement
 
-    path = build_xmf(str(tmp_path / "a.xmf"))
-    grid = findXMLgridElement(parse(path), MeshTag="NoSuchMesh")
-    assert grid.attrib["GridType"] == "Uniform"
+def test_a_step_can_be_selected_by_index():
+    from proteus.MeshTools import findArchiveGrid
 
+    domain = build_domain(n_steps=3)
+    first = findArchiveGrid(domain, id_in_collection=0)
+    last = findArchiveGrid(domain, id_in_collection=-1)
+    assert first["topology"]["data_item"]["data"].endswith("e0_p0")
+    assert last["topology"]["data_item"]["data"].endswith("e2_p0")
 
-def test_a_step_can_be_selected_by_index(tmp_path):
-    from proteus.MeshTools import findXMLgridElement
 
-    path = build_xmf(str(tmp_path / "a.xmf"), n_steps=3)
-    first = findXMLgridElement(parse(path), id_in_collection=0)
-    last = findXMLgridElement(parse(path), id_in_collection=-1)
-    assert first.find("Topology").find("DataItem").text.endswith("e0_p0")
-    assert last.find("Topology").find("DataItem").text.endswith("e2_p0")
+def test_a_spatial_collection_step_resolves_to_a_uniform_grid():
+    """A per-subdomain archive holds one grid per rank in each step."""
+    from proteus.MeshTools import findArchiveGrid
 
+    grid = findArchiveGrid(build_domain(spatial=True))
+    assert grid["name"] == "p0"
+    assert "spatial_collection" not in grid
 
-def test_a_spatial_collection_step_resolves_to_a_uniform_grid(tmp_path):
-    """A per-subdomain archive nests Grid(Spatial) > Grid(Uniform)."""
-    from proteus.MeshTools import findXMLgridElement
 
-    root = Element("Xdmf", {"Version": "2.0"})
-    domain = SubElement(root, "Domain")
-    coll = SubElement(domain, "Grid",
-                      {"Name": "Mesh Spatial_Domain", "GridType": "Collection",
-                       "CollectionType": "Temporal"})
-    spatial = SubElement(coll, "Grid", {"GridType": "Collection",
-                                        "CollectionType": "Spatial"})
-    SubElement(spatial, "Time", {"Value": "0.0", "Name": "0"})
-    for rank in range(2):
-        g = SubElement(spatial, "Grid",
-                       {"GridType": "Uniform", "Name": "p%d" % rank})
-        topo = SubElement(g, "Topology",
-                          {"Type": "Triangle", "NumberOfElements": "2"})
-        SubElement(topo, "DataItem",
-                   {"Format": "HDF", "DataType": "Int",
-                    "Dimensions": "2 3"}).text = "f.h5:/e_p%d" % rank
-    path = str(tmp_path / "s.xmf")
-    with open(path, "wb") as handle:
-        ElementTree(root).write(handle)
+def test_an_empty_archive_is_reported():
+    from proteus.MeshTools import findArchiveGrid
 
-    grid = findXMLgridElement(parse(path))
-    assert grid.attrib["GridType"] == "Uniform"
-    assert grid.attrib["Name"] == "p0"
+    with pytest.raises(AssertionError, match="no time collections"):
+        findArchiveGrid({"time_collections": []})
+    with pytest.raises(AssertionError, match="holds no grids"):
+        findArchiveGrid({"time_collections": [{"name": "Mesh", "data": []}]})
 
 
 # --------------------------------------------------------------------------
-# MeshTools.extractPropertiesFromXdmfGridNode and the DataItem reads
+# MeshTools.extractPropertiesFromGrid and the DataItem reads
 # --------------------------------------------------------------------------
 
 
-def test_grid_properties_are_found_regardless_of_child_order(tmp_path):
-    """`Grid[0]` assumed Time came first; nothing in XDMF says so."""
-    from proteus.MeshTools import extractPropertiesFromXdmfGridNode, findXMLgridElement
+def test_grid_properties_are_found_by_name():
+    from proteus.MeshTools import extractPropertiesFromGrid, findArchiveGrid
 
-    path = build_xmf(str(tmp_path / "a.xmf"), time_last=True)
-    grid = findXMLgridElement(parse(path))
+    grid = findArchiveGrid(build_domain())
     topology, geometry, node_materials, element_materials = \
-        extractPropertiesFromXdmfGridNode(grid)
-    assert topology is not None and topology.tag == "Topology"
-    assert geometry is not None and geometry.tag == "Geometry"
-    assert node_materials is not None
-    assert node_materials.attrib["Name"] == "nodeMaterialTypes"
+        extractPropertiesFromGrid(grid)
+    assert topology["type"] == "Triangle"
+    assert geometry["type"] == "XYZ"
+    assert node_materials["name"] == "nodeMaterialTypes"
     assert element_materials is None       # not written by this fixture
 
 
-def test_a_dataitem_reference_resolves_to_its_dataset_name(tmp_path):
-    from proteus.MeshTools import _dataitem_dataset, findXMLgridElement
+def test_a_dataitem_reference_resolves_to_its_dataset_name():
+    from proteus.MeshTools import _dataitem_dataset, findArchiveGrid
 
-    path = build_xmf(str(tmp_path / "a.xmf"))
-    grid = findXMLgridElement(parse(path))
-    assert _dataitem_dataset(grid.find("Geometry"), "Geometry") == "n1_p0"
+    grid = findArchiveGrid(build_domain())
+    assert _dataitem_dataset(grid["geometry"], "Geometry") == "n1_p0"
 
 
 def test_a_dataset_path_containing_a_colon_survives():
@@ -180,16 +187,17 @@ def test_a_dataset_path_containing_a_colon_survives():
     """
     from proteus.MeshTools import _dataitem_dataset
 
-    geometry = Element("Geometry")
-    SubElement(geometry, "DataItem").text = "out.h5:/group/od:d"
+    geometry = {"type": "XYZ", "data_item": {"data": "out.h5:/group/od:d"}}
     assert _dataitem_dataset(geometry, "Geometry") == "group/od:d"
 
 
-def test_a_missing_dataitem_is_reported(tmp_path):
+def test_a_missing_dataitem_is_reported():
     from proteus.MeshTools import _dataitem_dataset
 
-    with pytest.raises(AssertionError, match="no <DataItem>"):
-        _dataitem_dataset(Element("Geometry"), "Geometry")
+    with pytest.raises(AssertionError, match="no data_item"):
+        _dataitem_dataset({"type": "XYZ"}, "Geometry")
+    with pytest.raises(AssertionError, match="no HDF5 reference"):
+        _dataitem_dataset({"data_item": {"include": "./x.txt"}}, "Geometry")
 
 
 # --------------------------------------------------------------------------
diff --git a/test/test_archiver_write_field.py b/test/test_archiver_write_field.py
index f7ffef9db1..3582078e7a 100644
--- a/test/test_archiver_write_field.py
+++ b/test/test_archiver_write_field.py
@@ -1,20 +1,21 @@
-"""Tests for AR_base.write_field, the Attribute+DataItem chokepoint.
+"""Tests for AR_base.write_field, the attribute-and-data-item chokepoint.
 
 write_field replaces an eight-line block that was hand-written at ~50 call
-sites. These tests pin its XML against the block it replaces, so the
-conversion of each call site is a refactor and not a rewrite.
+sites. It used to build an XDMF ``Attribute``/``DataItem`` pair; it now
+builds the ymf attribute dict the archive stores, and these tests pin that
+dict -- the data reference, the declared dimensions and type, the sidecar
+path in text mode -- so the conversion of each call site is a refactor and
+not a rewrite.
 
 They build a stub archive rather than a real XdmfArchive: the real one opens
 an MPI-mode HDF5 file in its constructor, which is not what is under test
 here.
 """
 
-from xml.etree.ElementTree import Element, tostring
-
 import numpy as np
 import pytest
 
-from proteus.Archiver import AR_base
+from proteus.Archiver import AR_base, ArchiveCollection
 
 
 class StubComm(object):
@@ -65,49 +66,39 @@ def fake_savetxt(path, arr, *a, **kw):
 
 
 def grid():
-    return Element("Grid", {"GridType": "Uniform"})
+    """A grid as write_grid starts one: nothing on it yet."""
+    return ArchiveCollection("Mesh Spatial_Domain").new_grid(0.0, 0)
+
+
+def the_attribute(g):
+    (attr,) = g["attributes"]
+    return attr
+
+
+def the_item(g):
+    return the_attribute(g)["data_item"]
 
 
 # --------------------------------------------------------------------------
-# equivalence with the block write_field replaces
+# equivalence with what ymf's own builders produce
 # --------------------------------------------------------------------------
 
 
-def reference_block(ar, grid_elem, name, dof, tCount, n_declared):
-    """The hand-written block, verbatim in structure, for comparison.
+def test_write_field_builds_the_attribute_ymf_would():
+    """The dict is the one ymf.archive.attribute/data_item_for build.
 
-    This mirrors what FemTools.writeFunctionXdmf did before the conversion:
-    an Attribute, a DataItem with Format/DataType/Precision/Dimensions, and
-    the HDF5 reference as the DataItem's text.
+    That is what makes the archive one representation with one writer:
+    a field written by proteus is indistinguishable from one described
+    by hand with ymf's builders.
     """
-    from xml.etree.ElementTree import SubElement
-
-    attribute = SubElement(grid_elem, "Attribute",
-                           {"Name": name,
-                            "AttributeType": "Scalar",
-                            "Center": "Node"})
-    values = SubElement(attribute, "DataItem",
-                        {"Format": ar.dataItemFormat,
-                         "DataType": "Float",
-                         "Precision": "8",
-                         "Dimensions": "%i" % (n_declared,)})
-    values.text = (ar.hdfFilename + ":/" + name + "_p"
-                   + repr(ar.comm.rank()) + "_t{0:d}".format(tCount))
-    return values
-
-
-def test_write_field_emits_the_same_xml_as_the_block_it_replaces():
-    dof = np.zeros(17, dtype="float64")
-
-    old_ar = StubArchive()
-    old_grid = grid()
-    reference_block(old_ar, old_grid, "u", dof, 3, 17)
+    from ymf.archive import attribute, data_item_for
 
-    new_ar = StubArchive()
-    new_grid = grid()
-    new_ar.write_field(new_grid, "u", dof, 3, dimensions=[17])
-
-    assert tostring(new_grid) == tostring(old_grid)
+    dof = np.zeros(17, dtype="float64")
+    g = grid()
+    StubArchive().write_field(g, "u", dof, 3, dimensions=[17])
+    expected = attribute("u", data_item_for(dof, data="out.h5:/u_p0_t3", dimensions=[17]),
+                         center="Node", attribute_type="Scalar")
+    assert the_attribute(g) == expected
 
 
 def test_the_dataset_name_matches_the_old_convention():
@@ -121,23 +112,35 @@ def test_the_hdf5_reference_matches_the_dataset_name():
     ar = StubArchive(rank=1)
     g = grid()
     ar.write_field(g, "vof", np.zeros(4), 2, dimensions=[4])
-    item = g.find("Attribute").find("DataItem")
-    assert item.text == "out.h5:/vof_p1_t2"
+    assert the_item(g)["data"] == "out.h5:/vof_p1_t2"
     assert ar.async_calls[0][0] == "vof_p1_t2"
 
 
+def test_the_item_returned_is_the_one_attached():
+    g = grid()
+    item = StubArchive().write_field(g, "u", np.zeros(4), 0)
+    assert item is the_item(g)
+
+
+def test_fields_accumulate_on_the_grid_in_order():
+    g = grid()
+    ar = StubArchive()
+    ar.write_field(g, "p", np.zeros(4), 0)
+    ar.write_field(g, "u", np.zeros(4), 0)
+    assert [a["name"] for a in g["attributes"]] == ["p", "u"]
+
+
 # --------------------------------------------------------------------------
-# the Attribute itself
+# the attribute itself
 # --------------------------------------------------------------------------
 
 
 def test_defaults_are_node_centered_scalar():
     g = grid()
-    ar = StubArchive()
-    ar.write_field(g, "u", np.zeros(3), 0, dimensions=[3])
-    attr = g.find("Attribute")
-    assert attr.attrib["AttributeType"] == "Scalar"
-    assert attr.attrib["Center"] == "Node"
+    StubArchive().write_field(g, "u", np.zeros(3), 0, dimensions=[3])
+    attr = the_attribute(g)
+    assert attr["attribute_type"] == "Scalar"
+    assert attr["center"] == "Node"
 
 
 @pytest.mark.parametrize("center", ["Node", "Cell", "Face"])
@@ -146,20 +149,20 @@ def test_center_and_rank_reach_the_attribute(center, rank):
     g = grid()
     StubArchive().write_field(g, "u", np.zeros(3), 0, center=center,
                               rank=rank, dimensions=[3])
-    attr = g.find("Attribute")
-    assert (attr.attrib["Center"], attr.attrib["AttributeType"]) == (center, rank)
+    attr = the_attribute(g)
+    assert (attr["center"], attr["attribute_type"]) == (center, rank)
 
 
-def test_multi_dimensional_dimensions_are_space_separated():
+def test_multi_dimensional_dimensions_are_kept():
     g = grid()
     StubArchive().write_field(g, "vel", np.zeros((5, 3)), 0, rank="Vector")
-    assert g.find("Attribute").find("DataItem").attrib["Dimensions"] == "5 3"
+    assert list(the_item(g)["dimensions"]) == [5, 3]
 
 
 def test_dimensions_default_to_the_array_shape():
     g = grid()
     StubArchive().write_field(g, "u", np.zeros(11), 0)
-    assert g.find("Attribute").find("DataItem").attrib["Dimensions"] == "11"
+    assert list(the_item(g)["dimensions"]) == [11]
 
 
 # --------------------------------------------------------------------------
@@ -170,8 +173,8 @@ def test_dimensions_default_to_the_array_shape():
 def test_float64_is_described_as_eight_byte_float():
     g = grid()
     StubArchive().write_field(g, "u", np.zeros(3, dtype="float64"), 0)
-    item = g.find("Attribute").find("DataItem")
-    assert (item.attrib["DataType"], item.attrib["Precision"]) == ("Float", "8")
+    item = the_item(g)
+    assert (item["data_type"], item["precision"]) == ("Float", 8)
 
 
 def test_float32_is_described_as_four_byte_float():
@@ -179,15 +182,15 @@ def test_float32_is_described_as_four_byte_float():
     # therefore described to consumers as 8-byte, which reads as garbage.
     g = grid()
     StubArchive().write_field(g, "u", np.zeros(3, dtype="float32"), 0)
-    item = g.find("Attribute").find("DataItem")
-    assert (item.attrib["DataType"], item.attrib["Precision"]) == ("Float", "4")
+    item = the_item(g)
+    assert (item["data_type"], item["precision"]) == ("Float", 4)
 
 
 def test_integer_arrays_are_described_as_int():
     g = grid()
     StubArchive().write_field(g, "mat", np.zeros(3, dtype="int32"), 0, center="Cell")
-    item = g.find("Attribute").find("DataItem")
-    assert (item.attrib["DataType"], item.attrib["Precision"]) == ("Int", "4")
+    item = the_item(g)
+    assert (item["data_type"], item["precision"]) == ("Int", 4)
 
 
 # --------------------------------------------------------------------------
@@ -204,7 +207,7 @@ def test_global_sync_declares_the_global_size_not_the_slice_size():
     offsets = np.array([0, 4, 10])
     ar.write_field(g, "u", dof, 1, dimensions=[10], sync_offsets=offsets,
                    sync_data=dof[:4])
-    assert g.find("Attribute").find("DataItem").attrib["Dimensions"] == "10"
+    assert list(the_item(g)["dimensions"]) == [10]
     (name, got_offsets, got_data), = ar.sync_calls
     assert name == "u_t1"          # no _p<rank> in the sync convention
     np.testing.assert_array_equal(got_offsets, offsets)
@@ -216,7 +219,7 @@ def test_global_sync_uses_the_unranked_dataset_name():
     g = grid()
     ar.write_field(g, "p", np.zeros(4), 5, dimensions=[4], sync_offsets=np.array([0, 4]),
                    sync_data=np.zeros(4))
-    assert g.find("Attribute").find("DataItem").text == "out.h5:/p_t5"
+    assert the_item(g)["data"] == "out.h5:/p_t5"
     assert not ar.async_calls
 
 
@@ -225,16 +228,14 @@ def test_global_sync_uses_the_unranked_dataset_name():
 # --------------------------------------------------------------------------
 
 
-def test_text_mode_attaches_an_xi_include_to_its_own_data_item(_no_real_savetxt):
+def test_text_mode_references_its_own_sidecar(_no_real_savetxt):
     ar = StubArchive(hdf=False)
     g = grid()
     ar.write_field(g, "u", np.zeros(4), 2)
-    item = g.find("Attribute").find("DataItem")
-    assert item.attrib["Format"] == "XML"
-    (include,) = list(item)
-    assert include.tag == "xi:include"
-    assert include.attrib == {"parse": "text", "href": "./out_Data/u2.txt"}
-    assert not (item.text or "").strip()
+    item = the_item(g)
+    assert item["format"] == "XML"
+    assert item["include"] == "./out_Data/u2.txt"
+    assert "data" not in item
 
 
 def test_text_mode_writes_the_sidecar_file(_no_real_savetxt):
@@ -246,23 +247,18 @@ def test_text_mode_writes_the_sidecar_file(_no_real_savetxt):
     np.testing.assert_array_equal(arr, dof)
 
 
-def test_two_fields_in_text_mode_each_get_exactly_one_include(_no_real_savetxt):
+def test_two_fields_in_text_mode_each_reference_their_own_sidecar(_no_real_savetxt):
     # This is the bug at Archiver.py:1475 made unrepresentable: there, the
     # include for a second field was attached to the first field's
     # DataItem, leaving one field with no data reference and the other with
-    # two. write_field always attaches to the DataItem it just created.
+    # two. The reference is now part of the data item write_field creates.
     ar = StubArchive(hdf=False)
     g = grid()
     ar.write_field(g, "u", np.zeros(4), 0)
     ar.write_field(g, "u_dof", np.zeros(4), 0, center="Face")
-    items = [a.find("DataItem") for a in g.findall("Attribute")]
-    assert len(items) == 2
-    hrefs = []
-    for item in items:
-        includes = [c for c in item if c.tag == "xi:include"]
-        assert len(includes) == 1, "expected exactly one xi:include per DataItem"
-        hrefs.append(includes[0].attrib["href"])
-    assert hrefs == ["./out_Data/u0.txt", "./out_Data/u_dof0.txt"]
+    items = [a["data_item"] for a in g["attributes"]]
+    assert [item["include"] for item in items] == ["./out_Data/u0.txt",
+                                                   "./out_Data/u_dof0.txt"]
 
 
 def test_text_mode_with_global_sync_is_refused():
@@ -290,7 +286,7 @@ def test_hdf_mode_never_touches_text_data_dir():
     del ar.textDataDir
     g = grid()
     ar.write_field(g, "u", np.zeros(4), 0, dimensions=[4])
-    assert g.find("Attribute").find("DataItem").text == "out.h5:/u_p0_t0"
+    assert the_item(g)["data"] == "out.h5:/u_p0_t0"
 
 
 def test_hdf_mode_with_global_sync_never_touches_text_data_dir():
@@ -299,11 +295,11 @@ def test_hdf_mode_with_global_sync_never_touches_text_data_dir():
     g = grid()
     ar.write_field(g, "u", np.zeros(4), 0, dimensions=[4],
                    sync_offsets=np.array([0, 4]), sync_data=np.zeros(4))
-    assert g.find("Attribute").find("DataItem").text == "out.h5:/u_t0"
+    assert the_item(g)["data"] == "out.h5:/u_t0"
 
 
 # --------------------------------------------------------------------------
-# dataset / text_stem overrides
+# dataset / text_stem overrides, and a reference without a write
 # --------------------------------------------------------------------------
 
 
@@ -314,7 +310,7 @@ def test_an_explicit_dataset_name_overrides_the_convention():
     ar = StubArchive(rank=2)
     g = grid()
     ar.write_field(g, "u", np.zeros(4), 7, dimensions=[4], dataset="u7")
-    assert g.find("Attribute").find("DataItem").text == "out.h5:/u7"
+    assert the_item(g)["data"] == "out.h5:/u7"
     assert ar.async_calls[0][0] == "u7"
 
 
@@ -330,11 +326,21 @@ def test_text_stem_overrides_the_sidecar_filename(_no_real_savetxt):
     ar = StubArchive(hdf=False)
     g = grid()
     ar.write_field(g, "u", np.zeros(4), 3, text_stem="custom")
-    include = list(g.find("Attribute").find("DataItem"))[0]
-    assert include.attrib["href"] == "./out_Data/custom.txt"
+    assert the_item(g)["include"] == "./out_Data/custom.txt"
     assert _no_real_savetxt[0][0] == "out_Data/custom.txt"
 
 
+def test_a_reference_can_be_attached_without_writing(_no_real_savetxt):
+    # the mesh writers reference the local-to-global maps on every pass
+    # but write them only when the mesh changed
+    ar = StubArchive(rank=1)
+    g = grid()
+    ar.write_field(g, "NodeMapL2G", np.zeros(4, dtype="int32"), 2,
+                   dataset="nodeMapL2G1Spatial_Domain2", write_data=False)
+    assert the_item(g)["data"] == "out.h5:/nodeMapL2G1Spatial_Domain2"
+    assert not ar.async_calls and not ar.sync_calls and not _no_real_savetxt
+
+
 # --------------------------------------------------------------------------
 # the declared-vs-actual dimension check, and where it must not apply
 # --------------------------------------------------------------------------
@@ -368,7 +374,7 @@ def test_a_larger_declared_size_is_accepted_on_a_parallel_collective_write():
     ar.write_field(g, "u", np.zeros(6), 0, dimensions=[26],
                    sync_offsets=np.array([0, 6, 13, 19, 26]),
                    sync_data=np.zeros(6))
-    assert g.find("Attribute").find("DataItem").attrib["Dimensions"] == "26"
+    assert list(the_item(g)["dimensions"]) == [26]
 
 
 def test_the_per_rank_path_keeps_the_check_at_any_size():
@@ -379,3 +385,85 @@ def test_the_per_rank_path_keeps_the_check_at_any_size():
     ar = StubArchive(global_sync=False, size=8)
     with pytest.raises(YmfArchiveError, match="dimensions"):
         ar.write_field(grid(), "u", np.zeros(10), 0, dimensions=[99])
+
+
+# --------------------------------------------------------------------------
+# the grid builders: write_grid, write_topology, write_geometry
+# --------------------------------------------------------------------------
+
+
+def test_write_grid_starts_the_collections_current_grid():
+    ar = StubArchive()
+    collection = ArchiveCollection("Mesh_c0p2_Lagrange")
+    g = ar.write_grid(collection, "Grid_c0p2_Lagrange", 0.5, 3)
+    assert collection.grid is g
+    assert collection.is_current(0.5) and not collection.is_current(0.75)
+    assert g["name"] == "Grid_c0p2_Lagrange"
+    assert g["attributes"] == [] and g["topology"] is None
+
+
+def test_the_collection_forgets_its_grid_when_the_step_is_stored():
+    collection = ArchiveCollection("Mesh Spatial_Domain")
+    collection.new_grid(0.0, 0)
+    collection.clear()
+    assert collection.grid is None and not collection.is_current(0.0)
+
+
+def test_topology_and_geometry_reference_the_named_datasets():
+    ar = StubArchive(global_sync=True)
+    g = grid()
+    topo = ar.write_topology(g, "Tri_6", 2, [2, 6], "elements_c0p2_Lagrange0",
+                             "elements_c0p2_Lagrange0", dtype=np.zeros(1, "int32"))
+    geo = ar.write_geometry(g, [9, 3], "nodes_c0p2_Lagrange0",
+                            "nodes_c0p2_Lagrange0", dtype=np.zeros(1, "float64"))
+    assert g["topology"]["type"] == "Tri_6"
+    assert g["topology"]["number_of_elements"] == 2
+    assert g["topology"]["data_item"] is topo
+    assert (topo["data"], topo["data_type"], topo["precision"]) == \
+        ("out.h5:/elements_c0p2_Lagrange0", "Int", 4)
+    assert g["geometry"]["type"] == "XYZ" and g["geometry"]["data_item"] is geo
+    assert (geo["data"], geo["data_type"], geo["precision"]) == \
+        ("out.h5:/nodes_c0p2_Lagrange0", "Float", 8)
+    # the builders describe; the caller writes the arrays
+    assert not ar.sync_calls and not ar.async_calls
+
+
+def test_a_polyvertex_topology_says_how_many_nodes_per_element():
+    g = grid()
+    StubArchive().write_topology(g, "Polyvertex", 4, [4, 8], "elements_q0", "elements_q0",
+                                 nodes_per_element=8, dtype=np.dtype("int32"))
+    assert g["topology"]["nodes_per_element"] == 8
+
+
+def test_mesh_data_items_in_text_mode_reference_the_sidecar():
+    ar = StubArchive(hdf=False)
+    g = grid()
+    ar.write_topology(g, "Triangle", 2, [2, 3], "elements0", "elementsSpatial_Domain",
+                      dtype=np.dtype("int32"))
+    assert g["topology"]["data_item"]["include"] == "./out_Data/elementsSpatial_Domain.txt"
+
+
+def test_a_computed_connectivity_is_declared_as_numpy_types_the_expression():
+    # offset + l2g, as the DG writers compute their elements
+    g = grid()
+    StubArchive().write_topology(g, "Triangle", 2, [2, 3], "e", "e",
+                                 dtype=(np.zeros(1, "int64"), np.zeros(1, "int32")))
+    assert g["topology"]["data_item"]["precision"] == 8
+
+
+def test_a_grid_the_writers_build_is_a_grid_ymf_accepts():
+    """The whole point: the dict goes into the archive as it is."""
+    from ymf.archive import canonicalize_domain, new_domain, validate_domain
+
+    ar = StubArchive()
+    collection = ArchiveCollection("Mesh Spatial_Domain")
+    g = ar.write_grid(collection, None, 0.0, 0)
+    ar.write_topology(g, "Triangle", 2, [2, 3], "elements0", "elements",
+                      dtype=np.dtype("int32"))
+    ar.write_geometry(g, [4, 3], "nodes0", "nodes", dtype=np.dtype("float64"))
+    ar.write_field(g, "u", np.zeros(4), 0)
+    domain = new_domain("Mesh Spatial_Domain")
+    domain["time_collections"][0]["data"].append(
+        {"time": collection.time, "topology": g["topology"],
+         "geometry": g["geometry"], "attributes": g["attributes"]})
+    validate_domain(canonicalize_domain(domain))