Merge branch 'add-vtk-reader' into 'master'

Add reader for VTK files with vertex data

See merge request !6

.gitlab-ci.yml

@@ -64,7 +64,7 @@ test:all:
   before_script:
     - cd $CI_PROJECT_DIR/build-cmake
   script:
-    - ctest -R test --output-on-failure
+    - ctest -R test -E python --output-on-failure
 
 # Run simulations
 test:run_simulations:
@@ -76,3 +76,11 @@ test:run_simulations:
     paths:
       - build-cmake/experiments/hansbo_reproduction/
     expire_in: 1 day
+
+# Test Python code
+test:python:
+  stage: test
+  needs: []
+  script:
+    - dunecontrol --current configure
+    - dunecontrol --current bexec make test_python

docker/build-env.dockerfile

@@ -7,7 +7,7 @@ ARG DUNECI_CMAKE_FLAGS="-DDUNE_PYTHON_VIRTUALENV_SETUP=True"
 USER root
 RUN apt-get update \
     && apt-get -y upgrade \
-    && apt-get -y install gmsh libyaml-cpp-dev python3-pip python3-venv \
+    && apt-get -y install gmsh libgl1-mesa-dev libyaml-cpp-dev python3-pip python3-venv \
     && apt-get clean
 USER duneci
 
@@ -15,6 +15,22 @@ USER duneci
 WORKDIR /duneci/modules
 RUN cmake -DDUNE_PYTHON_VIRTUALENV_SETUP=True dune-common/build-cmake
 
+# BIG: Build VTK wheels
+WORKDIR /duneci/vtk
+RUN git clone https://github.com/Kitware/VTK -b release \
+    && mkdir VTK/build && cd VTK/build \
+    && cmake -DCMAKE_BUILD_TYPE=Release -DVTK_BUILD_TESTING=OFF \
+        -DVTK_BUILD_DOCUMENTATION=OFF -DVTK_BUILD_EXAMPLES=OFF -DVTK_DATA_EXCLUDE_FROM_ALL:BOOL=ON \
+        -DVTK_MODULE_ENABLE_VTK_PythonInterpreter:STRING=NO \
+        -DVTK_WHEEL_BUILD=ON -DVTK_PYTHON_VERSION=3 -DVTK_WRAP_PYTHON=ON \
+        ../ \
+    && make -j$DUNECI_PARALLEL \
+    && /duneci/modules/dune-common/build-cmake/dune-env/bin/python3 -m pip install setuptools wheel \
+    && /duneci/modules/dune-common/build-cmake/dune-env/bin/python3 setup.py bdist_wheel \
+    && /duneci/modules/dune-common/build-cmake/dune-env/bin/python3 -m pip install dist/vtk-*.whl \
+    && cd /duneci && rm -rf vtk/
+WORKDIR /duneci/modules
+
 # alugrid, testtools, pdelab
 RUN duneci-install-module --recursive https://gitlab.dune-project.org/extensions/dune-alugrid.git
 RUN duneci-install-module --recursive https://gitlab.dune-project.org/quality/dune-testtools.git

python/CMakeLists.txt

@@ -7,3 +7,10 @@ dune_python_install_package(PATH pygmsh)
 
 # Install the Python package dune.structures into the Dune virtual environment
 dune_python_install_package(PATH .)
+
+# Add Python tests
+dune_python_add_test(NAME python_tests
+    INTERPRETER ${DUNE_PYTHON_VIRTUALENV_EXECUTABLE}
+    MODULE pytest -v
+    WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/dune/structures/test
+)

python/dune/structures/__init__.py

@@ -14,3 +14,5 @@ from dune.structures.material import (
     NeoHookeanMaterial,
     StVenantKirchhoffMaterial,
 )
+
+from dune.structures.vtk import VTKVertexReader

python/dune/structures/test/test.vtu

+<?xml version="1.0"?>
+<VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">
+  <UnstructuredGrid>
+    <Piece NumberOfCells="2" NumberOfPoints="6">
+      <PointData Scalars="scalars" Vectors="vectors">
+        <DataArray type="Float32" Name="vectors" NumberOfComponents="3" format="ascii">
+          0 0 0 1 1 1 1 1 1 0 0 0 1 1 1 1 1 1
+        </DataArray>
+        <DataArray type="Float32" Name="scalars" NumberOfComponents="1" format="ascii">
+          0 1 1 0 1 1
+        </DataArray>
+      </PointData>
+      <Points>
+        <DataArray type="Float32" Name="Coordinates" NumberOfComponents="3" format="ascii">
+          0 0 0 1 0 0 0 1 0 1 1 0 0 1 0 1 0 0
+        </DataArray>
+      </Points>
+      <Cells>
+        <DataArray type="Int32" Name="connectivity" NumberOfComponents="1" format="ascii">
+          0 1 2 3 4 5
+        </DataArray>
+        <DataArray type="Int32" Name="offsets" NumberOfComponents="1" format="ascii">
+          3 6
+        </DataArray>
+        <DataArray type="UInt8" Name="types" NumberOfComponents="1" format="ascii">
+          5 5
+        </DataArray>
+      </Cells>
+    </Piece>
+  </UnstructuredGrid>
+</VTKFile>

python/dune/structures/test/test_vtk.py

+import pytest
+
+import math
+import numpy as np
+
+from dune.structures import VTKVertexReader
+
+FILE = "test.vtu"
+
+# --- Fixtures ---
+
+
+@pytest.fixture
+def vtk_reader():
+    """Create a VTKReader for test file"""
+    return VTKVertexReader(FILE)
+
+
+# --- Tests ---
+
+
+def test_init(vtk_reader):
+    """Test initialization of the reader"""
+    assert len(vtk_reader) == 2
+    assert "vectors" in vtk_reader
+    assert "scalars" in vtk_reader
+
+    with pytest.raises(KeyError):
+        vtk_reader["does-not-exist"]
+
+    with pytest.raises(FileNotFoundError):
+        VTKVertexReader("does-not-exist.vtu")
+
+
+def test_points(vtk_reader):
+    """Check if unique points are correctly detected"""
+    assert len(vtk_reader.points) == 4
+    assert vtk_reader.connectivity.shape == (2, 3)
+    assert len(vtk_reader.unique_idx) == 4
+    assert len(vtk_reader.unique_idx_inv) == 6
+
+
+def test_data(vtk_reader):
+    """Verify data"""
+    scalars = vtk_reader["scalars"]
+    assert scalars.shape == (4,)
+    assert math.isclose(float(np.mean(scalars)), 0.5)
+
+    vectors = vtk_reader["vectors"]
+    assert vectors.shape == (4, 3)
+    assert math.isclose(float(np.mean(vectors)), 0.5)

python/dune/structures/vtk.py

+import os
+from collections.abc import Mapping
+
+import numpy as np
+from vtk import vtkXMLGenericDataObjectReader, vtkUnstructuredGrid
+
+
+class VTKVertexReader(Mapping):
+    def __init__(self, filepath):
+        # Check file path
+        filepath = os.path.abspath(filepath)
+        if not os.path.isfile(filepath):
+            raise FileNotFoundError("File not found: {}".format(filepath))
+
+        # Set up the reader
+        reader = vtkXMLGenericDataObjectReader()
+        reader.SetFileName(filepath)
+        reader.Update()
+        self._grid = vtkUnstructuredGrid.SafeDownCast(reader.GetOutput())
+
+        # Retrieve the points and determine unique ones
+        self._points = np.array(
+            [self._grid.GetPoint(idx) for idx in range(self._grid.GetNumberOfPoints())]
+        )
+        self._points, self._unique_idx, self._unique_idx_inv = np.unique(
+            self._points, axis=0, return_index=True, return_inverse=True
+        )
+
+        # Get connectivity information
+        cell_array = self._grid.GetCells()
+        connect_array = cell_array.GetConnectivityArray()
+        self._connectivity = np.array(
+            [
+                connect_array.GetTuple(idx)
+                for idx in range(cell_array.GetNumberOfConnectivityIds())
+            ]
+        )[self._unique_idx_inv]
+
+        # NOTE: Assuming triangles!
+        self._connectivity = np.reshape(self._connectivity, (-1, 3))
+
+        # Prepare access to data
+        self._point_data = self._grid.GetPointData()
+
+    def __len__(self):
+        return self._point_data.GetNumberOfArrays()
+
+    def __iter__(self):
+        for idx in range(self.__len__()):
+            yield self._point_data.GetArrayName(idx)
+
+    def __getitem__(self, key):
+        if self._point_data.HasArray(key) == 1:
+            data_array = self._point_data.GetArray(key)
+            return np.squeeze(
+                np.array(
+                    [
+                        data_array.GetTuple(idx)
+                        for idx in range(self._grid.GetNumberOfPoints())
+                    ]
+                )
+            )[self._unique_idx]
+        else:
+            raise KeyError("No dataset found with name: {}".format(key))
+
+    @property
+    def points(self):
+        return self._points
+
+    @property
+    def connectivity(self):
+        return self._connectivity
+
+    @property
+    def unique_idx(self):
+        return self._unique_idx
+
+    @property
+    def unique_idx_inv(self):
+        return self._unique_idx_inv

python/setup.py

@@ -17,6 +17,9 @@ setup(
         "pyaml",
         "pygmsh",
         "meshio",
+        "vtk",
+        "numpy",
+        "pytest",
     ],
     entry_points={
         "console_scripts": [