Repository: RangamaniLabUCSD/smart (language: python)

Issue:
Update implicit curve function and fix diffusion units check
Two minor changes:
1. Add diffusion units check to initialization rather than during variational form assembly.
2. Allow user to choose a certain number of points used to interpolate along the implicit expression for the boundary of a mesh.


Discussion / hints:

Fix patch:
diff --git a/examples/example2/example2.ipynb b/examples/example2/example2.ipynb
index 0ca95b5..46b7a3c 100644
--- a/examples/example2/example2.ipynb
+++ b/examples/example2/example2.ipynb
@@ -528,7 +528,7 @@
    "notebook_metadata_filter": "-all"
   },
   "kernelspec": {
-   "display_name": "Python 3 (ipykernel)",
+   "display_name": "Python 3",
    "language": "python",
    "name": "python3"
   },
@@ -543,11 +543,6 @@
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
    "version": "3.10.12"
-  },
-  "vscode": {
-   "interpreter": {
-    "hash": "916dbcbb3f70747c44a77c7bcd40155683ae19c65e1c03b4aa3499c5328201f1"
-   }
   }
  },
  "nbformat": 4,
diff --git a/examples/example5/example5.ipynb b/examples/example5/example5.ipynb
index 6204ac8..db75e44 100644
--- a/examples/example5/example5.ipynb
+++ b/examples/example5/example5.ipynb
@@ -414,6 +414,35 @@
     "plt.plot(tArray, pulse_vals)"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "id": "111dae29",
+   "metadata": {},
+   "source": [
+    "We could also create a space- and time-dependent parameter by reading data from an xdmf file. This is illustrated below, introducing a new parameter, \"j1pulse_fromxdmf\", which we use in an additional reaction that can be turned on or off below."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "da13025d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "j1pulse_fromxdmf = Parameter.from_xdmf(\"j1pulse_fromxdmf\", \"test_data.xdmf\", flux_unit, compartment = \"PM\")\n",
+    "add_reaction = True\n",
+    "if add_reaction:\n",
+    "        jmag = Parameter(\"jmag\", 0.0, unit.dimensionless)\n",
+    "        r1add = Reaction(\n",
+    "        \"r1add\",\n",
+    "        [],\n",
+    "        [\"B\"],\n",
+    "        param_map={\"J\": \"j1pulse_fromxdmf\", \"mag\": \"jmag\"},\n",
+    "        eqn_f_str=\"mag*J\",\n",
+    "        explicit_restriction_to_domain=\"PM\",\n",
+    "    )"
+   ]
+  },
   {
    "attachments": {},
    "cell_type": "markdown",
@@ -431,10 +460,13 @@
    "outputs": [],
    "source": [
     "pc = ParameterContainer()\n",
-    "pc.add([k4Vmax, k3r, k3f, k2f, j1pulse])\n",
-    "\n",
     "rc = ReactionContainer()\n",
-    "rc.add([r1, r2, r3, r4])"
+    "if add_reaction:\n",
+    "    pc.add([k4Vmax, k3r, k3f, k2f, j1pulse, j1pulse_fromxdmf, jmag])\n",
+    "    rc.add([r1, r2, r3, r4, r1add])\n",
+    "else:\n",
+    "    pc.add([k4Vmax, k3r, k3f, k2f, j1pulse])\n",
+    "    rc.add([r1, r2, r3, r4])"
    ]
   },
   {
@@ -599,6 +631,26 @@
     "model_cur.initialize()"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "id": "82e82d51",
+   "metadata": {},
+   "source": [
+    "If parameter was loaded from xdmf, plot here."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "901118db",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "if add_reaction:\n",
+    "    visualization.plot(model_cur.pc['j1pulse_fromxdmf'].dolfin_function,\n",
+    "                           clip_plane=(1, 1, 0), clip_origin=(0.5, 0.5, 0.5))"
+   ]
+  },
   {
    "attachments": {},
    "cell_type": "markdown",
@@ -666,7 +718,7 @@
     "    if model_cur.t >= 0.2 and not displayed:\n",
     "        visualization.plot(model_cur.sc['A'].u['u'],\n",
     "                           clip_plane=(1, 1, 0), clip_origin=(0.5, 0.5, 0.5))\n",
-    "        displayed = True\n"
+    "        displayed = True"
    ]
   },
   {
@@ -723,7 +775,7 @@
    "notebook_metadata_filter": "-all"
   },
   "kernelspec": {
-   "display_name": "Python 3 (ipykernel)",
+   "display_name": "Python 3",
    "language": "python",
    "name": "python3"
   },
@@ -738,11 +790,6 @@
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
    "version": "3.10.12"
-  },
-  "vscode": {
-   "interpreter": {
-    "hash": "916dbcbb3f70747c44a77c7bcd40155683ae19c65e1c03b4aa3499c5328201f1"
-   }
   }
  },
  "nbformat": 4,
diff --git a/smart/mesh_tools.py b/smart/mesh_tools.py
index 576a1be..e1a700a 100644
--- a/smart/mesh_tools.py
+++ b/smart/mesh_tools.py
@@ -38,7 +38,7 @@
 ]
 
 
-def implicit_curve(boundExpr):
+def implicit_curve(boundExpr, num_points=51):
     """
     Output (r,z) coordinates given a string expression in the r-z plane
     Currently the function only outputs cell coordinates in the right half
@@ -47,6 +47,8 @@ def implicit_curve(boundExpr):
     which corresponds to a substrate in most cases.
     The (r,z) coordinates start at the maximum z value for r = 0
     and end when either z = 0 or r = 0.
+    num_points can be increased for more accuracy in generating the contour,
+    but this function will run more slowly for very high num_points.
 
     Examples:
     * boundExpr = "r**2 + z**2 - 1": defines a quarter circle from (0, 1) to (1, 0)
@@ -57,13 +59,16 @@ def implicit_curve(boundExpr):
     z = sym.Symbol("z", real=True)
     outerExpr0 = outerExprRef.subs({"r": 0, "z": z})
     z0 = solveset_real(outerExpr0, z)
+    z0Array = np.array([float(z0Val) for z0Val in z0])
     rVals = [0.0]
     zVals = [float(max(z0))]
     rMax = solveset_real(outerExprRef.subs({"r": r, "z": 0.0}), r)
     if len(rMax) > 0:
-        sGap = max([float(max(z0)), float(max(rMax))]) / 50
+        sGap = max([float(max(z0)), float(max(rMax))]) / (num_points - 1)
+    elif len(z0) == 2 and np.all(z0Array > 0):  # then probably a nuclear contour
+        sGap = float(max(z0) - min(z0)) / (num_points - 1)
     else:
-        sGap = float(max(z0)) / 50
+        sGap = float(max(z0)) / (num_points - 1)
     curTan = [1, 0]
     while rVals[-1] >= 0 and zVals[-1] >= 0:
         rNext = rVals[-1] + curTan[0] * sGap
diff --git a/smart/model.py b/smart/model.py
index 5f00867..8f51475 100644
--- a/smart/model.py
+++ b/smart/model.py
@@ -19,7 +19,7 @@
 from scipy import integrate
 from sympy.utilities.lambdify import lambdify
 from pathlib import Path
-import re
+import xml.etree.ElementTree as ET
 
 try:
     from ufl_legacy.algorithms.ad import expand_derivatives
@@ -632,7 +632,8 @@ def _init_2_4_check_for_unused_parameters_species_compartments(self):
                 raise ValueError(print_str)
 
     def _init_2_5_link_compartments_to_species(self):
-        """Linking compartments and compartment dimensionality to species"""
+        """Linking compartments and compartment dimensionality to species,
+        check for consistency of diffusion coefficient definition"""
         logger.debug(
             "Linking compartments and compartment dimensionality to species",
             extra=dict(format_type="log"),
@@ -640,6 +641,12 @@ def _init_2_5_link_compartments_to_species(self):
         for species in self.sc:
             species.compartment = self.cc[species.compartment_name]
             species.dimensionality = self.cc[species.compartment_name].dimensionality
+            # convert diffusion coeff to units consistent with mesh
+            diffusion_conversion = species.diffusion_units.to(
+                species.compartment.compartment_units**2 / unit.s
+            )
+            species.diffusion_units = species.compartment.compartment_units**2 / unit.s
+            species.D *= diffusion_conversion.magnitude
 
     def _init_2_6_link_species_to_compartments(self):
         """Links species to compartments - a species is considered to be
@@ -998,6 +1005,37 @@ def _init_4_7_set_initial_conditions(self):
                     u_cur.vector().set_local(values)
                     u_cur.vector().apply("insert")
 
+        for parameter in self.pc:
+            if parameter.type == ParameterType.from_xdmf:
+                # load the time vec from xdmf file
+                tVec = self.load_timesteps_from_xdmf(parameter.xdmf_file)
+                parameter.tVec = np.array(tVec)
+                assert np.all(
+                    np.diff(parameter.tVec) >= 0.0
+                ), f"t values are not strictly increasing in {str(parameter.xdmf_file)}"
+
+                # define function space
+                if parameter.compartment not in self.cc.keys:
+                    raise ValueError(
+                        f"Compartment name {parameter.compartment} for parameter"
+                        f"{parameter.name} does not match a known compartment"
+                    )
+                V_cur = d.FunctionSpace(self.cc[parameter.compartment].dolfin_mesh, "P", 1)
+                parameter.dolfin_function = d.Function(V_cur)
+
+                vec_new = self.load_vector(parameter.h5_file, parameter.tVec)
+                vec = parameter.dolfin_function.vector()
+                mesh_map = d.dof_to_vertex_map(parameter.dolfin_function.function_space())[:]
+                if len(vec_new) != len(mesh_map):
+                    raise ValueError(
+                        f"Vector from {str(parameter.h5_file)} "
+                        f"does not match function space for {parameter.name}"
+                    )
+                else:
+                    new_vals = vec_new[mesh_map]  # reorder to match dof ordering
+                vec.set_local(new_vals)
+                vec.apply("insert")
+
     def _init_5_1_reactions_to_fluxes(self):
         """Convert reactions to flux objects"""
         logger.debug("Convert reactions to flux objects", extra=dict(format_type="log"))
@@ -1072,11 +1110,6 @@ def _init_5_2_create_variational_forms(self):
                     extra=dict(format_type="log"),
                 )
             else:
-                if Dform_units != mass_form_units:  # unit conversion for consistency
-                    diffusion_conversion = species.diffusion_units.to(
-                        species.compartment.compartment_units**2 / unit.s
-                    )
-                    D *= diffusion_conversion.magnitude
                 D_constant = d.Constant(D, name=f"D_{species.name}")
                 if self.config.flags["axisymmetric_model"]:
                     Dform = x[0] * D_constant * d.inner(d.grad(u), d.grad(v)) * dx
@@ -1825,6 +1858,19 @@ def update_time_dependent_parameters(self):
         # Update time dependent parameters
         for parameter_name, parameter in self.pc.items:
             new_value = None
+            if parameter.type == ParameterType.from_xdmf:
+                vec_new = self.load_vector(parameter.h5_file, parameter.tVec)
+                vec = parameter.dolfin_function.vector()
+                mesh_map = d.dof_to_vertex_map(parameter.dolfin_function.function_space())[:]
+                if len(vec_new) != len(mesh_map):
+                    raise ValueError(
+                        f"Vector from {str(parameter.h5_file)} "
+                        f"does not match function space for {parameter.name}"
+                    )
+                new_vals = vec_new[mesh_map]  # reorder to match dof ordering
+                vec.set_local(new_vals)
+                vec.apply("insert")
+                continue
             if not parameter.is_time_dependent:
                 continue
             if not parameter.use_preintegration:
@@ -1986,36 +2032,47 @@ def dolfin_set_function_values(self, sp, ukey, unew):
                 )
 
             # load the time vec from xdmf file
-            xdmf_file = open(xdmfCur, "r")
-            xdmf_string = xdmf_file.read()
-            found_pattern = re.findall(r"Time Value=\"?[^\s]+", xdmf_string)
-            tVec = []
-            for i in range(len(found_pattern)):
-                tVec.append(float(found_pattern[i][12:-1]))
+            tVec = self.load_timesteps_from_xdmf(xdmfCur)
+            assert np.all(
+                np.diff(np.array(tVec)) >= 0.0
+            ), f"t values are not strictly increasing in {str(xdmfCur)}"
             if self.load_init_time is None:
                 if len(tVec) == 0:
                     raise ValueError(f"Could not load time from {xdmfCur}")
                 elif len(tVec) == 1:
                     self.load_init_time = tVec[0]
                     self.load_init_idx = 0
+                    self.tvec = [tVec[0]]
                 else:
                     self.load_init_time = tVec[-2]
                     self.load_init_idx = len(tVec) - 2
+                    self.set_dt(tVec[-1] - tVec[-2])
+                    self.tvec = [tVec[-2]]
                 # set to current time
                 self.t = self.rounded_decimal(self.load_init_time)
+                if self.load_init_idx > 0:  # then store prev time
+                    self.tn = self.rounded_decimal(tVec[-3])
                 self.T.assign(self.t)
 
-            assert np.isclose(
-                tVec[self.load_init_idx], self.load_init_time
-            ), "Time value does not match load in value"
+            if tVec[self.load_init_idx] != self.load_init_time:
+                if self.load_init_time in tVec:
+                    init_idx = np.nonzero(np.array(tVec) == self.load_init_time)[0]
+                    assert len(init_idx) == 1, "t values must be unique in xdmf file"
+                    init_idx = init_idx[0]
+                else:
+                    assert np.isclose(
+                        tVec[self.load_init_idx], self.load_init_time
+                    ), "Time vector does not contain load in value"
+            else:
+                init_idx = self.load_init_idx
 
             cur_file = d.HDF5File(self.parent_mesh.mpi_comm, h5Cur, "r")
-            if not cur_file.has_dataset(f"VisualisationVector/{self.load_init_idx}"):
+            if not cur_file.has_dataset(f"VisualisationVector/{init_idx}"):
                 raise TypeError(
                     f"Unable to read initial condition for {sp.name} from file {str(unew)}"
                 )
             start_vec = d.Vector()
-            cur_file.read(start_vec, f"VisualisationVector/{self.load_init_idx}", True)
+            cur_file.read(start_vec, f"VisualisationVector/{init_idx}", True)
             cur_file.close()
             vec = self.cc[sp.compartment_name].u[ukey].vector()
             orig_vals = vec.get_local()
@@ -2199,3 +2256,49 @@ def adjust_dt(self):
         #     dt_scale = min(dt_scale * 0.8, 0.8)
         dt_cur = float(self.dt) * dt_scale
         self.set_dt(dt_cur)
+
+    def load_timesteps_from_xdmf(self, xdmffile):
+        times = []
+        tree = ET.parse(xdmffile)
+        for elem in tree.iter():
+            if elem.tag == "Time":
+                times.append(float(elem.get("Value")))
+        return times
+
+    def load_vector(self, h5_file, tVec):
+        cur_level = d.get_log_level()
+        d.set_log_level(40)  # suppress warning about rank of h5 data
+        with d.HDF5File(self.parent_mesh.mpi_comm, h5_file, "r") as cur_file:
+            # Find index associated with the starting time
+            has_timestamp = np.flatnonzero(np.isclose(tVec, float(self.t)))
+            if len(has_timestamp) > 0:
+                vec_new = d.Vector()
+                idx1 = has_timestamp[0]
+                cur_file.read(vec_new, f"VisualisationVector/{idx1}", True)
+            elif self.t > tVec[-1]:  # then starting after final time in the xdmf file
+                logger.warning(
+                    f"File {str(h5_file)} ends before current time {self.t}"
+                    "Using final time point in file instead."
+                )
+                vec_new = d.Vector()
+                idx1 = len(tVec) - 1
+                cur_file.read(vec_new, f"VisualisationVector/{idx1}", True)
+            elif self.t < tVec[0]:  # then starting before initial time in xdmf file
+                logger.warning(
+                    f"File {str(h5_file)} starts after current time {self.t}"
+                    "Using initial time point in file instead."
+                )
+                vec_new = d.Vector()
+                idx1 = 0
+                cur_file.read(vec_new, f"VisualisationVector/{idx1}", True)
+            else:  # then in between two times in the xdmf file
+                vec1 = d.Vector()
+                vec2 = d.Vector()
+                idx1 = np.flatnonzero(tVec < float(self.t))[-1]
+                idx2 = np.flatnonzero(tVec > float(self.t))[0]
+                cur_file.read(vec1, f"VisualisationVector/{idx1}", True)
+                cur_file.read(vec2, f"VisualisationVector/{idx2}", True)
+                vec_new = (vec1 + vec2) / 2
+            cur_file.close()
+        d.set_log_level(cur_level)  # restore log level
+        return vec_new
diff --git a/smart/model_assembly.py b/smart/model_assembly.py
index 73b9997..b8c6ef0 100644
--- a/smart/model_assembly.py
+++ b/smart/model_assembly.py
@@ -63,6 +63,7 @@ class ParameterType(str, Enum):
     from_file = "from_file"
     constant = "constant"
     expression = "expression"
+    from_xdmf = "from_xdmf"
 
 
 class InvalidObjectException(Exception):
@@ -565,6 +566,28 @@ class Parameter(ObjectInstance):
         use_preintegration (optional):  use preintegration in solution process if
                                      "use_preintegration" is true (defaults to false),
                                      uses sci.integrate.cumtrapz for numerical integration
+
+    To load a space-dependent parameter over time from an .xdmf file, call:
+
+    .. code:: python
+
+        param_var = Parameter.from_file(
+            name, xdmf_file, unit, compartment, group (opt),
+            notes (opt), use_preintegration (opt)
+        )
+        from_xdmf(
+        cls, name, xdmf_file, unit, compartment, group="", notes="", use_preintegration=False
+    ):
+
+    Inputs are the same as described above, except:
+
+    Args:
+        xdmf_file: name of the xdmf file with parameters saved at multiple time points
+                    (linked to an hdf5 file), saved using dolfin.XDMFFile.write()
+        compartment: string matching the compartment associated with the saved xdmf.
+                     In the current implementation, dimensions must match exactly.
+        use_preintegration: will always be set to false, not implemented for this case yet
+
     """
 
     name: str
@@ -574,6 +597,11 @@ class Parameter(ObjectInstance):
     notes: str = ""
     use_preintegration: bool = False
     sym_expr: Union[str, sym.core.Expr] = ""
+    xdmf_file: Union[str, Path] = ""
+    h5_file: Union[str, Path] = ""
+    is_time_dependent: bool = False
+    is_space_dependent: bool = False
+    compartment: str = ""
 
     def to_dict(self):
         """Convert to a dict that can be used to recreate the object."""
@@ -644,7 +672,7 @@ def from_file(cls, name, sampling_file, unit, group="", notes="", use_preintegra
         parameter.sampling_file = sampling_file
         parameter.sampling_data = sampling_data
         parameter.is_time_dependent = True
-        parameter.is_space_dependent = False  # not supported yet
+        parameter.is_space_dependent = False
         parameter.type = ParameterType.from_file
         parameter.__post_init__()
         logger.info(
@@ -654,6 +682,55 @@ def from_file(cls, name, sampling_file, unit, group="", notes="", use_preintegra
 
         return parameter
 
+    @classmethod
+    def from_xdmf(
+        cls, name, xdmf_file, unit, compartment, group="", notes="", use_preintegration=False
+    ):
+        """ "
+        Data read in from an xdmf/h5 file pairing
+        """
+        xdmf_file = Path(xdmf_file)
+        assert xdmf_file.is_file(), f"{str(xdmf_file)} could not be found to load parameter"
+
+        logger.debug(f"Loading initial condition for {name} from file")
+        if xdmf_file.suffix == ".xdmf" and xdmf_file.with_suffix(".h5").exists():
+            xdmfCur = str(xdmf_file)
+            h5Cur = xdmfCur[0:-4] + "h5"
+        else:
+            raise TypeError(f"{str(xdmf_file)} cannot be used for loading parameter, must be xdmf")
+
+        # load in xdmf file
+        logger.info(f"Loading in data for parameter {name}", extra=dict(format_type="log"))
+
+        if use_preintegration:
+            logger.warning(
+                f"Setting use_preintegration to False for parameter {name}."
+                "Not currently implemented for parameters loaded from xdmf"
+            )
+            use_preintegration = False
+        parameter = cls(
+            name,
+            0.0,
+            unit,
+            group=group,
+            notes=notes,
+            use_preintegration=use_preintegration,
+        )
+        parameter.compartment = compartment
+        # initialize instance
+        parameter.xdmf_file = xdmfCur
+        parameter.h5_file = h5Cur
+        parameter.is_time_dependent = True
+        parameter.is_space_dependent = True
+        parameter.type = ParameterType.from_xdmf
+        parameter.__post_init__()
+        logger.info(
+            f"Parameter {name} linked to xdmf file.",
+            extra=dict(format_type="log"),
+        )
+
+        return parameter
+
     @classmethod
     def from_expression(
         cls,
@@ -697,8 +774,7 @@ def from_expression(
 
         if is_time_dependent and not is_space_dependent:
             value = float(sym_expr.subs({"t": 0.0}))
-
-        if is_space_dependent:
+        else:
             dolfin_expression = d.Expression(sym.printing.ccode(sym_expr), t=0.0, degree=3)
             value = float(sym_expr.subs({"t": 0.0, "x[0]": 0.0, "x[1]": 0.0, "x[2]": 0.0}))
 
@@ -770,7 +846,9 @@ def __post_init__(self):
 
     @property
     def dolfin_quantity(self):
-        if hasattr(self, "dolfin_expression"):
+        if self.type == ParameterType.from_xdmf:
+            return self.dolfin_function * self.unit
+        elif hasattr(self, "dolfin_expression"):
             return self.dolfin_expression * self.unit
         else:
             return self.dolfin_constant * self.unit
@@ -782,7 +860,9 @@ def quantity(self):
 
     @property
     def print_val(self):
-        if self.sym_expr == "":
+        if self.type == ParameterType.from_xdmf:
+            return str(self.xdmf_file) * self.unit
+        elif self.sym_expr == "":
             self._print_val = self.value * self.unit
         else:
             self._print_val = str(self.sym_expr) * self.unit
@@ -792,7 +872,15 @@ def check_validity(self):
         """Confirm that time-dependent parameter is defined in terms of time"""
         if self.is_time_dependent:
             if all(
-                [x in ("", None) for x in [self.sampling_file, self.sym_expr, self.preint_sym_expr]]
+                [
+                    x in ("", None)
+                    for x in [
+                        self.sampling_file,
+                        self.sym_expr,
+                        self.preint_sym_expr,
+                        self.xdmf_file,
+                    ]
+                ]
             ):
                 raise ValueError(
                     f"Parameter {self.name} is marked as time dependent "
@@ -920,7 +1008,7 @@ def __post_init__(self):
         elif isinstance(self.initial_condition, Path):
             pass  # keep as path
         else:
-            raise TypeError("initial_condition must be a float or string.")
+            raise TypeError("initial_condition must be a float or string or path.")
 
         self._convert_pint_quantity_to_unit()
         self._check_input_type_validity()
@@ -965,7 +1053,10 @@ def dolfin_quantity(self):
 
     @property
     def initial_condition_quantity(self):
-        self._Initial_Concentration = self.initial_condition * self.concentration_units
+        if isinstance(self.initial_condition, Path):
+            self._Initial_Concentration = "from file"
+        else:
+            self._Initial_Concentration = self.initial_condition * self.concentration_units
         return self._Initial_Concentration
 
     @property
@@ -1576,7 +1667,7 @@ def _post_init_get_flux_units(self):
         # The expected units
         if self.is_boundary_condition:
             self._expected_flux_units = (
-                1.0 * concentration_units / compartment_units * diffusion_units
+                1.0 * (concentration_units / compartment_units) * diffusion_units
             )  # ~D*du/dn
         else:
             self._expected_flux_units = 1.0 * concentration_units / unit.s  # rhs term. ~du/dt
@@ -1669,7 +1760,7 @@ def equation_variables(self):
         variables.update({"unit_scale_factor": self.unit_scale_factor})
         free_symbols = [str(x) for x in self.equation.free_symbols]
         if "curv" in free_symbols:
-            self.curv = self.surface.curv_func * unit.dimensionless
+            self.curv = self.surface.curv_func / self.surface.compartment_units
             variables.update({"curv": self.curv})
         return variables
 


Test patch:
diff --git a/examples/example5/test_data.h5 b/examples/example5/test_data.h5
new file mode 100644
index 0000000..8b18d27
Binary files /dev/null and b/examples/example5/test_data.h5 differ
diff --git a/examples/example5/test_data.xdmf b/examples/example5/test_data.xdmf
new file mode 100644
index 0000000..ece98e0
--- /dev/null
+++ b/examples/example5/test_data.xdmf
@@ -0,0 +1,143 @@
+<?xml version="1.0"?>
+<!DOCTYPE Xdmf SYSTEM "Xdmf.dtd" []>
+<Xdmf Version="3.0" xmlns:xi="http://www.w3.org/2001/XInclude">
+  <Domain>
+    <Grid Name="TimeSeries_f_146" GridType="Collection" CollectionType="Temporal">
+      <Grid Name="mesh" GridType="Uniform">
+        <Topology NumberOfElements="512" TopologyType="Triangle" NodesPerElement="3">
+          <DataItem Dimensions="512 3" NumberType="UInt" Format="HDF">test_data.h5:/Mesh/0/mesh/topology</DataItem>
+        </Topology>
+        <Geometry GeometryType="XYZ">
+          <DataItem Dimensions="289 3" Format="HDF">test_data.h5:/Mesh/0/mesh/geometry</DataItem>
+        </Geometry>
+        <Time Value="0.000000000000000e+00" />
+        <Attribute Name="f_146" AttributeType="Scalar" Center="Node">
+          <DataItem Dimensions="289 1" Format="HDF">test_data.h5:/VisualisationVector/0</DataItem>
+        </Attribute>
+      </Grid>
+      <Grid Name="mesh" GridType="Uniform">
+        <Topology NumberOfElements="512" TopologyType="Triangle" NodesPerElement="3">
+          <DataItem Dimensions="512 3" NumberType="UInt" Format="HDF">test_data.h5:/Mesh/1/mesh/topology</DataItem>
+        </Topology>
+        <Geometry GeometryType="XYZ">
+          <DataItem Dimensions="289 3" Format="HDF">test_data.h5:/Mesh/1/mesh/geometry</DataItem>
+        </Geometry>
+        <Time Value="1.000000000000000e-01" />
+        <Attribute Name="f_146" AttributeType="Scalar" Center="Node">
+          <DataItem Dimensions="289 1" Format="HDF">test_data.h5:/VisualisationVector/1</DataItem>
+        </Attribute>
+      </Grid>
+      <Grid Name="mesh" GridType="Uniform">
+        <Topology NumberOfElements="512" TopologyType="Triangle" NodesPerElement="3">
+          <DataItem Dimensions="512 3" NumberType="UInt" Format="HDF">test_data.h5:/Mesh/2/mesh/topology</DataItem>
+        </Topology>
+        <Geometry GeometryType="XYZ">
+          <DataItem Dimensions="289 3" Format="HDF">test_data.h5:/Mesh/2/mesh/geometry</DataItem>
+        </Geometry>
+        <Time Value="2.000000000000000e-01" />
+        <Attribute Name="f_146" AttributeType="Scalar" Center="Node">
+          <DataItem Dimensions="289 1" Format="HDF">test_data.h5:/VisualisationVector/2</DataItem>
+        </Attribute>
+      </Grid>
+      <Grid Name="mesh" GridType="Uniform">
+        <Topology NumberOfElements="512" TopologyType="Triangle" NodesPerElement="3">
+          <DataItem Dimensions="512 3" NumberType="UInt" Format="HDF">test_data.h5:/Mesh/3/mesh/topology</DataItem>
+        </Topology>
+        <Geometry GeometryType="XYZ">
+          <DataItem Dimensions="289 3" Format="HDF">test_data.h5:/Mesh/3/mesh/geometry</DataItem>
+        </Geometry>
+        <Time Value="3.000000000000000e-01" />
+        <Attribute Name="f_146" AttributeType="Scalar" Center="Node">
+          <DataItem Dimensions="289 1" Format="HDF">test_data.h5:/VisualisationVector/3</DataItem>
+        </Attribute>
+      </Grid>
+      <Grid Name="mesh" GridType="Uniform">
+        <Topology NumberOfElements="512" TopologyType="Triangle" NodesPerElement="3">
+          <DataItem Dimensions="512 3" NumberType="UInt" Format="HDF">test_data.h5:/Mesh/4/mesh/topology</DataItem>
+        </Topology>
+        <Geometry GeometryType="XYZ">
+          <DataItem Dimensions="289 3" Format="HDF">test_data.h5:/Mesh/4/mesh/geometry</DataItem>
+        </Geometry>
+        <Time Value="4.000000000000000e-01" />
+        <Attribute Name="f_146" AttributeType="Scalar" Center="Node">
+          <DataItem Dimensions="289 1" Format="HDF">test_data.h5:/VisualisationVector/4</DataItem>
+        </Attribute>
+      </Grid>
+      <Grid Name="mesh" GridType="Uniform">
+        <Topology NumberOfElements="512" TopologyType="Triangle" NodesPerElement="3">
+          <DataItem Dimensions="512 3" NumberType="UInt" Format="HDF">test_data.h5:/Mesh/5/mesh/topology</DataItem>
+        </Topology>
+        <Geometry GeometryType="XYZ">
+          <DataItem Dimensions="289 3" Format="HDF">test_data.h5:/Mesh/5/mesh/geometry</DataItem>
+        </Geometry>
+        <Time Value="5.000000000000000e-01" />
+        <Attribute Name="f_146" AttributeType="Scalar" Center="Node">
+          <DataItem Dimensions="289 1" Format="HDF">test_data.h5:/VisualisationVector/5</DataItem>
+        </Attribute>
+      </Grid>
+      <Grid Name="mesh" GridType="Uniform">
+        <Topology NumberOfElements="512" TopologyType="Triangle" NodesPerElement="3">
+          <DataItem Dimensions="512 3" NumberType="UInt" Format="HDF">test_data.h5:/Mesh/6/mesh/topology</DataItem>
+        </Topology>
+        <Geometry GeometryType="XYZ">
+          <DataItem Dimensions="289 3" Format="HDF">test_data.h5:/Mesh/6/mesh/geometry</DataItem>
+        </Geometry>
+        <Time Value="6.000000000000000e-01" />
+        <Attribute Name="f_146" AttributeType="Scalar" Center="Node">
+          <DataItem Dimensions="289 1" Format="HDF">test_data.h5:/VisualisationVector/6</DataItem>
+        </Attribute>
+      </Grid>
+      <Grid Name="mesh" GridType="Uniform">
+        <Topology NumberOfElements="512" TopologyType="Triangle" NodesPerElement="3">
+          <DataItem Dimensions="512 3" NumberType="UInt" Format="HDF">test_data.h5:/Mesh/7/mesh/topology</DataItem>
+        </Topology>
+        <Geometry GeometryType="XYZ">
+          <DataItem Dimensions="289 3" Format="HDF">test_data.h5:/Mesh/7/mesh/geometry</DataItem>
+        </Geometry>
+        <Time Value="7.000000000000000e-01" />
+        <Attribute Name="f_146" AttributeType="Scalar" Center="Node">
+          <DataItem Dimensions="289 1" Format="HDF">test_data.h5:/VisualisationVector/7</DataItem>
+        </Attribute>
+      </Grid>
+      <Grid Name="mesh" GridType="Uniform">
+        <Topology NumberOfElements="512" TopologyType="Triangle" NodesPerElement="3">
+          <DataItem Dimensions="512 3" NumberType="UInt" Format="HDF">test_data.h5:/Mesh/8/mesh/topology</DataItem>
+        </Topology>
+        <Geometry GeometryType="XYZ">
+          <DataItem Dimensions="289 3" Format="HDF">test_data.h5:/Mesh/8/mesh/geometry</DataItem>
+        </Geometry>
+        <Time Value="8.000000000000000e-01" />
+        <Attribute Name="f_146" AttributeType="Scalar" Center="Node">
+          <DataItem Dimensions="289 1" Format="HDF">test_data.h5:/VisualisationVector/8</DataItem>
+        </Attribute>
+      </Grid>
+      <Grid Name="mesh" GridType="Uniform">
+        <Topology NumberOfElements="512" TopologyType="Triangle" NodesPerElement="3">
+          <DataItem Dimensions="512 3" NumberType="UInt" Format="HDF">test_data.h5:/Mesh/9/mesh/topology</DataItem>
+        </Topology>
+        <Geometry GeometryType="XYZ">
+          <DataItem Dimensions="289 3" Format="HDF">test_data.h5:/Mesh/9/mesh/geometry</DataItem>
+        </Geometry>
+        <Time Value="9.000000000000000e-01" />
+        <Attribute Name="f_146" AttributeType="Scalar" Center="Node">
+          <DataItem Dimensions="289 1" Format="HDF">test_data.h5:/VisualisationVector/9</DataItem>
+        </Attribute>
+      </Grid>
+      <Grid Name="mesh" GridType="Uniform">
+        <Topology NumberOfElements="512" TopologyType="Triangle" NodesPerElement="3">
+          <DataItem Dimensions="512 3" NumberType="UInt" Format="HDF">test_data.h5:/Mesh/10/mesh/topology</DataItem>
+        </Topology>
+        <Geometry GeometryType="XYZ">
+          <DataItem Dimensions="289 3" Format="HDF">test_data.h5:/Mesh/10/mesh/geometry</DataItem>
+        </Geometry>
+        <Time Value="1.000000000000000e+00" />
+        <Attribute Name="f_146" AttributeType="Scalar" Center="Node">
+          <DataItem Dimensions="289 1" Format="HDF">test_data.h5:/VisualisationVector/10</DataItem>
+        </Attribute>
+        <Attribute Name="f_146" AttributeType="Scalar" Center="Node">
+          <DataItem Dimensions="289 1" Format="HDF">test_data.h5:/VisualisationVector/11</DataItem>
+        </Attribute>
+      </Grid>
+    </Grid>
+  </Domain>
+</Xdmf>