Commit b0290153 authored by Lukas Riedel's avatar Lukas Riedel

Merge branch 'update-vtkreader' into 'master'

Implement new VTK file reader

See merge request !143
parents ebbc64b0 305a9e28
......@@ -10,7 +10,7 @@ variables:
BASE_IMAGE: dorie/dune-env
# Use semantic versioning (not the version of DUNE) and bump according to
# to whether changes are backwards-compatible or not.
IMAGE_VERSION: "1.0"
IMAGE_VERSION: "1.1"
DUNE_ENV_IMAGE: ${BASE_IMAGE}:img-v${IMAGE_VERSION}
CMAKE_FLAGS:
......@@ -184,6 +184,13 @@ test:unit-tests:
- $CI_PROJECT_DIR/build-cmake/dune/dorie/test
expire_in: 1 day
test:python-tests:
<<: *test
dependencies: []
script:
- $DUNECONTROL --only=dorie configure
- $DUNECONTROL --only=dorie make test_python
# --- Deploy jobs ---
deploy:dockerhub-devel: &deploy
stage: deploy
......
......@@ -39,6 +39,7 @@
[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
### Changed
* `Simulation` is renamed `RichardsSimulation` and moved to
......
......@@ -78,6 +78,7 @@ by CI tests.
| [yaml-cpp](https://github.com/jbeder/yaml-cpp) | >= 5.2.0 |
| [spdlog](https://github.com/gabime/spdlog) | 1.1.0 | Included as Git Submodule
| [muparser](http://beltoforion.de/article.php?a=muparser) | master |
| [VTK](https://vtk.org/) | >= 7.1.1 | For the Python module only
| [dune-common](https://gitlab.dune-project.org/core/dune-common) | releases/2.6
| [dune-geometry](https://gitlab.dune-project.org/core/dune-geometry) | releases/2.6
| [dune-grid](https://gitlab.dune-project.org/core/dune-grid) | releases/2.6
......@@ -95,8 +96,8 @@ by CI tests.
| ---------| -------------- | -------- |
| [dune-testtools](https://gitlab.dune-project.org/quality/dune-testtools) | releases/2.6 | Handles system tests
| [doxygen](http://www.stack.nl/~dimitri/doxygen/) | 1.8.13 | Builds documentation
| [METIS](http://glaros.dtc.umn.edu/gkhome/views/metis) | 5 | For parallel runs
| [ParMETIS](http://glaros.dtc.umn.edu/gkhome/views/metis) | 4 | For parallel runs
| [METIS](http://glaros.dtc.umn.edu/gkhome/views/metis) | 5 | For parallel runs
| [ParMETIS](http://glaros.dtc.umn.edu/gkhome/views/metis) | 4 | For parallel runs
#### Step-by-step Instructions
......@@ -119,7 +120,7 @@ If you installed [Anaconda](https://conda.io/docs/user-guide/install/download.ht
git libatlas-base-dev libfftw3-dev libfftw3-mpi-dev \
libfreetype6-dev libhdf5-mpi-dev libmuparser-dev \
libopenmpi-dev libpng-dev libsuperlu-dev libyaml-cpp-dev \
libxft-dev python3-dev python3-pip
libxft-dev python3-dev python3-pip python3-vtk7
**macOS:**
......@@ -302,12 +303,14 @@ or former versions of itself. This ensures that DORiE is running correctly and
producing the expected results. We distinguish _unit tests_ for testing certain
features of the code, and _system tests_ for verifying the results of the
final application. As system tests require executing the DUNE solvers,
it is recommended to build them in a `Release` environment.
it is recommended to build them in a `Release` environment. Additionaly, there
is a set of tests for the Python module.
| Test category | Build tests | Execute tests | Recommended build type |
| ------------- | ----------- | ------------- | ---------------------- |
| Unit tests | `make build_unit_tests` | `make unit_tests` | `Debug`
| System tests | `make build_system_tests` | `make system_tests` | `Release`
| Python tests | N/A | `make test_python` | any
The `make` commands are to be executed from within the `build-cmake` directory.
......
......@@ -49,6 +49,17 @@ 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.
.. autoclass:: dorie.utilities.vtktools.vtkreader.VTKReader
:members:
.. automethod:: __init__
dorie.utilities.check_path
--------------------------
......
......@@ -6,30 +6,35 @@ ARG PROCNUM=1
ARG CC=gcc
ARG CXX=g++
RUN apt-get clean && apt-get update && apt-get install -y \
clang \
cmake \
doxygen \
gcc \
g++ \
gfortran \
git \
libatlas-base-dev \
libfftw3-dev \
libfftw3-mpi-dev \
libfreetype6-dev \
libhdf5-mpi-dev \
libmetis-dev \
libmuparser-dev \
libopenmpi-dev \
libpng-dev \
libparmetis-dev \
libsuperlu-dev \
libxft-dev \
libyaml-cpp-dev \
locales \
python3-dev \
python3-pip \
# disable any prompts while installing packages
ENV DEBIAN_FRONTEND=noninteractive
RUN apt-get clean \
&& apt-get update \
&& apt-get install -y \
clang \
cmake \
doxygen \
gcc \
g++ \
gfortran \
git \
libatlas-base-dev \
libfftw3-dev \
libfftw3-mpi-dev \
libfreetype6-dev \
libhdf5-mpi-dev \
libmetis-dev \
libmuparser-dev \
libopenmpi-dev \
libpng-dev \
libparmetis-dev \
libsuperlu-dev \
libxft-dev \
libyaml-cpp-dev \
locales \
python3-dev \
python3-pip \
python3-vtk7 \
&& apt-get clean
RUN rm -rf /var/lib/apt/lists/* \
......
......@@ -8,3 +8,12 @@ dune_python_install_package(PATH dorie)
# set cache variable
set(DORIE_PYTHON_MODULES ${python_paths} CACHE PATH "Paths to the python modules")
# add Python tests
dune_python_add_test(NAME python_tests
COMMAND ${CMAKE_BINARY_DIR}/run-in-dune-env
python -m pytest -v
${CMAKE_CURRENT_SOURCE_DIR}/dorie/dorie/test/
WORKING_DIRECTORY
${CMAKE_CURRENT_SOURCE_DIR}/dorie/dorie/test/
)
import pytest
import numpy as np
from dorie.utilities.vtktools.vtkreader import VTKReader
VTK_CELL = "vtkreader-cell.vtu"
VTK_VERTEX = "vtkreader-vertex.vtu"
DEFAULT_ARRAY = "head"
ARRAYS = ["head", "flux", "K_0", "Theta", "theta_w"]
# Fixtures ---------------------------------------------------------------
@pytest.fixture
def vtk_cell ():
"""Create a VTKReader for the cell VTK test file"""
return VTKReader(VTK_CELL, DEFAULT_ARRAY)
@pytest.fixture
def vtk_vertex ():
"""Create a VTKReader for the vertex VTK test file"""
return VTKReader(VTK_VERTEX, DEFAULT_ARRAY)
# Tests ------------------------------------------------------------------
def test_init ():
"""Test if initialization works"""
VTKReader(VTK_CELL, DEFAULT_ARRAY)
VTKReader(VTK_VERTEX, DEFAULT_ARRAY)
def test_array_set (vtk_cell, vtk_vertex):
"""Test if all included arrays are found"""
for array in ARRAYS:
vtk_cell.set_data_array(array)
vtk_vertex.set_data_array(array)
with pytest.raises(RuntimeError):
vtk_cell.set_data_array("other")
def test_properties (vtk_cell):
"""Test the VTKReader properties"""
assert np.allclose(vtk_cell.bounds, [[0, 1],
[0, 1],
[0, 0]])
assert vtk_cell.number_of_components == 1
vtk_cell.set_data_array("flux")
assert vtk_cell.number_of_components == 3
def test_evaluate (vtk_cell, vtk_vertex):
"""Test the evaluation function"""
center = np.mean(vtk_cell.bounds, axis=1)
# should return scalar
with pytest.raises(TypeError):
len(vtk_cell.evaluate(center))
# check return of vector
vtk_vertex.set_data_array("flux")
assert len(vtk_vertex.evaluate(center)) == 3
# check values are equal
for array in ["head", "K_0", "flux"]:
vtk_cell.set_data_array(array)
vtk_vertex.set_data_array(array)
assert np.allclose(vtk_cell.evaluate(center),
vtk_vertex.evaluate(center))
# check evaluation errors
# position has wrong shape
with pytest.raises(ValueError):
vtk_vertex.evaluate([0.0, 0.0])
# position out of bounds
pos = vtk_vertex.bounds[:, 1]
pos = pos + [1.0, 0.0, 1.0]
with pytest.raises(ValueError):
vtk_vertex.evaluate(pos)
<?xml version="1.0"?>
<VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">
<UnstructuredGrid>
<Piece NumberOfCells="1" NumberOfPoints="4">
<CellData Scalars="head" Vectors="flux">
<DataArray type="Float32" Name="head" NumberOfComponents="1" format="binary">
BAAAAA==AAAAvw==
</DataArray>
<DataArray type="Float32" Name="flux" NumberOfComponents="3" format="binary">
DAAAAA==AAAAgAAAAIAAAAAA
</DataArray>
<DataArray type="Float32" Name="K_0" NumberOfComponents="1" format="binary">
BAAAAA==pIy4N4==
</DataArray>
<DataArray type="Float32" Name="theta_w" NumberOfComponents="1" format="binary">
BAAAAA==Jh4iPh==
</DataArray>
<DataArray type="Float32" Name="Theta" NumberOfComponents="1" format="binary">
BAAAAA==jKPqPq==
</DataArray>
</CellData>
<Points>
<DataArray type="Float32" Name="Coordinates" NumberOfComponents="3" format="binary">
MAAAAA==AAAAAAAAAAAAAAAAAACAPwAAAAAAAAAAAAAAAAAAgD8AAAAAAACAPwAAgD8AAAAA
</DataArray>
</Points>
<Cells>
<DataArray type="Int32" Name="connectivity" NumberOfComponents="1" format="binary">
EAAAAA==AAAAAAEAAAADAAAAAgAAAA==
</DataArray>
<DataArray type="Int32" Name="offsets" NumberOfComponents="1" format="binary">
BAAAAA==BAAAAA==
</DataArray>
<DataArray type="UInt8" Name="types" NumberOfComponents="1" format="binary">
AQAAAA==CQ==
</DataArray>
</Cells>
</Piece>
</UnstructuredGrid>
</VTKFile>
<?xml version="1.0"?>
<VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">
<UnstructuredGrid>
<Piece NumberOfCells="1" NumberOfPoints="4">
<PointData Scalars="head" Vectors="flux">
<DataArray type="Float32" Name="head" NumberOfComponents="1" format="binary">
EAAAAA==AAAAAAAAAAAAAIC/AACAvw==
</DataArray>
<DataArray type="Float32" Name="flux" NumberOfComponents="3" format="binary">
MAAAAA==AAAAgAAAAIAAAAAAAAAAgAAAAIAAAAAAAAAAgAAAAIAAAAAAAAAAgAAAAIAAAAAA
</DataArray>
<DataArray type="Float32" Name="K_0" NumberOfComponents="1" format="binary">
EAAAAA==pIy4N6SMuDekjLg3pIy4N4==
</DataArray>
<DataArray type="Float32" Name="theta_w" NumberOfComponents="1" format="binary">
EAAAAA==UriePlK4nj6N0Uo9jdFKPd==
</DataArray>
<DataArray type="Float32" Name="Theta" NumberOfComponents="1" format="binary">
EAAAAA==AACAPwAAgD9Xv449V7+OPb==
</DataArray>
</PointData>
<Points>
<DataArray type="Float32" Name="Coordinates" NumberOfComponents="3" format="binary">
MAAAAA==AAAAAAAAAAAAAAAAAACAPwAAAAAAAAAAAAAAAAAAgD8AAAAAAACAPwAAgD8AAAAA
</DataArray>
</Points>
<Cells>
<DataArray type="Int32" Name="connectivity" NumberOfComponents="1" format="binary">
EAAAAA==AAAAAAEAAAADAAAAAgAAAA==
</DataArray>
<DataArray type="Int32" Name="offsets" NumberOfComponents="1" format="binary">
BAAAAA==BAAAAA==
</DataArray>
<DataArray type="UInt8" Name="types" NumberOfComponents="1" format="binary">
AQAAAA==CQ==
</DataArray>
</Cells>
</Piece>
</UnstructuredGrid>
</VTKFile>
from vtk import vtkXMLGenericDataObjectReader, \
vtkCellTreeLocator, \
vtkDataSet, \
vtkGenericCell, \
reference
import numpy as np
class VTKReader:
"""A generic VTK reader for any type of VTK XML files.
It is intended to evaluate cell or vertex (point) data arrays on any point
of the grid.
Attributes:
_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.
_cell: Generic grid cell implementation.
_cell_id: ID of the cell on the grid.
_eval: Evaluation method based on the selected data array type.
_data_array: The actual data to evaluate.
"""
def __init__(self, file_name, array_name, data_type=None):
"""Open the target VTK file, build a cell locator and open the data array.
See :func:`~set_data_array` for details on the arguments regarding the
data array.
Args:
file_name (str): Path to the VTK file.
array_name (str): Name of the data array to evaluate.
data_type (str): Type of the data array.
"""
self._eval = None
self._data_array = None
# set up the reader
self._reader = vtkXMLGenericDataObjectReader()
self._reader.SetFileName(file_name)
self._reader.Update()
self._dataset = vtkDataSet.SafeDownCast(self._reader.GetOutput())
# build the locator
self._locator = vtkCellTreeLocator()
self._locator.SetDataSet(self._dataset)
self._locator.BuildLocator()
# cache for cell and ID
self._cell = vtkGenericCell()
self._cell_id = -1 # invalid cell
# set the evaluation data array
self.set_data_array(array_name, data_type)
def set_data_array (self, array_name, data_type=None):
"""Switch the data array (dataset) to evaluate.
Args:
array_name (str): The name of the data array to evaluate.
data_type (str, optional): The type of the data array (either
``cell`` or ``vertex``). If this is ``None``, vertex arrays
take precedence.
Raises:
RuntimeError: Data array was not found.
"""
self._data_array = None
def retrieve_array (data, name):
"""Return the array if it exists in the data."""
if data.HasArray(name) == 1:
return data.GetArray(name)
else:
return None
if data_type == "cell" or data_type is None:
# retrieve cell data array...
self._data_array = retrieve_array(self._dataset.GetCellData(),
array_name)
self._eval = self.__eval_cell__
if self._data_array is None \
and (data_type == "vertex" or data_type is None):
# ...or vertex data array
self._data_array = retrieve_array(self._dataset.GetPointData(),
array_name)
self._eval = self.__eval_vertex__
if self._data_array is None:
if data_type is None:
raise RuntimeError("Data array '{}' not found in VTK file"
.format(array_name))
else:
raise RuntimeError("Data array '{}' of type '{}' not found in "
"VTK file. Supported types are: 'cell', "
"'vertex'".format(array_name, data_type))
self._num_comp = self._data_array.GetNumberOfComponents()
@property
def bounds (self):
"""Return the outer bounds of the grid.
Returns:
np.array: The rectangular bounds in any spatial direction.
Shape: ``(3, 2)``.
"""
bounds = self._dataset.GetBounds()
return np.reshape(bounds, (3, 2))
@property
def number_of_components (self):
"""Number of components of the evaluated quantity."""
return self._num_comp
def evaluate (self, position):
"""Evaluate the selected data array at a certain position.
Args:
position (list): 3D coordinate probing the data array.
Returns:
np.array: The value of the data array at the given position.
The output is squeezed: Scalar data arrays return scalar values,
vector arrays return numpy arrays with shape ``(3)``.
Raises:
ValueError: ``position`` argument not a 3D list or array.
ValueError: Position outside grid.
RuntimeError: Numerical error in cell locator algorithm.
"""
if len(position) != 3:
raise ValueError("Position must be a 3D lest/array")
# unused stubs
zeros = [0, 0, 0]
zero_ref = reference(0)
# interpolation weights
weights = np.zeros(10)
# check if we cached the correct cell
on_cell = self._cell.EvaluatePosition(position,
zeros,
zero_ref,
zeros,
zero_ref,
weights)
if on_cell == 0 or self._cell_id < 0:
# locate correct cell
self._cell_id = self._locator.FindCell(position,
0,
self._cell,
zeros,
weights)
if self._cell_id < 0:
raise ValueError("Position outside domain: {}".format(position))
elif on_cell < 0:
raise RuntimeError("Error while evaluating if position {} is on "
"cell {}".format(position, self._cell_id))
# evaluate the data
data = self._eval(weights)
return np.squeeze(data)
def __eval_cell__ (self, weights):
"""Evaluation function for a cell-type data array."""
return np.asarray(self._data_array.GetTuple(self._cell_id))
def __eval_vertex__ (self, weights):
"""Evaluation function for a vertex-type data array."""
data = np.zeros(self._num_comp)
for point in range(self._cell.GetNumberOfPoints()):
point_id = self._cell.GetPointId(point)
tup = np.asarray(self._data_array.GetTuple(point_id))
data += weights[point] * tup
return data
......@@ -9,7 +9,8 @@ setup(name='dorie',
author='Lukas Riedel',
author_email='dorieteam@iup.uni-heidelberg.de',
url='https://ts-gitlab.iup.uni-heidelberg.de/dorie/dorie',
packages=find_packages(),
packages=find_packages(exclude=["dorie.test"]),
test_suite='py.test',
install_requires= ['cycler',
'configparser',
'h5py',
......@@ -18,11 +19,13 @@ setup(name='dorie',
'Pillow',
'pygmsh',
'pyparsing',
'pytest',
'python-dateutil',
'pytz',
'pyyaml',
'scipy',
'six',
'vtk',
'wheel'
],
scripts=[
......
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