Commit 52288d0a authored by Lukas Riedel's avatar Lukas Riedel 📝 Committed by Santiago Ospina De Los Ríos

Add new PVD file reader and make VTKReader store time

* Add a PVD file reader which stores a list of VTKReaders for each VTK
  file referenced in the PVD file. Store the times in the VTK readers.
* Add a unit test case to check the capabilities of the PVD reader.
* Expose grid bounds in the VTKReader as well as the VTKDataArray.
parent e1ff8938
......@@ -19,10 +19,9 @@
[muparser]( !131
* Coupling between transient models for water flow and solute transport !96
* Initial conditions generated from H5 input data !130
* Generic Python VTK file reader !143
* Linear interpolator for initial conditions and scaling fields !145, !156
* Parameterizations for hydrodynamic dispersion in solute transport !141
* Generic Python VTK file reader !143, !150
* Generic Python VTK file reader !143, !150, !178
* Define compile-time settings via CPP macros !144
* [Google Test]( unit test framework
as Git Submodule !159
......@@ -46,7 +46,7 @@ configure_file(manual/
-T -b html
-d ${CMAKE_CURRENT_BINARY_DIR}/_doctrees # doctree pickles dir
......@@ -10,6 +10,8 @@ This wrapper script supplies input files, and manages the execution of these
routines. The CLI is available from within the virtual environment
.. _cli-venv:
The Virtual Environment
......@@ -49,22 +49,6 @@ Package reference
Generic VTK Reader
This is a flexible VTK XML file reader which can evaluate all data sets written
by DORiE. It is derived from the abstract base class ``Mapping`` and thus
supports accessing the data arrays stored inside like a Python ``dict``. Use
the data array names as keys.
.. autoclass:: dorie.utilities.vtktools.vtkreader.VTKReader
The ``VTKReader`` returns instances of ``VTKDataArray``, which is a wrapper
around the underlying VTK data array.
.. autoclass:: dorie.utilities.vtktools.vtkreader.VTKDataArray
Generic VTK File Readers
For an interactive inspection and analysis of the VTK files written by DORiE we
recommend using ParaView. For more information, learn
:doc:`how to use ParaView </cookbook/2-paraview/tutorial>`. More detailed analyses
typically require more complex operations on the data which is likely easier
in Python. The DORiE Python module provides facilities based on the
`VTK Python module <>`_.
.. contents::
:depth: 2
For using the modules and classes listed below, users need to install a set of
required software packages. If you followed the installation instructions in
the :doc:`ReadMe </introduction/readme>`, you are already set and can start
using them inside the :ref:`virtual environment <cli-venv>`.
If you install the DORiE Python package separately, you may need to install the
required VTK library, depending on your operating system. **Ubuntu** users need
to install ``python-vtk`` via
.. code-block:: console
apt-get install python-vtk
Example Code
This example code loads the resulting PVD file of a DORiE simulation, reports
the simulation times for which files are available and prints the water content
at the center of the domain for every file.
.. code-block:: python
from dorie.utilities.vtktools.vtkreader import PVDReader
# Load the PVD file written by a DORiE simulation run
pvd_reader = PVDReader("result.pvd")
# Report the times
# Iterate over the output files
# NOTE: `dataset` is an instance of `VTKReader`
for dataset in pvd_reader:
# Retrieve time from file directly
time = dataset.time
# Fetch the water content data array
# NOTE: `water_content` is an instance of `VTKDataArray`
water_content = dataset["theta_w"]
# Compute the center of the domain
center = np.mean(water_content.bounds, axis=1)
# Report the water content at this location
Source Code Documentation
This is a flexible VTK XML file reader which can evaluate all data sets written
by DORiE. It is derived from the abstract base class ``Mapping`` and thus
supports accessing the data arrays stored inside like a Python ``dict``. Use
the data array names as keys.
.. note:: VTK files do not store time inherently. When loading a file directly
using the ``VTKReader``, the ``time`` attribute will be set to
``None``. Use the
:class:`dorie.utilities.vtktools.vtkreader.PVDReader` to load all VTK
files of a simulation run and access their times.
.. autoclass:: dorie.utilities.vtktools.vtkreader.VTKReader
The ``VTKReader`` returns instances of ``VTKDataArray``, which is a wrapper
around the underlying VTK data array.
.. autoclass:: dorie.utilities.vtktools.vtkreader.VTKDataArray
DORiE writes PVD files containing references to the VTK file(s) stored for
every time step, and the associated simulation time. The ``PVDReader`` stores
a sequence of :class:`dorie.utilities.vtktools.vtkreader.VTKReader` for these
files and orders them by time.
.. autoclass:: dorie.utilities.vtktools.vtkreader.PVDReader
import pytest
import numpy as np
from ..utilities.vtktools.vtkreader import VTKReader, VTKDataArray
from ..utilities.vtktools.vtkreader import VTKReader, VTKDataArray, PVDReader
VTK_CELL = "vtkreader-cell.vtu"
VTK_VERTEX = "vtkreader-vertex.vtu"
PVD_FILE = "vtkreader.pvd"
ARRAYS = ["head", "flux", "K_0", "Theta", "theta_w"]
......@@ -19,6 +20,11 @@ def vtk_vertex ():
"""Create a VTKReader for the vertex VTK test file"""
return VTKReader(VTK_VERTEX)
def pvd_reader():
return PVDReader(PVD_FILE)
# Tests ------------------------------------------------------------------
def test_init ():
......@@ -53,6 +59,7 @@ def test_array_properties (vtk_cell):
assert np.allclose(array.bounds, [[0, 1],
[0, 1],
[0, 0]])
assert np.allclose(vtk_cell.bounds, array.bounds)
assert vtk_cell["flux"].number_of_components == 3
......@@ -87,3 +94,19 @@ def test_evaluate (vtk_cell, vtk_vertex):
pos = pos + [1.0, 0.0, 1.0]
with pytest.raises(ValueError):
def test_pvd_reader(pvd_reader, vtk_cell, vtk_vertex):
# Check sorted times via property
times_real = [0.0, 1.0, 1.2]
times_read = list(pvd_reader.times)
assert times_read == times_real
# Check sorted times via access of the VTKReaders
for vtk, time in zip(pvd_reader, times_real):
assert vtk.time == time
# Check if the VTKReaders point to the same files
assert pvd_reader[0] == vtk_cell
assert pvd_reader[1] == vtk_cell
assert pvd_reader[2] == vtk_vertex
<?xml version="1.0"?>
<VTKFile type="Collection" version="0.1" byte_order="LittleEndian">
<DataSet timestep="1.2" group="" part="0" name="" file="vtkreader-vertex.vtu"/>
<DataSet timestep="0" group="" part="0" name="" file="vtkreader-cell.vtu"/>
<DataSet timestep="1e0" group="" part="0" name="" file="vtkreader-cell.vtu"/>
from import Mapping
import os
import math
import xml.etree.ElementTree as ET
from import Mapping, Sequence
from vtk import vtkXMLGenericDataObjectReader, \
vtkCellTreeLocator, \
vtkDataSet, \
......@@ -7,7 +11,12 @@ from vtk import vtkXMLGenericDataObjectReader, \
import numpy as np
class VTKDataArray:
"""Wrapper for evaluating a VTK data array"""
"""Wrapper for evaluating a VTK data array
bounds (:obj:`np.array`, shape ``(3, 2)``): Outer bounds of the grid
this data array lives on.
def __init__(self, dataset, locator, data_array, data_type):
"""Initialize this wrapper from a dataset, a cell locator,
......@@ -37,22 +46,12 @@ class VTKDataArray:
# store physical dataset bounds
self._bounds = np.reshape(dataset.GetBounds(), (3, 2))
self.bounds = np.reshape(dataset.GetBounds(), (3, 2))
# cache for cell and ID
self._cell = vtkGenericCell()
self._cell_id = -1 # invalid cell
def bounds (self):
"""Return the outer bounds of the grid this data array lives on.
np.array: The rectangular bounds in any spatial direction.
Shape: ``(3, 2)``.
return self._bounds
def number_of_components (self):
"""Number of components of the evaluated quantity."""
......@@ -141,23 +140,32 @@ class VTKReader (Mapping):
vertex data set in both cases.
file_name (str): Name of the VTK file to open. Any file in XML
file_name (str): Name of the VTK file to open. Any VTK file in XML
format is supported.
filepath (os.path): The absolute path to the VTK file. Used for
comparing two instances of ``VTKReader``.
time (float): The time associated with the VTK file. This will be
``None`` if the file is loaded directly and not from a
PVD file, see :py:class:`PVDReader`.
_reader: Generic VTK XML data file reader.
_dataset: The complete data set loaded by the reader.
_locator: Fast locator of grid cells based on global coordinates.
_point_data: Pointer to the VTK point data arrays.
_cell_data: Pointer to the VTK cell data arrays.
time = None
filepath = None
def __init__(self, file_name):
"""Open the target VTK file, and build a cell locator."""
self.filepath = os.path.abspath(file_name)
if not os.path.isfile(self.filepath):
raise RuntimeError("File not found: {}".format(self.filepath))
# set up the reader
self._reader = vtkXMLGenericDataObjectReader()
self._dataset = vtkDataSet.SafeDownCast(self._reader.GetOutput())
......@@ -170,6 +178,24 @@ class VTKReader (Mapping):
self._point_data = self._dataset.GetPointData()
self._cell_data = self._dataset.GetCellData()
def bounds(self):
""":obj:`np.array`, shape ``(3, 2)``: The rectangular bounds of the
dataset grid in any spatial direction. This value is shared among all
:py:class:`VTKDataArray` objects retrieved from this instance.
return np.reshape(self._dataset.GetBounds(), (3, 2))
def __eq__(self, other):
"""Compare to another VTKReader instance based on the VTK file path.
This does not check the stored times, meaning that two instances
might compare equal in terms of the file path but may store
different times.
return self.filepath == other.filepath
def __len__(self):
return self._point_data.GetNumberOfArrays() \
+ self._cell_data.GetNumberOfArrays()
......@@ -221,3 +247,66 @@ class VTKReader (Mapping):
data = self._cell_data
return data.GetArrayName(array_id)
class PVDReader (Sequence):
"""A PVD file reader that stores one instance of :py:class:`VTKReader` for
every referenced VTK file. In addition, the respective times are stored
inside the readers and can be accessed independently through this instance.
Use this object like a list of :py:class:`VTKReader`.
The VTK files are guaranteed to be ordered with respect to time, regardless
of their actual order in the PVD file.
file_path (str): Path to the PVD file to open.
_datasets (list): :py:class:`VTKReader` instances.
def __init__(self, file_path):
# Open the PVD file
file_path = os.path.abspath(file_path)
tree = ET.parse(file_path)
root = tree.getroot()
dirname = os.path.dirname(file_path)
# Check structure
if not root.get("type") == "Collection" \
or root.get("version") != "0.1":
raise IOError("Input PVD file must have type 'Collection' and "
"version 0.1")
collection = root.find("Collection")
# Load the files
self._datasets = []
for vtk_file in collection.findall("DataSet"):
# PVD file contains local file paths
vtk_file_path = os.path.join(dirname, vtk_file.get("file"))
vtkreader = VTKReader(vtk_file_path)
# Extract time
time = float(vtk_file.get("timestep"))
vtkreader.time = time
# Store the readers
# Sort the elements according to the time
self._datasets.sort(key=lambda file: file.time)
def __getitem__(self, index):
return self._datasets[index]
def __len__(self):
return len(self._datasets)
def times(self):
""":obj:`iterable` of :obj:`float`: The times associated with the
stored VTK files.
for dataset in self._datasets:
yield dataset.time
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment