Merge branch 'expand_mesh_diff' into 'main'

Expand mesh diff to accept lists or array of meshes as input

See merge request !128

docs/examples/howto_meshlib/plot_calculate_diff.py

+"""
+Differences between meshes
+==========================
+
+.. sectionauthor:: Feliks Kiszkurno (Helmholtz Centre for Environmental Research GmbH - UFZ)
+
+This example explains how to use functions from meshlib to calculate differences between meshes.
+"""
+
+# %%
+
+# sphinx_gallery_start_ignore
+
+# sphinx_gallery_thumbnail_path = '_static/examples/meshlib/plot_calculate_diff_thumbnail.png'
+
+# sphinx_gallery_end_ignore
+
+from ogstools.meshlib import difference, difference_matrix, difference_pairwise
+from ogstools.meshplotlib.examples import meshseries_THM_2D
+from ogstools.propertylib import presets
+
+mesh_property = presets.temperature
+
+# %%
+# 0. Introduction
+# ---------------
+# This example shows how to calculate differences between meshes.
+# It is possible to calculate the difference between multiple meshes at the same time. Multiple meshes can be provided either as list or Numpy arrays.
+# 4 ways of calculating the difference are presented here.
+
+# %%
+# 1. Difference between two meshes
+# --------------------------------
+# The simplest case is calculating the difference between two meshes. For this example, we read two different timesteps from a MeshSeries. It is not required that they belong to the same MeshSeries object. As long, as the meshes share the same topology and contain the mesh_property of interest, the difference will work fine.
+mesh1 = meshseries_THM_2D.read(0)
+mesh2 = meshseries_THM_2D.read(-1)
+
+# %% [markdown]
+# The following call will return a mesh containing the difference of the mesh_property between the two provided meshes:
+#
+# .. math::
+#
+#   \text{Mesh}_1 - \text{Mesh}_2
+#
+
+mesh_diff = difference(mesh_property, mesh1, mesh2)
+
+# %% [markdown]
+# This returned object will be a PyVista UnstructuredGrid object:
+
+# %%
+print(f"Type of mesh_diff: {type(mesh_diff)}")
+
+# %% [markdown]
+# 2. Pairwise difference
+# ----------------------
+# In order to calculate pairwise differences, two lists or arrays of equal length have to be provided as the input. They can contain an arbitrary number of different meshes, as long as their length is equal.
+#
+# Consider the two following lists:
+#
+# .. math::
+#
+#   \text{List}_{1} = \begin{bmatrix} A_{1} & B_{1} & C_{1}\\ \end{bmatrix}
+#
+# and
+#
+# .. math::
+#
+#   \text{List}_{2} = \begin{bmatrix} A_{2} & B_{2} & C_{2}\\ \end{bmatrix}
+#
+# The output array will contain following differences between meshes:
+#
+# .. math::
+#
+#   \begin{bmatrix} A_{1}-A_{2} & B_{1}-B_{2} & C_{1}-C_{2}\\ \end{bmatrix}
+#
+# and will have the same shape as the input lists. As in the example below:
+
+meshes_1 = [mesh1] * 3
+meshes_2 = [mesh2] * 3
+
+mesh_diff_pair_wise = difference_pairwise(mesh_property, meshes_1, meshes_2)
+
+# %%
+print(f"Length of mesh_list1: {len(meshes_1)}")
+
+# %%
+print(f"Length of mesh_list2: {len(meshes_2)}")
+
+# %%
+print(f"Shape of mesh_diff_pair_wise: {mesh_diff_pair_wise.shape}")
+
+# %% [markdown]
+# 3. Matrix difference - one array
+# --------------------------------
+# If only one list or array is provided, the differences between every possible pair of objects in the array will be returned.
+#
+# Consider following list:
+#
+# .. math::
+#
+#   \text{List} = \begin{bmatrix} A & B & C\\ \end{bmatrix}
+#
+# The output array will contain following differences between meshes:
+#
+# .. math::
+#
+#   \begin{bmatrix} A-A & B-A & C-A\\ A-B & B-B & C-B \\ A-C & B-C & C-C \\ \end{bmatrix}
+#
+# and will have shape of (len(mesh_list), len(mesh_list)). As in the following example:
+
+mesh_list = [mesh1, mesh2, mesh1, mesh2]
+
+mesh_diff_matrix = difference_matrix(mesh_property, mesh_list)
+
+# %%
+print(f"Length of mesh_list1: {len(mesh_list)}")
+
+# %%
+print(f"Shape of mesh_list1: {mesh_diff_matrix.shape}")
+
+# %% [markdown]
+# 4. Matrix difference - two arrays of different length
+# -----------------------------------------------------
+# Unlike difference_pairwise(), difference_matrix() can accept two lists/arrays of different lengths. As in Section 3, the differences between all possible pairs of meshes in both lists/arrays will be calculated.
+#
+# Consider following lists:
+#
+# .. math::
+#
+#   \text{List}_1 = \begin{bmatrix} A_1 & B_1 & C_1\\ \end{bmatrix}
+#
+# .. math::
+#
+#   \text{List}_2 = \begin{bmatrix} A_2 & B_2 \\ \end{bmatrix}
+#
+# The output array will contain following differences between meshes:
+#
+# .. math::
+#
+#   \begin{bmatrix} A_1-A_2 & A_1-B_2 \\ B_1-A_2 & B_1-B_2 \\ C_1-A_2 & C_1-B_2 \\ \end{bmatrix}
+#
+# and will have a shape of (len(mesh_list_matrix_1), len(mesh_list_matrix_2)). As in the following example:
+
+mesh_list_matrix_1 = [mesh1, mesh2, mesh1]
+mesh_list_matrix_2 = [mesh2, mesh1]
+
+mesh_diff_matrix = difference_matrix(
+    mesh_property, mesh_list_matrix_1, mesh_list_matrix_2
+)
+
+# %%
+print(f"Length of mesh_list_matrix_1: {len(mesh_list_matrix_1)}")
+
+# %%
+print(f"Length of mesh_list_matrix_2: {len(mesh_list_matrix_2)}")
+
+# %%
+print(f"Shape of mesh_diff_matrix: {mesh_diff_matrix.shape}")

docs/examples/howto_meshplotlib/plot_special.py

@@ -44,5 +44,5 @@ fig = plot_limit(mesh_series, si, "min")
 # Difference between the last and he first timestep:
 
 # %%
-diff_mesh = difference(mesh_series.read(-1), mesh_series.read(0), si)
+diff_mesh = difference(si, mesh_series.read(-1), mesh_series.read(0))
 fig = plot(diff_mesh, si.delta)

docs/examples/howto_meshplotlib/plot_with_custom_fig_ax.py

@@ -50,7 +50,7 @@ ax[0].set_title(r"$T(\mathrm{t}_{0})$")
 plot(meshseries.read(1), presets.temperature, fig=fig, ax=ax[1])
 ax[1].set_title(r"$T(\mathrm{t}_{end})$")
 diff_mesh = difference(
-    meshseries.read(0), meshseries.read(1), presets.temperature
+    presets.temperature, meshseries.read(1), meshseries.read(0)
 )
 plot(diff_mesh, presets.temperature.delta, fig=fig, ax=ax[2])
 ax[2].set_title(r"$T(\mathrm{t}_{end})$-$T(\mathrm{t}_{0})$")

docs/examples/howto_meshplotlib/plot_xy_labels_with_shared_axes.py

@@ -54,8 +54,8 @@ fig = plot([mesh_a, mesh_b], temperature)
 fig, ax = plt.subplots(2, 2, figsize=(40, 20), sharex=True, sharey=True)
 plot(mesh_a, temperature, fig=fig, ax=ax[0][0])
 plot(mesh_b, temperature, fig=fig, ax=ax[1][0])
-diff_ab = difference(mesh_b, mesh_a, temperature)
-diff_ba = difference(mesh_b, mesh_a, temperature)
+diff_ab = difference(temperature, mesh_b, mesh_a)
+diff_ba = difference(temperature, mesh_b, mesh_a)
 plot(diff_ab, temperature.delta, fig=fig, ax=ax[0][1])
 plot(diff_ba, temperature.delta, fig=fig, ax=ax[1][1])
 fig.tight_layout()

ogstools/meshlib/__init__.py

 from .boundary import Boundary, Layer, LocationFrame, Raster
 from .boundary_set import LayerSet
 from .boundary_subset import Surface
-from .data_processing import difference
+from .data_processing import difference, difference_matrix, difference_pairwise
 from .gmsh_meshing import cuboid, rect
 from .mesh_series import MeshSeries
 
@@ -14,6 +14,8 @@ __all__ = [
     "MeshSeries",
     "LayerSet",
     "difference",
+    "difference_pairwise",
+    "difference_matrix",
     "rect",
     "cuboid",
 ]

ogstools/meshlib/data_processing.py

-from typing import Union
+from itertools import product
+from typing import Optional, Union
 
+import numpy as np
 import pyvista as pv
 
 from ogstools.propertylib import Property, presets
 
 
 def difference(
+    mesh_property: Union[Property, str],
     mesh1: pv.UnstructuredGrid,
     mesh2: pv.UnstructuredGrid,
-    mesh_property: Union[Property, str],
 ) -> pv.UnstructuredGrid:
     """
     Compute the difference between two meshes based on a specified property.
 
+    :param mesh_property: The property to of interest.
     :param mesh1: The first mesh to be subtracted from.
     :param mesh2: The second mesh whose data is subtracted from the first mesh.
-    :param mesh_property: The property to of interest.
     :returns: A new mesh representing the difference between mesh1 and mesh2.
     """
     mesh_property = presets.get_preset(mesh_property, mesh1)
     diff_mesh = mesh1.copy(deep=True)
     diff_mesh[mesh_property.data_name] -= mesh2[mesh_property.data_name]
     return diff_mesh
+
+
+def difference_pairwise(
+    mesh_property: Union[Property, str],
+    meshes_1: Union[list, np.ndarray],
+    meshes_2: Union[list, np.ndarray],
+) -> np.ndarray:
+    """
+    Compute pairwise difference between meshes from two lists/arrays
+    (they have to be of the same length).
+
+    :param mesh_property: The property to of interest.
+    :param meshes_1: The first list/array of meshes to be subtracted from.
+    :param meshes_2: The second list/array of meshes whose data is subtracted
+                     from the first list/array of meshes - meshes_1.
+    :returns: An array of differences between meshes from meshes_1 and meshes_2.
+    """
+    meshes_1 = np.asarray(meshes_1).flatten()
+    meshes_2 = np.asarray(meshes_2).flatten()
+    if len(meshes_1) != len(meshes_2):
+        msg = "Mismatch in length of provided lists/arrays. \
+              Their length has to be identical to calculate pairwise \
+              difference. Did you intend to use difference_matrix()?"
+        raise RuntimeError(msg)
+    diff_mesh = [
+        difference(mesh_property, m1, m2) for m1, m2 in zip(meshes_1, meshes_2)
+    ]
+    return np.array(diff_mesh)
+
+
+def difference_matrix(
+    mesh_property: Union[Property, str],
+    meshes_1: Union[list, np.ndarray],
+    meshes_2: Optional[Union[list, np.ndarray]] = None,
+) -> np.ndarray:
+    """
+    Compute the difference between all combinations of two meshes
+    from one or two arrays based on a specified property.
+
+    :param mesh_property: The property to of interest.
+    :param meshes_1: The first list/array of meshes to be subtracted from.
+    :param meshes_2: The second list/array of meshes, it is subtracted from
+                     the first list/array of meshes - meshes_1 (optional).
+    :returns: An array of differences between meshes from meshes_1 and meshes_2.
+    """
+    if not isinstance(meshes_1, (list, np.ndarray)):
+        msg = "mesh1 is neither of type list nor np.ndarray"  # type: ignore[unreachable]
+        raise TypeError(msg)
+    meshes_1 = np.asarray(meshes_1).flatten()
+    if not isinstance(meshes_2, (list, np.ndarray)) and meshes_2 is not None:
+        msg = "mesh2 is neither of type list nor np.ndarray."  # type: ignore[unreachable]
+        raise TypeError(msg)
+    if meshes_2 is None:
+        meshes_2 = meshes_1.copy()
+    meshes_2 = np.asarray(meshes_2).flatten()
+    diff_mesh = [
+        difference(mesh_property, m1, m2)
+        for m1, m2 in product(meshes_1, meshes_2)
+    ]
+    return np.array(diff_mesh).reshape((len(meshes_1), len(meshes_2)))

ogstools/studies/templates/convergence_study.py

@@ -77,7 +77,7 @@ fig = meshplotlib.plot(richardson, mesh_property)
 data_key = mesh_property.data_name
 if reference_solution_path is None:
     diff_mesh = meshlib.difference(
-        richardson, topology.sample(meshes[-1]), mesh_property
+        mesh_property, richardson, topology.sample(meshes[-1])
     )
     fig = meshplotlib.plot(diff_mesh, mesh_property.delta)
 else:
@@ -85,7 +85,7 @@ else:
         meshlib.MeshSeries(reference_solution_path).read_closest(timevalue)
     )
     diff_mesh = meshlib.difference(
-        reference_solution, richardson, mesh_property
+        mesh_property, reference_solution, richardson
     )
     fig = meshplotlib.plot(diff_mesh, mesh_property.delta)
 

tests/test_meshlib.py

@@ -3,8 +3,17 @@
 import unittest
 
 import numpy as np
+from pyvista import UnstructuredGrid
 
-from ogstools.meshlib import MeshSeries, examples
+from ogstools.meshlib import (
+    MeshSeries,
+    difference,
+    difference_matrix,
+    difference_pairwise,
+    examples,
+)
+from ogstools.meshplotlib import examples as examples_mpl
+from ogstools.propertylib import presets
 
 
 class UtilsTest(unittest.TestCase):
@@ -45,3 +54,66 @@ class UtilsTest(unittest.TestCase):
         for method in ["nearest", "linear", None]:
             values = mesh_series.probe(points, "temperature", method)
             self.assertTrue(not np.any(np.isnan(values)))
+
+    def test_diff_two_meshes(self):
+        meshseries = examples_mpl.meshseries_THM_2D
+        mesh_property = presets.temperature
+        mesh1 = meshseries.read(0)
+        mesh2 = meshseries.read(-1)
+        mesh_diff = difference(mesh_property, mesh1, mesh2)
+        self.assertTrue(isinstance(mesh_diff, UnstructuredGrid))
+
+    def test_diff_pairwise(self):
+        n = 5
+        meshseries = examples_mpl.meshseries_THM_2D
+        mesh_property = presets.temperature
+        meshes1 = [meshseries.read(0)] * n
+        meshes2 = [meshseries.read(-1)] * n
+        meshes_diff = difference_pairwise(mesh_property, meshes1, meshes2)
+        self.assertTrue(
+            isinstance(meshes_diff, np.ndarray) and len(meshes_diff) == n
+        )
+
+    def test_diff_matrix_single_list(self):
+        meshseries = examples_mpl.meshseries_THM_2D
+        mesh_property = presets.temperature
+        meshes1 = [meshseries.read(0), meshseries.read(-1)]
+        meshes_diff = difference_matrix(mesh_property, meshes1)
+        self.assertTrue(
+            isinstance(meshes_diff, np.ndarray)
+            and meshes_diff.shape == (len(meshes1), len(meshes1))
+        )
+
+    def test_diff_matrix_single_numpy(self):
+        meshseries = examples_mpl.meshseries_THM_2D
+        mesh_property = presets.temperature
+        meshes1 = np.array([meshseries.read(0), meshseries.read(-1)])
+        meshes_diff = difference_matrix(mesh_property, meshes1)
+        self.assertTrue(
+            isinstance(meshes_diff, np.ndarray)
+            and meshes_diff.shape == (len(meshes1), len(meshes1))
+        )
+
+    def test_diff_matrix_unequal_list(self):
+        meshseries = examples_mpl.meshseries_THM_2D
+        mesh_property = presets.temperature
+        meshes1 = [meshseries.read(0), meshseries.read(-1)]
+        meshes2 = [meshseries.read(0), meshseries.read(-1), meshseries.read(-1)]
+        meshes_diff = difference_matrix(mesh_property, meshes1, meshes2)
+        self.assertTrue(
+            isinstance(meshes_diff, np.ndarray)
+            and meshes_diff.shape == (len(meshes1), len(meshes2))
+        )
+
+    def test_diff_matrix_unequal_numpy(self):
+        meshseries = examples_mpl.meshseries_THM_2D
+        mesh_property = presets.temperature
+        meshes1 = np.array([meshseries.read(0), meshseries.read(-1)])
+        meshes2 = np.array(
+            [meshseries.read(0), meshseries.read(-1), meshseries.read(-1)]
+        )
+        meshes_diff = difference_matrix(mesh_property, meshes1, meshes2)
+        self.assertTrue(
+            isinstance(meshes_diff, np.ndarray)
+            and meshes_diff.shape == (len(meshes1), len(meshes2))
+        )

tests/test_meshplotlib.py

@@ -135,7 +135,7 @@ class MeshplotlibTest(unittest.TestCase):
     def test_diff_plots(self):
         """Test creation of difference plots."""
         meshseries = examples.meshseries_CT_2D
-        plot(difference(meshseries.read(0), meshseries.read(1), "Si"), "Si")
+        plot(difference("Si", meshseries.read(0), meshseries.read(1)), "Si")
 
     def test_user_defined_ax(self):
         """Test creating plot with subfigures and user provided ax"""
@@ -146,7 +146,7 @@ class MeshplotlibTest(unittest.TestCase):
         plot(meshseries.read(1), presets.temperature, fig=fig, ax=ax[1])
         ax[1].set_title(r"$T(\mathrm{t}_{end})$")
         diff_mesh = difference(
-            meshseries.read(0), meshseries.read(1), presets.temperature
+            presets.temperature, meshseries.read(0), meshseries.read(1)
         )
         plot(diff_mesh, presets.temperature.delta, fig=fig, ax=ax[2])
         ax[2].set_title(r"$T(\mathrm{t}_{end})$-$T(\mathrm{t}_{0})$")
@@ -197,8 +197,8 @@ class MeshplotlibTest(unittest.TestCase):
         ax = ax.flatten()
         plot(meshseries.read(0), presets.temperature, fig=fig, ax=ax[0])
         plot(meshseries.read(1), presets.temperature, fig=fig, ax=ax[1])
-        diff_ab = difference(mesh_a, mesh_b, presets.temperature)
-        diff_ba = difference(mesh_b, mesh_a, presets.temperature)
+        diff_ab = difference(presets.temperature, mesh_a, mesh_b)
+        diff_ba = difference(presets.temperature, mesh_b, mesh_a)
         plot(diff_ab, presets.temperature.delta, fig=fig, ax=ax[2])
         plot(diff_ba, presets.temperature.delta, fig=fig, ax=ax[3])
         fig.tight_layout()
@@ -211,8 +211,8 @@ class MeshplotlibTest(unittest.TestCase):
         fig, ax = plt.subplots(2, 2, sharex=True, sharey=True)
         plot(meshseries.read(0), presets.temperature, fig=fig, ax=ax[0][0])
         plot(meshseries.read(1), presets.temperature, fig=fig, ax=ax[1][0])
-        diff_ab = difference(mesh_a, mesh_b, presets.temperature)
-        diff_ba = difference(mesh_b, mesh_a, presets.temperature)
+        diff_ab = difference(presets.temperature, mesh_a, mesh_b)
+        diff_ba = difference(presets.temperature, mesh_b, mesh_a)
         plot(diff_ab, presets.temperature.delta, fig=fig, ax=ax[0][1])
         plot(diff_ba, presets.temperature.delta, fig=fig, ax=ax[1][1])
         ax = label_spatial_axes(ax, np.array([0, 1]))