diff --git a/.gitignore b/.gitignore
index f0aa73237582ffff44dc831a95971f0673e41590..3e4a206410866b236da1b1f761cba8be3881393f 100644
--- a/.gitignore
+++ b/.gitignore
@@ -162,6 +162,9 @@ cython_debug/
# be found at https://github.com/github/gitignore/blob/main/Global/JetBrains.gitignore
# and can be added to the global gitignore or merged into this file. For a more nuclear
# option (not recommended) you can uncomment the following to ignore the entire idea folder.
-#.idea/
+#.idea/S
/docs/auto_examples/
+
+# ignore jupyter notebooks in main folder (useful for testing, annoying for git)
+/**.ipynb
diff --git a/docs/examples/howto_meshlib/plot_calculate_diff.py b/docs/examples/howto_meshlib/plot_calculate_diff.py
index ea6828acf646ca0162e1450072415062d89256b2..294e6df557a4730bc8caecad1561c64102258529 100644
--- a/docs/examples/howto_meshlib/plot_calculate_diff.py
+++ b/docs/examples/howto_meshlib/plot_calculate_diff.py
@@ -4,7 +4,8 @@ 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.
+This example explains how to use functions from meshlib to calculate differences
+between meshes.
"""
# %%
@@ -25,25 +26,31 @@ 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.
+# 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.
+# 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:
+# 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)
+mesh_diff = difference(mesh1, mesh2, mesh_property)
# %% [markdown]
# This returned object will be a PyVista UnstructuredGrid object:
@@ -54,7 +61,9 @@ 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.
+# 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:
#
@@ -79,7 +88,7 @@ print(f"Type of mesh_diff: {type(mesh_diff)}")
meshes_1 = [mesh1] * 3
meshes_2 = [mesh2] * 3
-mesh_diff_pair_wise = difference_pairwise(mesh_property, meshes_1, meshes_2)
+mesh_diff_pair_wise = difference_pairwise(meshes_1, meshes_2, mesh_property)
# %%
print(f"Length of mesh_list1: {len(meshes_1)}")
@@ -93,7 +102,8 @@ 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.
+# 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:
#
@@ -107,11 +117,12 @@ print(f"Shape of mesh_diff_pair_wise: {mesh_diff_pair_wise.shape}")
#
# \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:
+# 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)
+mesh_diff_matrix = difference_matrix(mesh_list, mesh_property=mesh_property)
# %%
print(f"Length of mesh_list1: {len(mesh_list)}")
@@ -122,7 +133,9 @@ 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.
+# 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:
#
@@ -140,13 +153,14 @@ print(f"Shape of mesh_list1: {mesh_diff_matrix.shape}")
#
# \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:
+# 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
+ mesh_list_matrix_1, mesh_list_matrix_2, mesh_property
)
# %%
diff --git a/docs/examples/howto_meshplotlib/plot_animation.py b/docs/examples/howto_meshplotlib/plot_animation.py
index 47ff7a8612abdb2ab61966b8e69fb3079d248649..7c0ad10d26f2218f3fa371a3808617361224bc64 100644
--- a/docs/examples/howto_meshplotlib/plot_animation.py
+++ b/docs/examples/howto_meshplotlib/plot_animation.py
@@ -17,37 +17,59 @@ from ogstools.propertylib import Scalar
setup.reset()
mesh_series = examples.meshseries_CT_2D
-# alternatively:
-# from ogstools.meshlib import MeshSeries
-# mesh_series = MeshSeries("filepath/filename_pvd_or_xdmf")
-# %%
+mesh_property = Scalar(
+ data_name="Si", data_unit="", output_unit="%", output_name="Saturation"
+)
+
+# %% [markdown]
+# To read your own data as a mesh series you can do:
+#
+# .. code-block:: python
+#
+# from ogstools.meshlib import MeshSeries
+# mesh_series = MeshSeries("filepath/filename_pvd_or_xdmf")
+#
+# You can also use a property from the available presets instead of needing to
+# create your own:
+# :ref:`sphx_glr_auto_examples_howto_propertylib_plot_propertylib.py`
+
+# %% [markdown]
# Let's use fixed scale limits to prevent rescaling during the animation.
+
+# %%
setup.p_min = 0
setup.p_max = 100
setup.dpi = 50
-# %%
+# %% [markdown]
# You can choose which timesteps to render by passing either an int array
# corresponding to the indices, or a float array corresponding to the timevalues
# to render. If a requested timevalue is not part of the timeseries it will be
# interpolated. In this case every second frame will be interpolated.
+
+# %%
timevalues = np.linspace(
mesh_series.timevalues[0], mesh_series.timevalues[-1], num=25
)
-# %%
+# %% [markdown]
# Now, let's animate the saturation solution. A timescale at the top
# indicates existing timesteps and the position of the current timevalue.
# Note that rendering many frames in conjunction with large meshes might take
# a really long time.
+
+# %%
titles = [f"{tv/(365.25*86400):.1f} yrs" for tv in timevalues]
-si = Scalar(
- data_name="Si", data_unit="", output_unit="%", output_name="Saturation"
-)
-anim = animate(mesh_series, si, timevalues, titles)
-# the animation can be saved (as mp4) like so:
-# from ogstools.meshplotlib.animation import save_animation
-# save_animation(anim, "Saturation", fps=5)
+anim = animate(mesh_series, mesh_property, timevalues, titles)
+
+# %% [markdown]
+# The animation can be saved (as mp4) like so:
+#
+# .. code-block:: python
+#
+# from ogstools.meshplotlib.animation import save_animation
+# save_animation(anim, "Saturation", fps=5)
+#
# sphinx_gallery_start_ignore
# note for developers:
diff --git a/docs/examples/howto_meshplotlib/plot_meshplotlib_2d.py b/docs/examples/howto_meshplotlib/plot_meshplotlib_2d.py
index 363064503006f47f1256d4fce93449625443648b..862b88355125ab2b6b3269f9a5d52c90dcbeab16 100644
--- a/docs/examples/howto_meshplotlib/plot_meshplotlib_2d.py
+++ b/docs/examples/howto_meshplotlib/plot_meshplotlib_2d.py
@@ -4,10 +4,10 @@ Visualizing 2D model data
.. sectionauthor:: Florian Zill (Helmholtz Centre for Environmental Research GmbH - UFZ)
-For this example we load a 2D meshseries dataset from within the ``meshplotlib`` examples.
-In the ``meshplotlib.setup`` we can provide a dictionary to map names to material ids.
-First, let's plot the material ids (cell_data). Per default in the setup, this
-will automatically show the element edges.
+For this example we load a 2D meshseries from within the ``meshplotlib``
+examples. In the ``meshplotlib.setup`` we can provide a dictionary to map names
+to material ids. First, let's plot the material ids (cell_data). Per default in
+the setup, this will automatically show the element edges.
"""
# %%
@@ -19,34 +19,52 @@ mpl.setup.reset()
mpl.setup.length.output_unit = "km"
mpl.setup.material_names = {i + 1: f"Layer {i+1}" for i in range(26)}
mesh = meshseries_THM_2D.read(1)
-fig = mpl.plot(mesh, presets.material_id)
+
+# %% [markdown]
+# To read your own data as a mesh series you can do:
+#
+# .. code-block:: python
+#
+# from ogstools.meshlib import MeshSeries
+# mesh_series = MeshSeries("filepath/filename_pvd_or_xdmf")
+#
# %%
+fig = mpl.plot(mesh, presets.material_id)
+
+# %% [markdown]
# Now, let's plot the temperature field (point_data) at the first timestep.
# The default temperature property from the `propertylib` reads the temperature
# data as Kelvin and converts them to degrees Celsius.
+# %%
fig = mpl.plot(mesh, presets.temperature)
-# %%
+# %% [markdown]
# We can also plot components of vector properties:
+# %%
fig = mpl.plot(mesh, presets.displacement[0])
# %%
fig = mpl.plot(mesh, presets.displacement[1])
-# %%
+# %% [markdown]
# This example has hydraulically deactivated subdomains:
+# %%
fig = mpl.plot(mesh, presets.pressure.get_mask())
-# %%
+# %% [markdown]
# Let's plot the fluid velocity field.
-fig = mpl.plot(mesh, presets.velocity)
# %%
+fig = mpl.plot(mesh, presets.velocity)
+
+# %% [markdown]
# Let's plot it again, this time log-scaled.
+
+# %%
mpl.setup.log_scaled = True
mpl.setup.p_min = -8
fig = mpl.plot(mesh, presets.velocity)
diff --git a/docs/examples/howto_meshplotlib/plot_meshplotlib_3d.py b/docs/examples/howto_meshplotlib/plot_meshplotlib_3d.py
index 6ef0975e3674bd81df777c7b66dba76cff3f2580..8c9e6a70857770c15db9d5d6a757b0f260a082c1 100644
--- a/docs/examples/howto_meshplotlib/plot_meshplotlib_3d.py
+++ b/docs/examples/howto_meshplotlib/plot_meshplotlib_3d.py
@@ -9,7 +9,7 @@ To create them we use ``pyvista`` as it provides all the necessary functionality
If we want to plot multiple meshes, they have to be in a 2D ``numpy`` array.
The shape of this array determines the number of rows and columns in our plot.
First, let's load 3D example data from ``pyvista`` and plot the only available
-dataset: "facies" with a native pyvista plot.
+data: "facies" with a native pyvista plot.
"""
# %%
diff --git a/docs/examples/howto_meshplotlib/plot_observation_points.py b/docs/examples/howto_meshplotlib/plot_observation_points.py
index 70b6331cd2875d0deb1cc1173b66de2bdd714021..91773041794a19e04ba0dfd1cbcbe6336f7500d5 100644
--- a/docs/examples/howto_meshplotlib/plot_observation_points.py
+++ b/docs/examples/howto_meshplotlib/plot_observation_points.py
@@ -33,9 +33,18 @@ meshplotlib.setup.reset()
si = Scalar(
data_name="Si", data_unit="", output_unit="%", output_name="Saturation"
)
-# alternatively:
-# from ogstools.meshlib import MeshSeries
-# mesh_series = MeshSeries("filepath/filename_pvd_or_xdmf")
+
+# %% [markdown]
+# To read your own data as a mesh series you can do:
+#
+# .. code-block:: python
+#
+# from ogstools.meshlib import MeshSeries
+# mesh_series = MeshSeries("filepath/filename_pvd_or_xdmf")
+#
+# You can also use a property from the available presets instead of needing to
+# create your own:
+# :ref:`sphx_glr_auto_examples_howto_propertylib_plot_propertylib.py`
# %% [markdown]
# Let's define 4 observation points and plot them on the mesh.
diff --git a/docs/examples/howto_meshplotlib/plot_special.py b/docs/examples/howto_meshplotlib/plot_special.py
index 48ae7d864291993e0355e838ad01f85bb992117c..41e320deb993afb391cfda04753236da1d25d2e2 100644
--- a/docs/examples/howto_meshplotlib/plot_special.py
+++ b/docs/examples/howto_meshplotlib/plot_special.py
@@ -1,11 +1,11 @@
"""
-How to plot limits and differences
-==================================
+Analyzing Meshseries Data
+=========================
.. sectionauthor:: Florian Zill (Helmholtz Centre for Environmental Research GmbH - UFZ)
-In this example we plot the limits of all nodes in a model through all
-timesteps, as well as differences between to timesteps. For this purpose we use
+In this example we show how to aggregate data in a model over all timesteps as
+well as plot differences between to timesteps. For this purpose we use
a component transport example from the ogs benchmark gallery
(https://www.opengeosys.org/docs/benchmarks/hydro-component/elder/).
@@ -15,7 +15,7 @@ To see this benchmark results over all timesteps have a look at
# %%
from ogstools.meshlib import difference
-from ogstools.meshplotlib import examples, plot, plot_limit, setup
+from ogstools.meshplotlib import examples, plot, setup
from ogstools.propertylib import Scalar
setup.reset()
@@ -23,26 +23,44 @@ mesh_series = examples.meshseries_CT_2D
si = Scalar(
data_name="Si", data_unit="", output_unit="%", output_name="Saturation"
)
-# alternatively:
-# from ogstools.meshlib import MeshSeries
-# mesh_series = MeshSeries("filepath/filename_pvd_or_xdmf")
# %% [markdown]
-# Maximum though all timesteps
+# To read your own data as a mesh series you can do:
+#
+# .. code-block:: python
+#
+# from ogstools.meshlib import MeshSeries
+# mesh_series = MeshSeries("filepath/filename_pvd_or_xdmf")
+#
+# You can also use a property from the available presets instead of needing to
+# create your own:
+# :ref:`sphx_glr_auto_examples_howto_propertylib_plot_propertylib.py`
+
+# %% [markdown]
+# You aggregate the data in MeshSeries over all timesteps given some
+# aggregation function, e.g. "min", "max", "var"
+# (see: :meth:`~ogstools.meshlib.mesh_series.MeshSeries.aggregate`).
+# The following code gets the maximum saturation for each point in the mesh over
+# all timesteps and plots it. Note: the aggragation function returns the
+# resulting mesh and a new property corresponding to the aggregated data, to
+# ensure no duplicate unit conversion or data transformation is happening in the
+# plot function.
# %%
-fig = plot_limit(mesh_series, si, "max")
+mesh = mesh_series.aggregate(si, "max")
+fig = plot(mesh, si)
# %% [markdown]
-# Minimum though all timesteps
+# Likewise we can calculate and visualize the variance of the saturation:
-# %%
-fig = plot_limit(mesh_series, si, "min")
+# %%
+mesh = mesh_series.aggregate(si, "var")
+fig = plot(mesh, si)
# %% [markdown]
-# Difference between the last and he first timestep:
+# Difference between the last and the first timestep:
# %%
-diff_mesh = difference(si, mesh_series.read(-1), mesh_series.read(0))
-fig = plot(diff_mesh, si.delta)
+mesh = difference(mesh_series.read(-1), mesh_series.read(0), si)
+fig = plot(mesh, si)
diff --git a/docs/examples/howto_meshplotlib/plot_with_custom_fig_ax.py b/docs/examples/howto_meshplotlib/plot_with_custom_fig_ax.py
index f3e5b8a020c1f6f2ceb8c6dbbea22be44e9a7b53..0c1cc5c6a551841366599a242208db84da65d47a 100644
--- a/docs/examples/howto_meshplotlib/plot_with_custom_fig_ax.py
+++ b/docs/examples/howto_meshplotlib/plot_with_custom_fig_ax.py
@@ -4,7 +4,7 @@ Plotting different process variables on already existing Matplotlib figures / ax
.. sectionauthor:: Feliks Kiszkurno (Helmholtz Centre for Environmental Research GmbH - UFZ)
-For this example we load a 2D meshseries dataset from within the ``meshplotlib``
+For this example we load a 2D meshseries from within the ``meshplotlib``
examples. This tutorial covers using meshplotlib to plot meshseries data using
Matplotlib objects for Figure and / or Axis. This is useful if different
plotting functions from Meshplotlib are to be used on different subplots
@@ -50,9 +50,9 @@ 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(
- presets.temperature, meshseries.read(1), meshseries.read(0)
+ meshseries.read(1), meshseries.read(0), presets.temperature
)
-plot(diff_mesh, presets.temperature.delta, fig=fig, ax=ax[2])
+plot(diff_mesh, presets.temperature, fig=fig, ax=ax[2])
ax[2].set_title(r"$T(\mathrm{t}_{end})$-$T(\mathrm{t}_{0})$")
fig.suptitle("Plot two time steps and their difference - with colorbars")
fig.tight_layout()
diff --git a/docs/examples/howto_meshplotlib/plot_xy_labels_with_shared_axes.py b/docs/examples/howto_meshplotlib/plot_xy_labels_with_shared_axes.py
index 26773da7ac9d81b1781e4917dea35fde0ab39042..ded8afdfb74b569656f56e6a4f11efd81858c61c 100644
--- a/docs/examples/howto_meshplotlib/plot_xy_labels_with_shared_axes.py
+++ b/docs/examples/howto_meshplotlib/plot_xy_labels_with_shared_axes.py
@@ -6,7 +6,7 @@ Labeling directional shared axes
.. warning:: This example discusses functionality that may experience breaking changes in the near future!
-For this example we load a 2D meshseries dataset from within the ``meshplotlib``
+For this example we load a 2D meshseries from within the ``meshplotlib``
examples. This tutorial covers automatic labeling the directional axes (X and Y)
under various conditions (shared and nor shared X and Y axes).
"""
@@ -14,7 +14,6 @@ under various conditions (shared and nor shared X and Y axes).
# %%
# Import Python packages, change some settings and load example data set
import matplotlib.pyplot as plt
-import numpy as np
from ogstools.meshlib import difference
from ogstools.meshplotlib import (
@@ -54,10 +53,10 @@ 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(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])
+diff_ab = difference(mesh_a, mesh_b, temperature)
+diff_ba = difference(mesh_b, mesh_a, temperature)
+plot(diff_ab, temperature, fig=fig, ax=ax[0][1])
+plot(diff_ba, temperature, fig=fig, ax=ax[1][1])
fig.tight_layout()
# %%
@@ -69,8 +68,8 @@ fig.tight_layout()
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])
-plot(diff_ab, temperature.delta, fig=fig, ax=ax[0][1])
-plot(diff_ba, temperature.delta, fig=fig, ax=ax[1][1])
-ax = clear_labels(ax)
-ax = label_spatial_axes(ax, np.array([0, 1]))
+plot(difference(mesh_a, mesh_b, temperature), temperature, fig=fig, ax=ax[0][1])
+plot(difference(mesh_b, mesh_a, temperature), temperature, fig=fig, ax=ax[1][1])
+clear_labels(ax)
+label_spatial_axes(ax)
fig.tight_layout()
diff --git a/ogstools/meshlib/boundary.py b/ogstools/meshlib/boundary.py
index d298cc7e6ff93daea1865a6c9822f46c4d76f0a0..030e3fb9511b3555e8b55b5f35273306be9bc89a 100644
--- a/ogstools/meshlib/boundary.py
+++ b/ogstools/meshlib/boundary.py
@@ -21,7 +21,7 @@ class Boundary(ABC):
"""
@abstractmethod
- def dim(self):
+ def dim(self) -> int:
"""
Get the dimension of the boundary.
@@ -29,7 +29,7 @@ class Boundary(ABC):
int: The dimension of the boundary. For example, the dimension of a boundary
of a cube (3D) is 2.
"""
- return
+ return 0
@dataclass(frozen=True)
@@ -47,7 +47,7 @@ class Layer(Boundary):
of the layer in the GIS system.
"""
- def __post_init__(self):
+ def __post_init__(self) -> None:
if not self.material_id:
object.__setattr__(self, "material_id", self.bottom._material_id)
@@ -89,7 +89,7 @@ class Layer(Boundary):
rasters.append(bottom_raster)
return rasters
- def dim(self):
+ def dim(self) -> int:
return 3
@@ -100,7 +100,7 @@ class LocationFrame:
ymin: float
ymax: float
- def as_gml(self, filename: Path):
+ def as_gml(self, filename: Path) -> None:
"""
Generate GML representation of the location frame.
diff --git a/ogstools/meshlib/boundary_set.py b/ogstools/meshlib/boundary_set.py
index 17dd54b920616e3940536ecdd6ca417f6148a9fa..22de2ef82abdb80a0b727b610139f19819b0abc4 100644
--- a/ogstools/meshlib/boundary_set.py
+++ b/ogstools/meshlib/boundary_set.py
@@ -20,12 +20,12 @@ class BoundarySet(ABC):
"""
@abstractmethod
- def bounds(self):
- return
+ def bounds(self) -> list:
+ return []
@abstractmethod
- def filenames(self):
- return
+ def filenames(self) -> list[Path]:
+ return []
class LayerSet(BoundarySet):
@@ -49,16 +49,16 @@ class LayerSet(BoundarySet):
raise ValueError(msg)
self.layers = layers
- def bounds(self):
+ def bounds(self) -> list:
return list(self.layers[0].top.mesh.bounds)
- def filenames(self):
+ def filenames(self) -> list[Path]:
layer_filenames = [layer.bottom.filename for layer in self.layers]
layer_filenames.insert(0, self.layers[0].top.filename) # file interface
return layer_filenames
@classmethod
- def from_pandas(cls, df: pd.DataFrame):
+ def from_pandas(cls, df: pd.DataFrame) -> "LayerSet":
"""Create a LayerSet from a Pandas DataFrame."""
Row = namedtuple("Row", ["material_id", "mesh", "resolution"])
surfaces = [
@@ -79,7 +79,7 @@ class LayerSet(BoundarySet):
]
return cls(layers=base_layer)
- def create_raster(self, resolution):
+ def create_raster(self, resolution: float) -> tuple[Path, Path]:
"""
Create raster representations for the LayerSet.
@@ -103,7 +103,7 @@ class LayerSet(BoundarySet):
file.write("\n".join(str(item) for item in raster_set))
return raster_vtu, rastered_layers_txt
- def create_rasters(self, resolution: int) -> list[Path]:
+ def create_rasters(self, resolution: float) -> list[Path]:
"""
For each surface a (temporary) raster file with given resolution is created.
"""
diff --git a/ogstools/meshlib/boundary_subset.py b/ogstools/meshlib/boundary_subset.py
index 6a7f7b6a2c5d2bdd9569c227d61bdd7bcd596826..ef74fc93e93c7eb3323076fb364ecfda29e3fc39 100644
--- a/ogstools/meshlib/boundary_subset.py
+++ b/ogstools/meshlib/boundary_subset.py
@@ -5,6 +5,7 @@ from typing import Union
import numpy as np
import pyvista as pv
from ogs import cli
+from typeguard import typechecked
class Surface:
@@ -19,11 +20,8 @@ class Surface:
def material_id(self) -> int:
return self._material_id
- def __init__(
- self,
- input: Union[Path, pv.DataSet],
- material_id: int,
- ):
+ @typechecked
+ def __init__(self, input: Union[Path, pv.DataObject], material_id: int):
"""Initialize a surface mesh. Either from pyvista or from a file."""
self._material_id = material_id
@@ -33,21 +31,16 @@ class Surface:
msg = f"{self.filename} does not exist."
raise ValueError(msg)
self.mesh = pv.get_reader(self.filename).read()
- elif isinstance(input, pv.DataSet):
+ elif isinstance(input, pv.DataObject):
self.mesh = input
self.filename = Path(tempfile.mkstemp(".vtu", "surface")[1])
pv.save_meshio(self.filename, self.mesh, file_format="vtu")
- else:
- msg = "{} given, must be either Path or pyvista Dataset.".format(
- self.filename
- )
- raise ValueError(msg)
self.mesh.cell_data["MaterialIDs"] = (
np.ones(self.mesh.n_cells) * self.material_id
).astype(np.int32)
- def __eq__(self, other):
+ def __eq__(self, other: object) -> bool:
return self.__dict__ == other.__dict__
def create_raster_file(self, resolution: float) -> Path:
@@ -81,7 +74,7 @@ def Gaussian2D(
spread: float,
height_offset: float,
n: int,
-):
+) -> pv.DataSet:
"""
Generate a 2D Gaussian-like surface using the provided parameters.
diff --git a/ogstools/meshlib/data_processing.py b/ogstools/meshlib/data_processing.py
index 61f359d4ee389f19daf7e7a3bfae8020d4ecfcbb..9b65e77d4290694aaa75c9f83f2e701d81afbaff 100644
--- a/ogstools/meshlib/data_processing.py
+++ b/ogstools/meshlib/data_processing.py
@@ -3,43 +3,73 @@ from typing import Optional, Union
import numpy as np
import pyvista as pv
+from typeguard import typechecked
-from ogstools.propertylib import Property, presets
+from ogstools.propertylib import Property
+
+
+def _raw_differences_all_data(
+ mesh1: pv.UnstructuredGrid, mesh2: pv.UnstructuredGrid
+) -> pv.UnstructuredGrid:
+ diff_mesh = mesh1.copy(deep=True)
+ for point_data_key in mesh1.point_data:
+ diff_mesh.point_data[point_data_key] -= mesh2.point_data[point_data_key]
+ for cell_data_key in mesh1.cell_data:
+ if cell_data_key == "MaterialIDs":
+ continue
+ diff_mesh.cell_data[cell_data_key] -= mesh2.cell_data[cell_data_key]
+ return diff_mesh
def difference(
- mesh_property: Union[Property, str],
mesh1: pv.UnstructuredGrid,
mesh2: pv.UnstructuredGrid,
+ mesh_property: Optional[Union[Property, str]] = None,
) -> pv.UnstructuredGrid:
"""
- Compute the difference between two meshes based on a specified property.
+ Compute the difference of properties between two meshes.
- :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.
- :returns: A new mesh representing the difference between mesh1 and mesh2.
+ :param mesh_property: The property of interest. If not given, all point
+ and cell_data will be processed raw.
+ :returns: A new mesh containing the difference of `mesh_property` or all
+ datasets between mesh1 and mesh2.
"""
- mesh_property = presets.get_preset(mesh_property, mesh1)
+ if mesh_property is None:
+ return _raw_differences_all_data(mesh1, mesh2)
+ if isinstance(mesh_property, Property):
+ vals = np.asarray(
+ [mesh_property.transform(mesh) for mesh in [mesh1, mesh2]]
+ )
+ outname = mesh_property.output_name + "_difference"
+ else:
+ vals = np.asarray([mesh[mesh_property] for mesh in [mesh1, mesh2]])
+ outname = mesh_property + "_difference"
+
diff_mesh = mesh1.copy(deep=True)
- diff_mesh[mesh_property.data_name] -= mesh2[mesh_property.data_name]
+ diff_mesh.clear_data()
+ diff_mesh[outname] = np.empty(vals.shape[1:])
+ diff_mesh[outname] = vals[0] - vals[1]
return diff_mesh
def difference_pairwise(
- mesh_property: Union[Property, str],
meshes_1: Union[list, np.ndarray],
meshes_2: Union[list, np.ndarray],
+ mesh_property: Optional[Union[Property, str]] = None,
) -> 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.
+ :param mesh_property: The property of interest. If not given, all point
+ and cell_data will be processed raw.
+ :returns: An array of meshes containing the differences of `mesh_property`
+ or all datasets between meshes_1 and meshes_2.
"""
meshes_1 = np.asarray(meshes_1).flatten()
meshes_2 = np.asarray(meshes_2).flatten()
@@ -48,39 +78,39 @@ def difference_pairwise(
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)
+ return np.asarray(
+ [
+ difference(m1, m2, mesh_property)
+ for m1, m2 in zip(meshes_1, meshes_2)
+ ]
+ )
+@typechecked
def difference_matrix(
- mesh_property: Union[Property, str],
meshes_1: Union[list, np.ndarray],
meshes_2: Optional[Union[list, np.ndarray]] = None,
+ mesh_property: Optional[Union[Property, str]] = 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.
+ :param mesh_property: The property of interest. If not given, all point
+ and cell_data will be processed raw.
+ :returns: An array of meshes containing the differences of `mesh_property`
+ or all datasets between meshes_1 and meshes_2 for all possible
+ combinations.
"""
- 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)
+ diff_meshes = [
+ difference(m1, m2, mesh_property)
for m1, m2 in product(meshes_1, meshes_2)
]
- return np.array(diff_mesh).reshape((len(meshes_1), len(meshes_2)))
+ return np.asarray(diff_meshes).reshape((len(meshes_1), len(meshes_2)))
diff --git a/ogstools/meshlib/gmsh_meshing.py b/ogstools/meshlib/gmsh_meshing.py
index c58cb076da802f1cbb6ab3a3c85a94de9ec3fa93..c985a5000fd5119c0d45e048635b745c3c77b75d 100644
--- a/ogstools/meshlib/gmsh_meshing.py
+++ b/ogstools/meshlib/gmsh_meshing.py
@@ -9,7 +9,7 @@ def _geo_square(
lengths: Union[float, list[float]],
n_edge_cells: Union[int, list[int]],
structured: bool,
-):
+) -> None:
_lengths = lengths if isinstance(lengths, list) else [lengths] * 2
_n = n_edge_cells if isinstance(n_edge_cells, list) else [n_edge_cells] * 2
geo.addPoint(0, 0, 0, tag=1)
@@ -41,7 +41,7 @@ def rect(
structured_grid: bool = True,
order: int = 1,
out_name: Path = Path("unit_square.msh"),
-):
+) -> None:
gmsh.initialize()
gmsh.option.set_number("General.Verbosity", 0)
gmsh.model.add("unit_square")
@@ -75,7 +75,7 @@ def cuboid(
structured_grid: bool = True,
order: int = 1,
out_name: Path = Path("unit_cube.msh"),
-):
+) -> None:
gmsh.initialize()
gmsh.option.set_number("General.Verbosity", 0)
gmsh.model.add("unit_cube")
@@ -118,7 +118,7 @@ def bhe_mesh(
bhe_depth: float = 20,
order: int = 1,
out_name: Path = Path("bhe_mesh.msh"),
-):
+) -> None:
gmsh.initialize()
model, geo = (gmsh.model, gmsh.model.geo)
model.add(Path(out_name).stem)
diff --git a/ogstools/meshlib/mesh_series.py b/ogstools/meshlib/mesh_series.py
index d43f396541e919545b9263f3cf3650a2d2c620bb..c7dd4a5d929794225ba21ec0ca9c5415c8cedabc 100644
--- a/ogstools/meshlib/mesh_series.py
+++ b/ogstools/meshlib/mesh_series.py
@@ -7,9 +7,12 @@ import meshio
import numpy as np
import pyvista as pv
import vtuIO
+from h5py import File
from scipy.interpolate import LinearNDInterpolator, NearestNDInterpolator
from tqdm.auto import tqdm
+from ogstools.propertylib import Property
+
from .xdmf_reader import XDMFReader
@@ -53,7 +56,7 @@ class MeshSeries:
raise TypeError(msg)
@property
- def hdf5(self):
+ def hdf5(self) -> File:
# We assume there is only one h5 file
return next(iter(self._xdmf_reader.hdf5_files.values()))
@@ -164,13 +167,66 @@ class MeshSeries:
)
return mesh[data_name]
+ def aggregate(
+ self,
+ mesh_property: Union[Property, str],
+ func: Literal["min", "max", "mean", "median", "sum", "std", "var"],
+ ) -> pv.UnstructuredGrid:
+ """Aggregate data over all timesteps using a specified function.
+
+ :param mesh_property:
+ The mesh property to be aggregated. If given as type `Property`, the
+ :meth:`~ogstools.propertylib.property.Property.transform` function
+ will be applied on each timestep and aggregation afterwards.
+ :param func:
+ The aggregation function to apply. It must be one of "min", "max",
+ "mean", "median", "sum", "std", "var". The equally named numpy
+ function will be used to aggregate over all timesteps.
+ :returns: A mesh with aggregated data according to the given function.
+
+ """
+ np_func = {
+ "min": np.min,
+ "max": np.max,
+ "mean": np.mean,
+ "median": np.median,
+ "sum": np.sum,
+ "std": np.std,
+ "var": np.var,
+ }[func]
+ mesh = self.read(0).copy(deep=True)
+ mesh.clear_data()
+ if isinstance(mesh_property, Property):
+ if mesh_property.mesh_dependent:
+ vals = np.asarray(
+ [
+ mesh_property.transform(self.read(t))
+ for t in tqdm(self.timesteps)
+ ]
+ )
+ else:
+ vals = mesh_property.transform(
+ self.values(mesh_property.data_name)
+ )
+ else:
+ vals = self.values(mesh_property)
+ output_name = (
+ f"{mesh_property.output_name}_{func}"
+ if isinstance(mesh_property, Property)
+ else f"{mesh_property}_{func}"
+ )
+ mesh[output_name] = np.empty(vals.shape[1:])
+ assert isinstance(np_func, type(np.max))
+ np_func(vals, out=mesh[output_name], axis=0)
+ return mesh
+
def _probe_pvd(
self,
points: np.ndarray,
data_name: str,
interp_method: Optional[Literal["nearest", "probefilter"]] = None,
interp_backend: Optional[Literal["vtk", "scipy"]] = None,
- ):
+ ) -> np.ndarray:
obs_pts_dict = {f"pt{j}": point for j, point in enumerate(points)}
dim = self.read(0).get_cell(0).dimension
pvd_path = self.filepath
@@ -187,7 +243,7 @@ class MeshSeries:
points: np.ndarray,
data_name: str,
interp_method: Optional[Literal["nearest", "linear"]] = None,
- ):
+ ) -> np.ndarray:
values = self.hdf5["meshes"][self.hdf5_bulk_name][data_name][:]
geom = self.hdf5["meshes"][self.hdf5_bulk_name]["geometry"][0]
values = np.swapaxes(values, 0, 1)
diff --git a/ogstools/meshlib/region/region.py b/ogstools/meshlib/region/region.py
index 840542d9011a04dee3b767932f32e1983561c446..ac6e4b8dc426d1a1fa759c4b8e3328bd9866f430 100644
--- a/ogstools/meshlib/region/region.py
+++ b/ogstools/meshlib/region/region.py
@@ -7,7 +7,7 @@ import numpy as np
import pyvista as pv
from ogs import cli
-from ..boundary_set import LayerSet
+from ..boundary_set import Layer, LayerSet
class RegionSet:
@@ -18,7 +18,7 @@ class RegionSet:
subsets. Each subset within a region is uniquely identified by "MaterialID".
"""
- def __init__(self, input: Union[Path, pv.DataSet]):
+ def __init__(self, input: Union[Path, pv.UnstructuredGrid]):
if type(input) is Path:
self.filename = input
self.mesh = None
@@ -26,7 +26,7 @@ class RegionSet:
self.filename = Path(tempfile.mkstemp(".vtu", "region_set")[1])
self.mesh = input
- def box_boundaries(self):
+ def box_boundaries(self) -> tuple[pv.UnstructuredGrid, ...]:
"""
Retrieve the boundaries of the mesh in local coordinate system (u, v, w).
@@ -46,6 +46,7 @@ class RegionSet:
mesh = ...
u_min, u_max, v_min, v_max, w_min, w_max = mesh.box_boundaries()
"""
+ assert isinstance(self.mesh, pv.UnstructuredGrid)
surface = self.mesh.extract_surface()
u_max = to_boundary(surface, lambda normals: normals[:, 0] > 0.5)
u_min = to_boundary(surface, lambda normals: normals[:, 0] < -0.5)
@@ -152,7 +153,7 @@ def to_region_prism(layer_set: LayerSet, resolution: float) -> RegionSet:
def layer_to_simplified_mesh(
- layer, resolution: float, rank: int, bounds: list[float]
+ layer: Layer, resolution: float, rank: int, bounds: list[float]
) -> pv.UnstructuredGrid:
"""Convert a geological layer to a simplified mesh.
diff --git a/ogstools/meshlib/xdmf_reader.py b/ogstools/meshlib/xdmf_reader.py
index 6f1b49b5797a73bc50b140484ae44eac20611adb..c82d490410affd5557017979d7180e3f09b68abc 100644
--- a/ogstools/meshlib/xdmf_reader.py
+++ b/ogstools/meshlib/xdmf_reader.py
@@ -15,6 +15,8 @@ to be read like::
"""
+from xml.etree.ElementTree import Element
+
import meshio
import numpy as np
from meshio.xdmf.time_series import (
@@ -25,10 +27,10 @@ from meshio.xdmf.time_series import (
class XDMFReader(meshio.xdmf.TimeSeriesReader):
- def __init__(self, filename):
+ def __init__(self, filename: str):
super().__init__(filename)
- def read_data(self, k: int):
+ def read_data(self, k: int) -> tuple[float, dict, dict, dict]:
point_data = {}
cell_data_raw = {}
other_data = {}
@@ -67,8 +69,8 @@ class XDMFReader(meshio.xdmf.TimeSeriesReader):
return t, point_data, cell_data, other_data
- def _read_data_item(self, data_item):
- dims = [int(d) for d in data_item.get("Dimensions").split()]
+ def _read_data_item(self, data_item: Element) -> np.ndarray:
+ dims = [int(d) for d in data_item.get("Dimensions", "").split()]
# Actually, `NumberType` is XDMF2 and `DataType` XDMF3, but many files out there
# use both keys interchangeably.
@@ -92,6 +94,7 @@ class XDMFReader(meshio.xdmf.TimeSeriesReader):
data_format = data_item.attrib["Format"]
+ assert isinstance(data_item.text, str)
if data_format == "XML":
return np.fromstring(
data_item.text,
diff --git a/ogstools/meshplotlib/__init__.py b/ogstools/meshplotlib/__init__.py
index 17498bba65329635341b4dfafbe50b3264caa8dc..74b7c05e1fe9335136b564962ef44591fbab2e0a 100644
--- a/ogstools/meshplotlib/__init__.py
+++ b/ogstools/meshplotlib/__init__.py
@@ -4,14 +4,12 @@
from .plot_setup import _setup as setup # noqa: I001: noqa
from .core import (
- plot_limit,
plot_probe,
plot,
subplot,
update_font_sizes,
label_spatial_axes,
clear_labels,
- clear_labels_ax,
)
from .plot_features import plot_contour, plot_on_top
@@ -19,12 +17,10 @@ __all__ = [
"setup",
"plot",
"plot_contour",
- "plot_limit",
"plot_on_top",
"plot_probe",
"subplot",
"update_font_sizes",
"label_spatial_axes",
"clear_labels",
- "clear_labels_ax",
]
diff --git a/ogstools/meshplotlib/animation.py b/ogstools/meshplotlib/animation.py
index d210699339c16074efe1444aaa9f6c8c94d9c997..e05db2f6a688efe34f83d49f9f6c8f6af3317a79 100644
--- a/ogstools/meshplotlib/animation.py
+++ b/ogstools/meshplotlib/animation.py
@@ -44,7 +44,7 @@ def animate(
index = np.argmin(np.abs(np.asarray(ts) - i))
fig.axes[-1].remove() # remove colorbar
- for ax in np.ravel(fig.axes):
+ for ax in np.ravel(np.asarray(fig.axes)):
ax.clear()
if titles is not None:
setup.title_center = titles[index]
@@ -54,18 +54,20 @@ def animate(
mesh = mesh_series.read_interp(i, True)
with warnings.catch_warnings():
warnings.simplefilter("ignore")
- fig = _draw_plot(mesh, property, fig=fig)
+ fig = _draw_plot(
+ mesh, property, fig=fig
+ ) # type: ignore[assignment]
_func = partial(animate_func, fig=fig)
return FuncAnimation(
- fig,
- _func,
+ fig, # type: ignore[arg-type]
+ _func, # type: ignore[arg-type]
frames=tqdm(ts),
blit=False,
interval=50,
repeat=False,
- init_func=init,
+ init_func=init, # type: ignore[arg-type]
)
@@ -83,7 +85,7 @@ def save_animation(anim: FuncAnimation, filename: str, fps: int) -> bool:
"-vf pad=ceil(iw/2)*2:ceil(ih/2)*2"
).split(" ")
- writer = None
+ writer: Optional[Union[FFMpegWriter, ImageMagickWriter]] = None
if FFMpegWriter.isAvailable():
writer = FFMpegWriter(fps=fps, codec="libx265", extra_args=codec_args)
filename += ".mp4"
@@ -91,7 +93,7 @@ def save_animation(anim: FuncAnimation, filename: str, fps: int) -> bool:
print("\nffmpeg not available. It is recommended for saving animation.")
filename += ".gif"
if ImageMagickWriter.isAvailable():
- writer = "imagemagick"
+ writer = ImageMagickWriter()
else:
print(
"ImageMagick also not available. Falling back to"
diff --git a/ogstools/meshplotlib/core.py b/ogstools/meshplotlib/core.py
index 691dcc35688de8021101726230a4d659c203b3ef..8fe6ac039f3f922ed383d04c7922ca6c262d72aa 100644
--- a/ogstools/meshplotlib/core.py
+++ b/ogstools/meshplotlib/core.py
@@ -2,11 +2,10 @@
import warnings
from math import nextafter
-from typing import Literal, Optional, Union
+from typing import Any, Literal, Optional, Union
import numpy as np
import pyvista as pv
-from matplotlib import axes as matplax
from matplotlib import cm as mcm
from matplotlib import colormaps, rcParams
from matplotlib import colors as mcolors
@@ -14,6 +13,7 @@ from matplotlib import figure as mfigure
from matplotlib import pyplot as plt
from matplotlib import ticker as mticker
from matplotlib.patches import Rectangle as Rect
+from typeguard import typechecked
from ogstools.meshlib import MeshSeries
from ogstools.propertylib import Property, Vector
@@ -28,13 +28,13 @@ from .utils import get_style_cycler
# TODO: define default data_name for regions in setup
-def _q_zero_line(mesh_property: Property, levels: np.ndarray):
+def _q_zero_line(mesh_property: Property, levels: np.ndarray) -> bool:
return mesh_property.bilinear_cmap or (
mesh_property.data_name == "temperature" and levels[0] < 0 < levels[-1]
)
-def get_level_boundaries(levels: np.ndarray):
+def get_level_boundaries(levels: np.ndarray) -> np.ndarray:
return np.array(
[
levels[0] - 0.5 * (levels[1] - levels[0]),
@@ -57,6 +57,7 @@ def get_cmap_norm(
continuous_cmap = colormaps[mesh_property.cmap]
else:
continuous_cmap = mesh_property.cmap
+ conti_norm: Union[mcolors.TwoSlopeNorm, mcolors.Normalize]
if mesh_property.bilinear_cmap:
if vmin <= 0.0 <= vmax:
vcenter = 0.0
@@ -258,9 +259,13 @@ def subplot(
cmap, norm = get_cmap_norm(levels, mesh_property)
if mesh_property.data_name in mesh.point_data:
- ax.tricontourf(x, y, tri, values, levels=levels, cmap=cmap, norm=norm)
+ ax.tricontourf( # type: ignore[call-overload]
+ x, y, tri, values, levels=levels, cmap=cmap, norm=norm
+ )
if _q_zero_line(mesh_property, levels):
- ax.tricontour(x, y, tri, values, levels=[0], colors="w")
+ ax.tricontour( # type: ignore[call-overload]
+ x, y, tri, values, levels=[0], colors="w"
+ )
else:
ax.tripcolor(x, y, tri, facecolors=values, cmap=cmap, norm=norm)
if mesh_property.is_mask():
@@ -299,77 +304,38 @@ def subplot(
sec_labels += [""]
# TODO: use a function to make this short
secax = ax.secondary_xaxis("top")
- secax.xaxis.set_major_locator(mticker.FixedLocator(ax.get_xticks()))
+ secax.xaxis.set_major_locator(
+ mticker.FixedLocator(list(ax.get_xticks()))
+ )
secax.set_xticklabels(sec_labels)
secax.set_xlabel(f'{"xyz"[projection]} / {setup.length.output_unit}')
-def clear_labels_ax(ax: plt.axes) -> plt.axes:
- ax.set_xlabel(None)
- ax.set_ylabel(None)
- return ax
-
-
-def clear_labels(ax: Union[plt.axes, np.array]) -> Union[plt.axes, np.array]:
- if isinstance(ax, np.ndarray):
- for id_n in range(ax.shape[0]):
- for id_m in range(ax.shape[1]):
- ax[id_n, id_m] = clear_labels_ax(ax[id_n, id_m])
- if isinstance(ax, matplax.Axes):
- # Wrap single axis in np.array
- ax = clear_labels_ax(ax)
- return ax
+def clear_labels(axes: Union[plt.Axes, np.ndarray]) -> None:
+ ax: plt.Axes
+ for ax in np.ravel(np.array(axes)):
+ ax.set_xlabel("")
+ ax.set_ylabel("")
+@typechecked
def label_spatial_axes(
- ax: Union[plt.axes, np.array],
- ax_ids: np.array,
-) -> plt.axes:
+ axes: Union[plt.Axes, np.ndarray], x_label: str = "x", y_label: str = "y"
+) -> None:
"""
- Add labels to X and Y axis
-
- Automatically selects correct pair of directions and unit.
- Respects sharex / sharey settings.
+ Add labels to x and y axis.
- :param ax: Matplotlib Axis object
- :param ax_ids: indices of axes label [0,1,2] for [x,y,z] for horizontal and vertical axis
+ If given an array of axes, only the outer axes will be labeled.
"""
- if isinstance(ax, np.ndarray):
- # Labels will be applied to shared axis
- # Shared axis = value in projection is not None
- if ax_ids[0] is not None:
- x_label = (
- setup.x_label
- or f'{"xyz"[ax_ids[0]]} / {setup.length.output_unit}'
- )
- for ax_temp in ax[-1, :]:
- ax_temp.set_xlabel(x_label)
- if ax_ids[1] is not None:
- y_label = (
- setup.y_label
- or f'{"xyz"[ax_ids[1]]} / {setup.length.output_unit}'
- )
- for ax_temp in ax[:, 0]:
- ax_temp.set_ylabel(y_label)
- elif isinstance(ax, matplax.Axes):
- # Labels will only be applied to non shared axis:
- # Non shared = value in projection is not None
- if ax_ids[0] is not None:
- x_label = (
- setup.x_label
- or f'{"xyz"[ax_ids[0]]} / {setup.length.output_unit}'
- )
- ax.set_xlabel(x_label)
- if ax_ids[1] is not None:
- y_label = (
- setup.y_label
- or f'{"xyz"[ax_ids[1]]} / {setup.length.output_unit}'
- )
- ax.set_ylabel(y_label)
+ if isinstance(axes, np.ndarray):
+ ax: plt.Axes
+ for ax in axes[-1, :]:
+ ax.set_xlabel(f"{x_label} / {setup.length.output_unit}")
+ for ax in axes[:, 0]:
+ ax.set_ylabel(f"{y_label} / {setup.length.output_unit}")
else:
- msg = "ax is neither Matplotlib axis nor Numpy array!"
- raise TypeError(msg)
- return ax
+ axes.set_xlabel(f"{x_label} / {setup.length.output_unit}")
+ axes.set_ylabel(f"{y_label} / {setup.length.output_unit}")
def _get_rows_cols(
@@ -428,7 +394,7 @@ def get_combined_levels(
"""
Calculate well spaced levels for the encompassing property range in meshes.
"""
- mesh_property = get_preset(mesh_property, meshes[0])
+ mesh_property = get_preset(mesh_property, meshes.ravel()[0])
p_min, p_max = np.inf, -np.inf
unique_vals = np.array([])
for mesh in np.ravel(meshes):
@@ -460,21 +426,21 @@ def _draw_plot(
meshes: Union[list[pv.UnstructuredGrid], np.ndarray, pv.UnstructuredGrid],
mesh_property: Property,
fig: Optional[mfigure.Figure] = None,
- ax: Optional[plt.Axes] = None,
-) -> mfigure.Figure:
+ axes: Optional[plt.Axes] = None,
+) -> Optional[mfigure.Figure]:
"""
Plot the property field of meshes on existing figure.
:param meshes: Singular mesh of 2D numpy array of meshes
:param property: the property field to be visualized on all meshes
- :param fig: Matplotlib Figure object to use for plotting (optional)
- :param ax: Matplotlib Axis object to use for plotting (optional)
+ :param fig: Matplotlib figure to use for plotting (optional)
+ :param axes: Matplotlib Axes to use for plotting (optional)
"""
shape = _get_rows_cols(meshes)
np_meshes = np.reshape(meshes, shape)
- if fig is not None and ax is not None:
- np_axs = np.reshape(np.array(ax), shape)
- elif fig is not None and ax is None:
+ if fig is not None and axes is not None:
+ np_axs = np.reshape(np.array(axes), shape)
+ elif fig is not None and axes is None:
# Only Fig is given
# Multiple meshes should be accepted
warnings.warn(
@@ -482,14 +448,14 @@ def _draw_plot(
Warning,
stacklevel=4,
)
- np_axs = np.reshape(fig.axes, shape)
- elif fig is None and ax is not None:
+ np_axs = np.reshape(np.asarray(fig.axes), shape)
+ elif fig is None and axes is not None:
# Only ax is given
# Only one mesh should be accepted
if shape != (1, 1):
msg = "You have provided only one Axis object but multiple meshes. Provide only one mesh per Axis object, or provide Figure object instead."
raise ValueError(msg)
- np_axs = np.reshape(np.array(ax), (1, 1))
+ np_axs = np.reshape(np.array(axes), (1, 1))
else:
msg = "Neither Figure nor Axis object was provided."
raise TypeError(msg)
@@ -500,14 +466,14 @@ def _draw_plot(
_levels = (
combined_levels
if setup.combined_colorbar
- else get_combined_levels([np_meshes[i, j]], mesh_property)
+ else get_combined_levels(np_meshes[i, j, None], mesh_property)
)
subplot(np_meshes[i, j], mesh_property, np_axs[i, j], _levels)
x_id, y_id = get_projection(
np_meshes[0, 0]
) # One mesh is sufficient, it should be the same for all of them
- np_axs = label_spatial_axes(np_axs, np.array([x_id, y_id]))
+ label_spatial_axes(np_axs, "xyz"[x_id], "xyz"[y_id])
np_axs[0, 0].set_title(setup.title_center, loc="center", y=1.02)
np_axs[0, 0].set_title(setup.title_left, loc="left", y=1.02)
np_axs[0, 0].set_title(setup.title_right, loc="right", y=1.02)
@@ -522,7 +488,7 @@ def _draw_plot(
stacklevel=4,
)
else:
- cb_axs = np.ravel(fig.axes).tolist()
+ cb_axs = np.ravel(np.asarray(fig.axes)).tolist()
add_colorbars(
fig, cb_axs, mesh_property, combined_levels, pad=0.05 / shape[1]
)
@@ -538,13 +504,13 @@ def _draw_plot(
for i in range(shape[0]):
for j in range(shape[1]):
_levels = get_combined_levels(
- [np_meshes[i, j]], mesh_property
+ np_meshes[i, j, None], mesh_property
)
add_colorbars(fig, np_axs[i, j], mesh_property, _levels)
return fig
-def get_data_aspect(mesh: pv.DataSet) -> float:
+def get_data_aspect(mesh: pv.UnstructuredGrid) -> float:
"""
Calculate the data aspect ratio of a 2D mesh.
"""
@@ -573,10 +539,14 @@ def update_font_sizes(
subax_xlim = subax.get_xlim()
subax_ylim = subax.get_ylim()
subax.set_xticks(
- subax.get_xticks(), subax.get_xticklabels(), fontsize=fontsize
+ subax.get_xticks(),
+ [label.get_text() for label in subax.get_xticklabels()],
+ fontsize=fontsize,
)
subax.set_yticks(
- subax.get_yticks(), subax.get_yticklabels(), fontsize=fontsize
+ subax.get_yticks(),
+ [label.get_text() for label in subax.get_yticklabels()],
+ fontsize=fontsize,
)
subax.set_xlim(subax_xlim)
subax.set_ylim(subax_ylim)
@@ -590,7 +560,7 @@ def plot(
mesh_property: Union[Property, str],
fig: Optional[mfigure.Figure] = None,
ax: Optional[plt.Axes] = None,
-) -> mfigure.Figure:
+) -> Optional[mfigure.Figure]:
"""
Plot the property field of meshes with default settings.
@@ -604,7 +574,7 @@ def plot(
"""
rcParams.update(setup.rcParams_scaled)
shape = _get_rows_cols(meshes)
- _meshes = np.reshape(meshes, shape).flatten()
+ _meshes = np.reshape(meshes, shape).ravel()
mesh_property = get_preset(mesh_property, _meshes[0])
data_aspects = np.asarray([get_data_aspect(mesh) for mesh in _meshes])
if setup.min_ax_aspect is None and setup.max_ax_aspect is None:
@@ -617,52 +587,25 @@ def plot(
n_axs = shape[0] * shape[1]
if ax is None and fig is None:
_fig, _ax = _fig_init(rows=shape[0], cols=shape[1], aspect=fig_aspect)
- fig = _draw_plot(meshes, mesh_property, fig=_fig, ax=_ax)
+ fig = _draw_plot(meshes, mesh_property, fig=_fig, axes=_ax)
+ assert isinstance(fig, plt.Figure)
for ax, aspect in zip(fig.axes[: n_axs + 1], ax_aspects):
ax.set_aspect(1.0 / aspect)
elif ax is not None and fig is None:
- _draw_plot(meshes, mesh_property, ax=ax)
+ _draw_plot(meshes, mesh_property, axes=ax)
ax.set_aspect(1.0 / ax_aspects[0])
elif ax is None and fig is not None:
fig = _draw_plot(meshes, mesh_property, fig=fig)
+ assert isinstance(fig, plt.Figure)
for ax, aspect in zip(fig.axes[: n_axs + 1], ax_aspects):
ax.set_aspect(1.0 / aspect)
elif ax is not None and fig is not None:
- _draw_plot(meshes, mesh_property, fig=fig, ax=ax)
+ _draw_plot(meshes, mesh_property, fig=fig, axes=ax)
for ax, aspect in zip(fig.axes[: n_axs + 1], ax_aspects):
ax.set_aspect(1.0 / aspect)
return fig
-def plot_limit(
- mesh_series: MeshSeries,
- mesh_property: Union[Property, str],
- limit: Literal["min", "max"],
- fig: Optional[mfigure.Figure] = None,
- ax: Optional[plt.Axes] = None,
-) -> mfigure.Figure:
- """
- Plot the property limits through all timesteps of a MeshSeries.
-
- :param mesh_series: MeshSeries object containing the data to be plotted
- :param property: The property field to be evaluated
- :param limit: Type of limit to be computed
- :param fig: Matplotlib Figure object to use for plotting (optional)
- :param ax: Matplotlib Axis object to use for plotting (optional)
-
- :returns: A matplotlib Figure
- """
- mesh = mesh_series.read(0)
- mesh_property = get_preset(mesh_property, mesh)
- func = {"min": np.min, "max": np.max}[limit]
- vals = mesh_series.values(mesh_property.data_name)
- func(vals, out=mesh[mesh_property.data_name], axis=0)
- limit_property = mesh_property.replace(
- output_name=limit + " " + mesh_property.output_name
- )
- return plot(mesh, limit_property, fig=fig, ax=ax)
-
-
def plot_probe(
mesh_series: MeshSeries,
points: np.ndarray,
@@ -676,8 +619,8 @@ def plot_probe(
linestyles: Optional[list] = None,
ax: Optional[plt.Axes] = None,
fill_between: bool = False,
- **kwargs,
-) -> mfigure.Figure:
+ **kwargs: Any,
+) -> Optional[mfigure.Figure]:
"""
Plot the transient property on the observation points in the MeshSeries.
diff --git a/ogstools/meshplotlib/plot_features.py b/ogstools/meshplotlib/plot_features.py
index 5b358b30cf2f743385bfb60bedc67690e3b53386..9a5ee911c2c4d6b11a90a506e3989e6b52b5c343 100644
--- a/ogstools/meshplotlib/plot_features.py
+++ b/ogstools/meshplotlib/plot_features.py
@@ -74,7 +74,9 @@ def plot_element_edges(ax: plt.Axes, surf: pv.DataSet, projection: int) -> None:
]
verts = setup.length.transform(np.delete(cell_pts, projection, -1))
lw = 0.5 * setup.rcParams_scaled["lines.linewidth"]
- pc = PolyCollection(verts, fc="None", ec="black", lw=lw)
+ pc = PolyCollection(
+ verts, fc="None", ec="black", lw=lw # type: ignore[arg-type]
+ )
ax.add_collection(pc)
@@ -131,9 +133,9 @@ def plot_streamlines(
lw = 2.5 * val_norm / max(1e-16, np.max(val_norm))
lw *= setup.rcParams_scaled["lines.linewidth"]
x_g, y_g = setup.length.transform(np.meshgrid(x, y))
- plot_args = [x_g, y_g, val[..., 0], val[..., 1]]
if plot_type == "streamlines":
- ax.streamplot(*plot_args, color="k", linewidth=lw, density=1.5)
+ ax.streamplot(x_g, y_g, val[..., 0], val[..., 1],
+ color="k", linewidth=lw, density=1.5) # fmt: skip
else:
line_args = (
dict( # noqa: C408
@@ -142,7 +144,8 @@ def plot_streamlines(
if plot_type == "lines"
else {}
)
- ax.quiver(*plot_args, **line_args, scale=1 / 0.03)
+ scale = 1.0 / 0.03
+ ax.quiver(x_g, y_g, val[..., 0], val[..., 1], **line_args, scale=scale)
def plot_on_top(
@@ -161,7 +164,7 @@ def plot_on_top(
x_vals = df_pts.groupby("x")["x"].agg(np.mean).to_numpy()
y_vals = df_pts.groupby("x")["y"].agg(np.max).to_numpy()
contour_vals = [y + scaling * contour(x) for y, x in zip(y_vals, x_vals)]
- ax.set_ylim(top=setup.length.transform(np.max(contour_vals)))
+ ax.set_ylim(top=float(setup.length.transform(np.max(contour_vals))))
ax.fill_between(
setup.length.transform(x_vals),
setup.length.transform(y_vals),
@@ -172,7 +175,7 @@ def plot_on_top(
def plot_contour(
ax: plt.Axes, mesh: pv.DataSet, style: str, lw: int, projection: int = 2
-):
+) -> None:
contour = mesh.extract_surface().strip(join=True)
x_id, y_id = np.delete([0, 1, 2], projection)
x, y = 1e-3 * contour.points[contour.lines[1:]].T[[x_id, y_id]]
diff --git a/ogstools/meshplotlib/plot_setup.py b/ogstools/meshplotlib/plot_setup.py
index d36e7f327061557b55cabcaa6cdfdbb35a62acac..6c933c6b779641275be17059ae374aefb6b8d255 100644
--- a/ogstools/meshplotlib/plot_setup.py
+++ b/ogstools/meshplotlib/plot_setup.py
@@ -1,7 +1,7 @@
"""Plot configuration setup."""
from dataclasses import dataclass
-from typing import Union
+from typing import Optional, Union
from ogstools.propertylib.property import Scalar
@@ -23,9 +23,9 @@ class PlotSetup:
"The resolution (dots per inch) for the figure."
fig_scale: float
"A scaling factor for the figure."
- min_ax_aspect: float
+ min_ax_aspect: Optional[float]
"Minimum aspect ratio of subplots."
- max_ax_aspect: float
+ max_ax_aspect: Optional[float]
"Maximum aspect ratio of subplots."
invert_colorbar: bool
"A boolean indicating whether to invert the colorbar."
@@ -38,11 +38,11 @@ class PlotSetup:
num_levels: int
"""The aimed number of levels / bins of the colorbar. See
:obj:`ogstools.meshplotlib.levels`"""
- num_streamline_interp_pts: int
+ num_streamline_interp_pts: Optional[int]
"The number of interpolation points for streamlines."
- p_max: float
+ p_max: Optional[float]
"The fixed upper limit for the current scale."
- p_min: float
+ p_min: Optional[float]
"The fixed lower limit for the current scale."
rcParams: dict
"""Matplotlib runtime configuration. See
@@ -78,7 +78,7 @@ class PlotSetup:
return params
@classmethod
- def from_dict(cls: type["PlotSetup"], obj: dict):
+ def from_dict(cls: type["PlotSetup"], obj: dict) -> "PlotSetup":
"""Create a PlotSetup instance from a dictionary."""
return cls(
fig_scale=obj["fig_scale"],
diff --git a/ogstools/meshplotlib/utils.py b/ogstools/meshplotlib/utils.py
index e31ce760c94b50f10df3f2ac01913d0ebc62718d..8cdc448c9cd4a69727ace1a78db66ab077766488 100644
--- a/ogstools/meshplotlib/utils.py
+++ b/ogstools/meshplotlib/utils.py
@@ -2,12 +2,13 @@ from typing import Optional
import matplotlib.pyplot as plt
import numpy as np
+from cycler import Cycler
def justified_labels(points: np.ndarray) -> list[str]:
"Formats an array of points to a list of aligned str."
- def fmt(val: float):
+ def fmt(val: float) -> str:
return f"{val:.2f}".rstrip("0").rstrip(".")
col_lens = np.max(
@@ -24,7 +25,7 @@ def get_style_cycler(
min_number_of_styles: int,
colors: Optional[Optional[list]] = None,
linestyles: Optional[list] = None,
-) -> plt.cycler:
+) -> Cycler:
if colors is None:
colors = plt.rcParams["axes.prop_cycle"].by_key()["color"]
if linestyles is None:
diff --git a/ogstools/msh2vtu/__init__.py b/ogstools/msh2vtu/__init__.py
index ac2a24d6e224523cd716e46e8f5d53183f802288..251124ed3e6b399185a6d2244a504e16105d2087 100644
--- a/ogstools/msh2vtu/__init__.py
+++ b/ogstools/msh2vtu/__init__.py
@@ -1,7 +1,7 @@
# Author: Dominik Kern (TU Bergakademie Freiberg)
import logging
from pathlib import Path
-from typing import Union
+from typing import Any, Union
import meshio
import numpy as np
@@ -18,7 +18,7 @@ __all__ = [
]
-def my_remove_orphaned_nodes(my_mesh: meshio.Mesh):
+def my_remove_orphaned_nodes(my_mesh: meshio.Mesh) -> None:
"""Auxiliary function to remove points not belonging to any cell"""
# find connected points and derive mapping from all points to them
@@ -51,11 +51,9 @@ def my_remove_orphaned_nodes(my_mesh: meshio.Mesh):
# cell data are not affected by point changes
my_mesh.cells = output_cell_blocks # update cells
- return
-
# print info for mesh: statistics and data field names
-def print_info(mesh: meshio.Mesh):
+def print_info(mesh: meshio.Mesh) -> None:
N, D = mesh.points.shape
logging.info("%d points in %d dimensions", N, D)
cell_info = "cells: "
@@ -70,13 +68,15 @@ def print_info(mesh: meshio.Mesh):
# function to create node connectivity list, i.e. store for each node (point) to
# which element (cell) it belongs
-def find_cells_at_nodes(cells, node_count, cell_start_index):
+def find_cells_at_nodes(
+ cells: Any, node_count: int, cell_start_index: int
+) -> list[set]:
# depending on the numbering of mixed meshes in OGS one may think of an
# object-oriented way to add elements (of different type) to node
# connectivity
# initialize list of sets
- node_connectivity = [set() for _ in range(node_count)]
+ node_connectivity: list[set] = [set() for _ in range(node_count)]
cell_index = cell_start_index
for cell in cells:
for node in cell:
@@ -95,7 +95,7 @@ def find_cells_at_nodes(cells, node_count, cell_start_index):
# function to find out to which domain elements a boundary element belongs
def find_connected_domain_cells(
- boundary_cells_values, domain_cells_at_node
+ boundary_cells_values: Any, domain_cells_at_node: list[set[int]]
) -> tuple[np.ndarray, np.ndarray]:
warned_gt1 = False # to avoid flood of warnings
warned_lt1 = False # to avoid flood of warnings
@@ -160,7 +160,7 @@ def msh2vtu(
keep_ids: bool = False,
ascii: bool = False,
log_level: Union[int, str] = "DEBUG",
-):
+) -> int:
"""
Convert a gmsh mesh (.msh) to an unstructured grid file (.vtu).
diff --git a/ogstools/propertylib/matrix.py b/ogstools/propertylib/matrix.py
index e37c88117119b9fc5730d910943239323244fa44..96b3a9a9eb5ce1ffa44a04657f4266eb669647c9 100644
--- a/ogstools/propertylib/matrix.py
+++ b/ogstools/propertylib/matrix.py
@@ -8,7 +8,7 @@ from ogstools.propertylib.vector import Vector, VectorList
@dataclass
class Matrix(Property):
- """Represent a matrix property of a dataset.
+ """Represent a matrix property.
Matrix properties should contain either 4 (2D) or 6 (3D) components.
Matrix components can be accesses with brackets e.g. stress[0]
diff --git a/ogstools/propertylib/mesh_dependent.py b/ogstools/propertylib/mesh_dependent.py
index 373fa25a0245721ba4cd09f1996628e67d1dd2a8..64bf6324f31d500cea64a61d7daea98229bb1ba5 100644
--- a/ogstools/propertylib/mesh_dependent.py
+++ b/ogstools/propertylib/mesh_dependent.py
@@ -7,7 +7,7 @@ import pyvista as pv
from pint.facets.plain import PlainQuantity
from .property import Property
-from .tensor_math import eigenvalues, mean, octahedral_shear
+from .tensor_math import _split_quantity, eigenvalues, mean, octahedral_shear
from .unit_registry import u_reg
ValType = Union[PlainQuantity, np.ndarray]
@@ -48,11 +48,16 @@ def depth(mesh: pv.UnstructuredGrid, use_coords: bool = False) -> np.ndarray:
]
edge_centers = edges.cell_centers().points
adj_cells = [mesh.find_closest_cell(point) for point in edge_centers]
- adj_cell_centers = mesh.extract_cells(adj_cells).cell_centers().points
- are_above = (
- edge_centers[..., vertical_dim]
- > adj_cell_centers[..., vertical_dim]
+ adj_centers = np.vstack(
+ [
+ mesh.extract_cells(adj_cell).cell_centers().points
+ for adj_cell in adj_cells
+ ]
)
+ are_above = [
+ edge_center[vertical_dim] > adj_center[vertical_dim]
+ for edge_center, adj_center in zip(edge_centers, adj_centers)
+ ]
are_non_vertical = np.asarray(edge_horizontal_extent) > 1e-12
top_cells = are_above & are_non_vertical
top = edges.extract_cells(top_cells)
@@ -110,8 +115,9 @@ def fluid_pressure_criterion(
"""
Qty = u_reg.Quantity
- sigma = -Qty(mesh[mesh_property.data_name], mesh_property.data_unit)
- return p_fluid(mesh) - eigenvalues(sigma)[..., 0]
+ sigma = mesh[mesh_property.data_name]
+ sig_min = _split_quantity(eigenvalues(-sigma))[0][..., 0]
+ return p_fluid(mesh) - Qty(sig_min, mesh_property.data_unit)
def dilatancy_critescu(
diff --git a/ogstools/propertylib/presets.py b/ogstools/propertylib/presets.py
index 7d07843b4d497e4be4fc746f0dce93c498354398..9da9a7f0262736a9c91fd40ca1abd85a3b47c7ac 100644
--- a/ogstools/propertylib/presets.py
+++ b/ogstools/propertylib/presets.py
@@ -14,6 +14,7 @@ from . import mesh_dependent, tensor_math
from .custom_colormaps import integrity_cmap, temperature_cmap
from .matrix import Matrix
from .property import Property, Scalar
+from .tensor_math import identity
from .vector import Vector
T_MASK = "temperature_active"
@@ -142,21 +143,45 @@ all_properties = [v for v in locals().values() if isinstance(v, Property)]
def get_preset(
- mesh_property: Union[Property, str], mesh: pv.DataSet
+ mesh_property: Union[Property, str], mesh: pv.UnstructuredGrid
) -> Property:
"""
Returns a Property preset or creates one with correct type.
Searches for presets by data_name and output_name and returns if found.
+ If 'mesh_property' is given as type Property this will also look for
+ derived properties (difference, aggregate).
Otherwise create Scalar, Vector, or Matrix Property depending on the shape
of data in mesh.
:param mesh_property: The property to retrieve or its name if a string.
- :param mesh: The mesh dataset containing the property data.
+ :param mesh: The mesh containing the property data.
:returns: A corresponding Property preset or a new Property of correct type.
"""
+ data_keys: list[str] = list(set().union(mesh.point_data, mesh.cell_data))
+ error_msg = (
+ f"Data not found in mesh. Available data names are {data_keys}. "
+ )
+
if isinstance(mesh_property, Property):
- return mesh_property
+ if mesh_property.data_name in data_keys:
+ return mesh_property
+ matches = [
+ mesh_property.output_name in data_key for data_key in data_keys
+ ]
+ if not any(matches):
+ raise KeyError(error_msg)
+ data_key = data_keys[matches.index(True)]
+ if data_key == f"{mesh_property.output_name}_difference":
+ return mesh_property.difference
+ return mesh_property.replace(
+ data_name=data_key,
+ data_unit=mesh_property.output_unit,
+ output_unit=mesh_property.output_unit,
+ output_name=data_key,
+ func=identity,
+ mesh_dependent=False,
+ )
for prop in all_properties:
if prop.output_name == mesh_property:
@@ -164,9 +189,11 @@ def get_preset(
for prop in all_properties:
if prop.data_name == mesh_property:
return prop
- if mesh_property not in set().union(mesh.point_data, mesh.cell_data):
- msg = f"Property {mesh_property} not found in mesh."
- raise KeyError(msg)
+
+ matches = [mesh_property in data_key for data_key in data_keys]
+ if not any(matches):
+ raise KeyError(error_msg)
+
data_shape = mesh[mesh_property].shape
if len(data_shape) == 1:
return Scalar(mesh_property)
diff --git a/ogstools/propertylib/property.py b/ogstools/propertylib/property.py
index 50a9e84f027b42637b29265539bda5fee85429ef..b77eef48e26e77171ffb6505c93c69f87bd1f8bf 100644
--- a/ogstools/propertylib/property.py
+++ b/ogstools/propertylib/property.py
@@ -7,11 +7,12 @@ via pint.
from collections.abc import Sequence
from dataclasses import dataclass, replace
-from typing import Callable, Union
+from typing import Any, Callable, Union
import numpy as np
import pyvista as pv
from matplotlib.colors import Colormap
+from pint.facets.plain import PlainQuantity
from .custom_colormaps import mask_cmap
from .tensor_math import identity
@@ -20,18 +21,18 @@ from .unit_registry import u_reg
@dataclass
class Property:
- """Represent a property of a dataset."""
+ """Represent a generic mesh property."""
data_name: str
- """The name of the property data in the dataset."""
+ """The name of the property data in the mesh."""
data_unit: str = ""
- """The unit of the property data in the dataset."""
+ """The unit of the property data in the mesh."""
output_unit: str = ""
"""The output unit of the property."""
output_name: str = ""
"""The output name of the property."""
mask: str = ""
- """The name of the mask data in the dataset."""
+ """The name of the mask data in the mesh."""
func: Callable = identity
"""The function to be applied on the data.
.. seealso:: :meth:`~ogstools.propertylib.Property.transform`"""
@@ -46,15 +47,15 @@ class Property:
categoric: bool = False
"""Does this property only have categoric values?"""
- def __post_init__(self):
+ def __post_init__(self) -> None:
if not self.output_name:
self.output_name = self.data_name
@property
- def type_name(self):
+ def type_name(self) -> str:
return type(self).__name__
- def replace(self, **changes):
+ def replace(self: "Property", **changes: Any) -> "Property":
"""
Create a new Property object with modified attributes.
@@ -68,7 +69,9 @@ class Property:
return replace(self, **changes)
@classmethod
- def from_property(cls, new_property: "Property", **changes):
+ def from_property( # type: ignore[no-untyped-def]
+ cls, new_property: "Property", **changes: Any
+ ):
"Create a new Property object with modified attributes."
return cls(
data_name=new_property.data_name,
@@ -85,31 +88,31 @@ class Property:
def transform(
self,
- data: Union[int, float, np.ndarray, pv.DataSet, Sequence],
+ data: Union[int, float, np.ndarray, pv.UnstructuredGrid, Sequence],
strip_unit: bool = True,
) -> np.ndarray:
"""
Return the transformed data values.
- Converts the data from data_unit to output_unit and apply a function
- the transformation function of this property. The result is returned by
+ Converts the data from data_unit to output_unit and applies the
+ transformation function of this property. The result is returned by
default without units. if `strip_unit` is False, a quantity is returned.
Note:
If `self.mesh_dependent` is True, `self.func` is applied directly to the
- DataSet. Otherwise, it is determined by `self.process_with_units` if the
+ mesh. Otherwise, it is determined by `self.process_with_units` if the
data is passed to the function with units (i.e. as a pint quantity) or
without.
"""
Qty, d_u, o_u = u_reg.Quantity, self.data_unit, self.output_unit
if self.mesh_dependent:
- if isinstance(data, pv.DataSet):
+ if isinstance(data, (pv.DataSet, pv.UnstructuredGrid)):
result = Qty(self.func(data, self), o_u)
else:
msg = "This property can only be evaluated on a mesh."
raise TypeError(msg)
else:
- if isinstance(data, pv.DataSet):
+ if isinstance(data, (pv.DataSet, pv.UnstructuredGrid)):
result = Qty(self.func(Qty(self._get_data(data), d_u)), o_u)
elif self.process_with_units:
result = Qty(self.func(Qty(data, d_u)), o_u)
@@ -126,20 +129,22 @@ class Property:
return "%" if self.output_unit == "percent" else self.output_unit
@property
- def delta(self) -> "Property":
- "returns: A property relating to the difference in a quantity."
- data_property = self.replace(output_unit=self.data_unit)
- diff_unit = str(
- (
- data_property.transform(1, strip_unit=False)
- - data_property.transform(1, strip_unit=False)
- ).units
- )
+ def difference(self) -> "Property":
+ "returns: A property relating to differences in this quantity."
+ quantity = u_reg.Quantity(1, self.output_unit)
+ diff_quantity: PlainQuantity = quantity - quantity
+ diff_unit = str(diff_quantity.units)
+ if diff_unit == "delta_degC":
+ diff_unit = "kelvin"
+ outname = self.output_name + "_difference"
return self.replace(
+ data_name=outname,
data_unit=diff_unit,
output_unit=diff_unit,
- output_name=self.output_name + " difference",
+ output_name=outname,
bilinear_cmap=True,
+ func=identity,
+ mesh_dependent=False,
cmap=self.cmap if self.bilinear_cmap else "coolwarm",
)
@@ -151,7 +156,7 @@ class Property:
"""
return self.data_name == self.mask
- def get_mask(self):
+ def get_mask(self) -> "Property":
"""
:returns: A property representing this properties mask.
"""
@@ -174,9 +179,14 @@ class Property:
def _get_data(
self, mesh: pv.UnstructuredGrid, masked: bool = True
) -> np.ndarray:
- """Get the data associated with a scalar or vector property from a mesh."""
- if self.data_name not in set().union(mesh.point_data, mesh.cell_data):
- msg = f"Data name {self.data_name} not found in mesh."
+ "Get the data associated with a scalar or vector property from a mesh."
+ if self.data_name not in (
+ data_keys := ",".join(set().union(mesh.point_data, mesh.cell_data))
+ ):
+ msg = (
+ f"Data name {self.data_name} not found in mesh. "
+ f"Available data names are {data_keys}. "
+ )
raise KeyError(msg)
if masked and self.mask_used(mesh):
return mesh.ctp(True).threshold(value=[1, 1], scalars=self.mask)[
@@ -194,4 +204,4 @@ class Property:
@dataclass
class Scalar(Property):
- "Represent a scalar property of a dataset."
+ "Represent a scalar property."
diff --git a/ogstools/propertylib/vector.py b/ogstools/propertylib/vector.py
index 25abd336fa89bfe8f211626a6b98a9f9e39cc0f6..e7c0630a9198f805f4fe8f2b75c1e8fee2ec9d37 100644
--- a/ogstools/propertylib/vector.py
+++ b/ogstools/propertylib/vector.py
@@ -6,25 +6,21 @@ from pint.facets.plain import PlainQuantity
from ogstools.propertylib.property import Property, Scalar
-from .unit_registry import u_reg
+from .tensor_math import _split_quantity, _to_quantity
ValType = Union[PlainQuantity, np.ndarray]
-def vector_norm(vals: ValType) -> ValType:
+def vector_norm(values: ValType) -> ValType:
":returns: The norm of the vector."
- if isinstance(vals, PlainQuantity):
- unit = vals.units
- vals = vals.magnitude
- else:
- unit = None
+ vals, unit = _split_quantity(values)
result = np.linalg.norm(vals, axis=-1)
- return result if unit is None else u_reg.Quantity(result, unit)
+ return _to_quantity(result, unit)
@dataclass
class Vector(Property):
- """Represent a vector property of a dataset.
+ """Represent a vector property.
Vector properties should contain either 2 (2D) or 3 (3D) components.
Vector components can be accesses with brackets e.g. displacement[0]
@@ -58,7 +54,7 @@ class Vector(Property):
@dataclass
class VectorList(Property):
- """Represent a list of vector properties of a dataset."""
+ """Represent a list of vector properties."""
def __getitem__(self, index: int) -> Vector:
":returns: A vector property as a component of the vectorlist property."
diff --git a/ogstools/studies/convergence/convergence.py b/ogstools/studies/convergence/convergence.py
index da1ea4c23bb86956a4d86f29ad1df93889b99c07..c5282ecfa3d80a6dd0cc88b29cca46c071d790d0 100644
--- a/ogstools/studies/convergence/convergence.py
+++ b/ogstools/studies/convergence/convergence.py
@@ -15,8 +15,8 @@ u_reg.default_format = "~.3g"
def resample(
- topology: pv.DataSet, meshes: list[pv.DataSet]
-) -> list[pv.DataSet]:
+ topology: pv.UnstructuredGrid, meshes: list[pv.UnstructuredGrid]
+) -> list[pv.UnstructuredGrid]:
meshes_resampled = []
for mesh in meshes:
mesh_temp = deepcopy(topology)
@@ -27,7 +27,7 @@ def resample(
return meshes_resampled
-def add_grid_spacing(mesh: pv.DataSet) -> pv.DataSet:
+def add_grid_spacing(mesh: pv.UnstructuredGrid) -> pv.UnstructuredGrid:
dim = mesh.get_cell(0).dimension
key = ["Length", "Area", "Volume"][dim - 1]
_mesh = mesh.compute_cell_sizes()
@@ -36,11 +36,11 @@ def add_grid_spacing(mesh: pv.DataSet) -> pv.DataSet:
def grid_convergence(
- meshes: list[pv.DataSet],
+ meshes: list[pv.UnstructuredGrid],
mesh_property: propertylib.Property,
- topology: pv.DataSet,
+ topology: pv.UnstructuredGrid,
refinement_ratio: float,
-) -> pv.DataSet:
+) -> pv.UnstructuredGrid:
"""
Calculate the grid convergence field for the given meshes on the topology.
@@ -92,11 +92,11 @@ def grid_convergence(
def richardson_extrapolation(
- meshes: list[pv.DataSet],
+ meshes: list[pv.UnstructuredGrid],
mesh_property: propertylib.Property,
- topology: pv.DataSet,
+ topology: pv.UnstructuredGrid,
refinement_ratio: float,
-) -> pv.DataSet:
+) -> pv.UnstructuredGrid:
"""
Estimate a better approximation of a property on a mesh.
@@ -133,8 +133,8 @@ def richardson_extrapolation(
def convergence_metrics(
- meshes: list[pv.DataSet],
- reference: pv.DataSet,
+ meshes: list[pv.UnstructuredGrid],
+ reference: pv.UnstructuredGrid,
mesh_property: propertylib.Property,
timestep_sizes: list[float],
) -> pd.DataFrame:
@@ -142,13 +142,13 @@ def convergence_metrics(
Calculate convergence metrics for a given reference and property.
:param meshes: The List of meshes to be analyzed for convergence.
- :param reference: The reference Dataset to compare against.
+ :param reference: The reference mesh to compare against.
:param mesh_property: The property of interest.
:returns: A pandas Dataframe containing all metrics.
"""
- def _data(m: pv.DataSet):
+ def _data(m: pv.UnstructuredGrid) -> np.ndarray:
return mesh_property.magnitude.transform(
m.point_data[mesh_property.data_name]
)
@@ -304,7 +304,9 @@ def convergence_metrics_evolution(
.to(units[1])
.magnitude
)
- p_metrics_per_t = np.concatenate(([time_vals], p_metrics_per_t.T)).T
+ p_metrics_per_t = np.concatenate(
+ (np.asarray([time_vals]), p_metrics_per_t.T)
+ ).T
columns = ["timevalue"] + [
f"{t} ({x})"
for t in ["abs. error", "rel. error", "p"]
diff --git a/ogstools/studies/convergence/examples/steady_state_diffusion.py b/ogstools/studies/convergence/examples/steady_state_diffusion.py
index 1144b180a7936481dba93f68e63c5311fe3ee776..aba2173017032ae3d53e5e36a4620c84106ea8b3 100644
--- a/ogstools/studies/convergence/examples/steady_state_diffusion.py
+++ b/ogstools/studies/convergence/examples/steady_state_diffusion.py
@@ -6,15 +6,15 @@ import numpy as np
import pyvista as pv
-def _c_k(k):
+def _c_k(k: float) -> float:
return 0.5 * (2 * k - 1) * np.pi
-def _a_k(k):
+def _a_k(k: float) -> float:
return 2 / (_c_k(k) ** 2 * np.cosh(_c_k(k)))
-def _h(points):
+def _h(points: np.ndarray) -> np.ndarray:
result = np.ones(len(points))
for k in np.arange(1, 100):
c_k_val = _c_k(k)
@@ -24,8 +24,14 @@ def _h(points):
return result
-def analytical_solution(topology: Union[Path, pv.DataSet]) -> pv.DataSet:
- mesh = topology if isinstance(topology, pv.DataSet) else pv.read(topology)
+def analytical_solution(
+ topology: Union[Path, pv.UnstructuredGrid]
+) -> pv.UnstructuredGrid:
+ mesh = (
+ topology
+ if isinstance(topology, pv.UnstructuredGrid)
+ else pv.read(topology)
+ )
new_mesh = deepcopy(mesh)
new_mesh.clear_point_data()
points = new_mesh.points
diff --git a/ogstools/studies/templates/convergence_study.py b/ogstools/studies/templates/convergence_study.py
index 4cee10d9e7dd16e2447ae19da203961609e78b86..709067f56c5a527a8efceb0aa68c4cbca9ad43fc 100644
--- a/ogstools/studies/templates/convergence_study.py
+++ b/ogstools/studies/templates/convergence_study.py
@@ -40,7 +40,7 @@ meshplotlib.setup.combined_colorbar = False
mesh_series = [meshlib.MeshSeries(mesh_path) for mesh_path in mesh_paths]
timestep_sizes = [np.mean(np.diff(ms.timevalues)) for ms in mesh_series]
meshes = [ms.read_closest(timevalue) for ms in mesh_series]
-topology: pv.DataSet = meshes[-3]
+topology: pv.UnstructuredGrid = meshes[-3]
mesh_property = propertylib.presets.get_preset(property_name, meshes[0])
richardson = studies.convergence.richardson_extrapolation(
meshes, mesh_property, topology, refinement_ratio
@@ -77,17 +77,17 @@ fig = meshplotlib.plot(richardson, mesh_property)
data_key = mesh_property.data_name
if reference_solution_path is None:
diff_mesh = meshlib.difference(
- mesh_property, richardson, topology.sample(meshes[-1])
+ richardson, topology.sample(meshes[-1]), mesh_property
)
- fig = meshplotlib.plot(diff_mesh, mesh_property.delta)
+ fig = meshplotlib.plot(diff_mesh, mesh_property)
else:
reference_solution = topology.sample(
meshlib.MeshSeries(reference_solution_path).read_closest(timevalue)
)
diff_mesh = meshlib.difference(
- mesh_property, reference_solution, richardson
+ reference_solution, richardson, mesh_property
)
- fig = meshplotlib.plot(diff_mesh, mesh_property.delta)
+ fig = meshplotlib.plot(diff_mesh, mesh_property)
# %% [markdown]
# ## Convergence metrics
diff --git a/pyproject.toml b/pyproject.toml
index 3ff6c2b35b5de17b25187130b44fa71cb0efe18e..f1d5d22d80da6eaed0da382afff7be66300847ff 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -26,6 +26,7 @@ dependencies = [
"ogs>=6.5.1.dev0",
"ogs6py>=0.370",
"tqdm>=4.60",
+ "typeguard>=4.0.0",
"VTUinterface>=0.704",
]
diff --git a/tests/test_meshlib.py b/tests/test_meshlib.py
index 2f46c458aaf1d256f460fb34e0556a3fd9c16136..89aa577ef1ebfc5f6ca95f0b3e6c397282fd6793 100644
--- a/tests/test_meshlib.py
+++ b/tests/test_meshlib.py
@@ -39,6 +39,26 @@ class UtilsTest(unittest.TestCase):
)
mesh_series.clear()
+ def test_aggregate(self):
+ "Test aggregation of meshseries."
+ mesh_series = MeshSeries(examples.xdmf_file)
+ funcs = ["min", "max", "mean", "median", "sum", "std", "var"]
+ for func in funcs:
+ agg_mesh = mesh_series.aggregate("temperature", func)
+ self.assertTrue(
+ not np.any(np.isnan(agg_mesh["temperature_" + func]))
+ )
+
+ def test_aggregate_mesh_dependent(self):
+ "Test aggregation of mesh_dependent property on meshseries."
+ mesh_series = MeshSeries(examples.pvd_file)
+ agg_mesh = mesh_series.aggregate(presets.dilatancy_alkan, "max")
+ self.assertTrue(
+ not np.any(
+ np.isnan(agg_mesh[presets.dilatancy_alkan.output_name + "_max"])
+ )
+ )
+
def test_probe_pvd(self):
"Test point probing on pvd."
mesh_series = MeshSeries(examples.pvd_file)
@@ -57,63 +77,43 @@ class UtilsTest(unittest.TestCase):
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)
+ mesh_diff = difference(mesh1, mesh2, "temperature")
+ mesh_diff = difference(mesh1, mesh2, presets.temperature)
self.assertTrue(isinstance(mesh_diff, UnstructuredGrid))
+ mesh_diff = difference(mesh1, mesh2)
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)
+ meshes_diff = difference_pairwise(meshes1, meshes2, presets.temperature)
self.assertTrue(
isinstance(meshes_diff, np.ndarray) and len(meshes_diff) == n
)
+ meshes_diff = difference_pairwise(meshes1, meshes2)
- def test_diff_matrix_single_list(self):
+ def test_diff_matrix_single(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))
+ meshes_diff = difference_matrix(
+ meshes1, mesh_property=presets.temperature
)
-
- 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))
)
+ meshes_diff = difference_matrix(meshes1)
- def test_diff_matrix_unequal_list(self):
+ def test_diff_matrix_unequal(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)
+ meshes_diff = difference_matrix(meshes1, meshes2, presets.temperature)
self.assertTrue(
isinstance(meshes_diff, np.ndarray)
and meshes_diff.shape == (len(meshes1), len(meshes2))
)
+ meshes_diff = difference_matrix(meshes1, meshes2)
diff --git a/tests/test_meshplotlib.py b/tests/test_meshplotlib.py
index 98f7271bf8ac1e791df0f29b5717b3de2e4f602a..e0f187451e0eb1d2deba3d91ca2e40720737cd84 100644
--- a/tests/test_meshplotlib.py
+++ b/tests/test_meshplotlib.py
@@ -14,7 +14,6 @@ from ogstools.meshplotlib import (
examples,
label_spatial_axes,
plot,
- plot_limit,
plot_probe,
setup,
update_font_sizes,
@@ -32,7 +31,9 @@ assert_allclose = partial(
class MeshplotlibTest(unittest.TestCase):
- """Test case for meshplotlib."""
+ """Test case for meshplotlib.
+
+ Most of these tests only test for no-throw, currently."""
def test_pyvista_offscreen(self):
import pyvista as pv
@@ -131,11 +132,22 @@ class MeshplotlibTest(unittest.TestCase):
plot_on_top(
fig.axes[0], mesh, lambda x: min(max(0, 0.1 * (x - 3)), 100)
)
+ plt.close()
def test_diff_plots(self):
"""Test creation of difference plots."""
- meshseries = examples.meshseries_CT_2D
- plot(difference("Si", meshseries.read(0), meshseries.read(1)), "Si")
+ meshseries = examples.meshseries_THM_2D
+ mesh0 = meshseries.read(0)
+ mesh1 = meshseries.read(1)
+ plot(difference(mesh1, mesh0, "temperature"), "temperature_difference")
+ for prop in [
+ presets.temperature,
+ presets.displacement,
+ presets.stress,
+ presets.stress.von_Mises,
+ ]:
+ plot(difference(mesh1, mesh0, prop), prop)
+ plt.close()
def test_user_defined_ax(self):
"""Test creating plot with subfigures and user provided ax"""
@@ -146,12 +158,13 @@ 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(
- presets.temperature, meshseries.read(0), meshseries.read(1)
+ meshseries.read(0), meshseries.read(1), presets.temperature
)
- plot(diff_mesh, presets.temperature.delta, fig=fig, ax=ax[2])
+ plot(diff_mesh, presets.temperature, fig=fig, ax=ax[2])
ax[2].set_title(r"$T(\mathrm{t}_{end})$-$T(\mathrm{t}_{0})$")
# fig.suptitle("Test user defined ax")
fig.tight_layout()
+ plt.close()
def test_user_defined_ax_diff_vals(self):
"""Test creating plot with subfigures and user provided ax with different values plotted"""
@@ -161,7 +174,7 @@ class MeshplotlibTest(unittest.TestCase):
plot(meshseries.read(0), presets.temperature, fig=fig, ax=ax[0])
plot(meshseries.read(1), presets.displacement, fig=fig, ax=ax[1])
fig.suptitle("Test user defined ax")
- fig.tight_layout()
+ plt.close()
def test_user_defined_fig(self):
"""Test creating plot with subfigures and user provided fig"""
@@ -174,6 +187,7 @@ class MeshplotlibTest(unittest.TestCase):
fig=fig,
)
fig.suptitle("Test user defined fig")
+ plt.close()
def test_update_font_sizes(self):
"""Test creating plot with subfigures and user provided fig"""
@@ -187,6 +201,7 @@ class MeshplotlibTest(unittest.TestCase):
)
fig = update_font_sizes(fig, fontsize=25)
fig.suptitle("Test user defined fig")
+ plt.close()
def test_sharexy(self):
"""Test if labels are skipped if axis are shared"""
@@ -197,11 +212,11 @@ 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(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()
+ diff_ab = difference(mesh_a, mesh_b, presets.temperature)
+ diff_ba = difference(mesh_b, mesh_a, presets.temperature)
+ plot(diff_ab, presets.temperature, fig=fig, ax=ax[2])
+ plot(diff_ba, presets.temperature, fig=fig, ax=ax[3])
+ plt.close()
def test_label_sharedxy(self):
"""Test labeling shared x and y axes"""
@@ -211,29 +226,31 @@ 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(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]))
- fig.tight_layout()
+ diff_ab = difference(mesh_a, mesh_b, presets.temperature)
+ diff_ba = difference(mesh_b, mesh_a, presets.temperature)
+ plot(diff_ab, presets.temperature, fig=fig, ax=ax[0][1])
+ plot(diff_ba, presets.temperature, fig=fig, ax=ax[1][1])
+ label_spatial_axes(ax, "x", "y")
+ plt.close()
def test_spatial_label(self):
- """Test if labels are added to x and y axes"""
+ """Test axes labeling"""
fig, ax = plt.subplots(2, 2)
- ax = label_spatial_axes(ax, np.array([0, 1]))
+ label_spatial_axes(ax, "x", "y")
+ plt.close()
def test_spatial_label_clear(self):
- """Test if labels are added to x and y axes"""
+ """Test axes labels clearing"""
fig, ax = plt.subplots(2, 2)
- ax = label_spatial_axes(ax, np.array([0, 1]))
- ax = clear_labels(ax)
+ label_spatial_axes(ax, "x", "y")
+ clear_labels(ax)
+ plt.close()
def test_limit_plots(self):
"""Test creation of limit plots."""
- meshseries = examples.meshseries_CT_2D
- plot_limit(meshseries, "Si", "min")
- plot_limit(meshseries, "Si", "max")
+ mesh = examples.meshseries_CT_2D.aggregate("Si", "var")
+ plot(mesh, "Si_var")
+ plt.close()
def test_plot_probe(self):
"""Test creation of probe plots."""
@@ -250,6 +267,7 @@ class MeshplotlibTest(unittest.TestCase):
plot_probe(mesh_series, points, presets.temperature)
mesh_property = presets.velocity.replace(data_name="darcy_velocity")
plot_probe(mesh_series, points, mesh_property)
+ plt.close()
def test_animation(self):
"""Test creation of animation."""
@@ -258,6 +276,7 @@ class MeshplotlibTest(unittest.TestCase):
titles = [str(tv) for tv in timevalues]
anim = animate(meshseries, presets.temperature, timevalues, titles)
anim.to_jshtml()
+ plt.close()
def test_save_animation(self):
"""Test saving of an animation."""
@@ -266,6 +285,7 @@ class MeshplotlibTest(unittest.TestCase):
anim = animate(meshseries, presets.temperature, timevalues)
if not save_animation(anim, mkstemp()[1], 5):
self.skipTest("Saving animation failed.")
+ plt.close()
def test_plot_3D(self):
"""Test creation of slice plots for 3D mesh."""
@@ -274,13 +294,16 @@ class MeshplotlibTest(unittest.TestCase):
meshes = np.reshape(mesh.slice_along_axis(4, "x"), (2, 2))
plot(meshes, "Spatial Point Data")
plot(mesh.slice([1, -2, 0]), "Spatial Point Data")
+ plt.close()
def test_xdmf(self):
"""Test creation of 2D plots from xdmf data."""
mesh = examples.meshseries_CT_2D.read(0)
plot(mesh, Scalar("Si"))
+ plt.close()
def test_xdmf_with_slices(self):
"""Test creation of 2D plots from xdmf data."""
mesh = examples.meshseries_XDMF.read(0)
plot(mesh, presets.pressure)
+ plt.close()