Commit cedc75d3 authored by Santiago Ospina De Los Ríos's avatar Santiago Ospina De Los Ríos

Merge branch 'new-pvd-reader' into 'master'

Add new PVD Reader

See merge request !178
parents e1ff8938 52288d0a
......@@ -19,10 +19,9 @@
[muparser](http://beltoforion.de/article.php?a=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](https://github.com/google/googletest) unit test framework
as Git Submodule !159
......
......@@ -46,7 +46,7 @@ configure_file(manual/config-file.rst.in
configure_file(conf.py.in conf.py)
add_custom_target(sphinx_html
COMMAND ${SPHINX_EXECUTABLE}
COMMAND ${CMAKE_BINARY_DIR}/run-in-dune-env ${SPHINX_EXECUTABLE}
-T -b html
-c ${CMAKE_CURRENT_BINARY_DIR} # conf.py dir
-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
(``venv``).
.. _cli-venv:
The Virtual Environment
=======================
......
......@@ -49,22 +49,6 @@ Package reference
:members:
:undoc-members:
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
:members:
dorie.utilities.check_path
--------------------------
......
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 <https://pypi.org/project/vtk/>`_.
.. contents::
:depth: 2
:local:
Prerequisites
=============
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
print(pvd_reader.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
print(water_content.evaluate(center))
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
:members:
:inherited-members:
The ``VTKReader`` returns instances of ``VTKDataArray``, which is a wrapper
around the underlying VTK data array.
.. autoclass:: dorie.utilities.vtktools.vtkreader.VTKDataArray
:members:
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
:members:
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"
DEFAULT_ARRAY = "head"
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)
@pytest.fixture
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):
array_vertex.evaluate(pos)
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">
<Collection>
<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"/>
</Collection>
</VTKFile>
from collections.abc import Mapping
import os
import math
import xml.etree.ElementTree as ET
from collections.abc 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
Attributes:
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
@property
def bounds (self):
"""Return the outer bounds of the grid this data array lives on.
Returns:
np.array: The rectangular bounds in any spatial direction.
Shape: ``(3, 2)``.
"""
return self._bounds
@property
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.
Args:
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.
Attributes:
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._reader.SetFileName(file_name)
self._reader.SetFileName(self.filepath)
self._reader.Update()
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()
@property
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.
Note:
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.
Args:
file_path (str): Path to the PVD file to open.
Attributes:
_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
self._datasets.append(vtkreader)
# 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)
@property
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