diff --git a/docs/examples/howto_meshlib/plot_meshlib_pyvista_input.py b/docs/examples/howto_meshlib/plot_meshlib_pyvista_input.py
index a457de915b1bf4244999a6482a9b4c5f8e711981..6ad01daf3e7148639346ad69b446c91cd8f29382 100644
--- a/docs/examples/howto_meshlib/plot_meshlib_pyvista_input.py
+++ b/docs/examples/howto_meshlib/plot_meshlib_pyvista_input.py
@@ -8,15 +8,10 @@ For this example we create meshes from pyvista surfaces.
"""
# %%
-import numpy as np # For visualization only
-
-import ogstools.meshplotlib as mpl # For visualization only
from ogstools.meshlib.boundary import Layer
from ogstools.meshlib.boundary_set import LayerSet
from ogstools.meshlib.boundary_subset import Gaussian2D, Surface
-from ogstools.meshlib.region import (
- to_region_tetraeder,
-)
+from ogstools.meshlib.region import to_region_tetraeder
# See other examples for different meshing algorithms
@@ -30,13 +25,9 @@ surface3 = Surface(Gaussian2D(**args, height_offset=-200), material_id=2)
ls = LayerSet([Layer(surface1, surface2), Layer(surface2, surface3)])
mesh = to_region_tetraeder(ls, 40).mesh
+# interprets values as categories
+mesh["regions"] = [str(m) for m in mesh["MaterialIDs"]]
# %%
# Visualize the prism mesh
-
-slices = np.reshape(mesh.slice_along_axis(n=4, axis="y"), (1, -1))
-mpl.setup.aspect_limits = [0.2, 5.0]
-fig = mpl.plot(slices, "MaterialIDs")
-for ax, slice in zip(fig.axes, np.ravel(slices)):
- ax.set_title(f"z = {slice.center[2]:.1f} {mpl.setup.length.data_unit}")
-_ = fig
+mesh.plot(scalars="regions", show_edges=True)
diff --git a/docs/examples/howto_meshlib/plot_meshlib_vtu_input.py b/docs/examples/howto_meshlib/plot_meshlib_vtu_input.py
index a2956693ca917b3d2cbeee87f221bcf561c9bb91..f8b60287e28c1fd6d04f83e8e3c89b5aafcab318 100644
--- a/docs/examples/howto_meshlib/plot_meshlib_vtu_input.py
+++ b/docs/examples/howto_meshlib/plot_meshlib_vtu_input.py
@@ -10,9 +10,6 @@ For this example we create meshes from surface layers.
# %%
from pathlib import Path # To get example vtu files
-import numpy as np # For visualization only
-
-import ogstools.meshplotlib as mpl # For visualization only
from ogstools.definitions import TESTS_DIR # To get example vtu files
from ogstools.meshlib.boundary import Layer
from ogstools.meshlib.boundary_set import LayerSet
@@ -24,10 +21,6 @@ from ogstools.meshlib.region import (
to_region_voxel,
)
-mpl.setup.reset()
-mpl.setup.length.output_unit = "km"
-mpl.setup.aspect_limits = [0.2, 5.0]
-
# %%
# The loaded surfaces are defined within VTU files and adhere to properties such as non-intersecting boundaries with consistent x and y bounds. Alternatively, surfaces can also be created using PyVista with the same properties.
@@ -54,28 +47,30 @@ pm = to_region_prism(layer_set1, resolution=200).mesh
vm = to_region_voxel(layer_set1, resolution=[200, 200, 50]).mesh
tm = to_region_tetraeder(layer_set1, resolution=200).mesh
-
+# %% [markdown]
+# Simplified mesh
+# ---------------
# %%
-# Visualize the prism mesh
-
-mesh = pm
-slices = np.reshape(mesh.slice_along_axis(n=4, axis="y"), (-1, 1))
-fig = mpl.plot(slices, "MaterialIDs")
-for ax, slice in zip(fig.axes, np.ravel(slices)):
- ax.set_title(f"z = {slice.center[2]:.1f} {mpl.setup.length.data_unit}")
+sm["regions"] = [str(m) for m in sm["MaterialIDs"]]
+sm.scale([1, 1, 5]).plot(scalars="regions", show_edges=True)
+# %% [markdown]
+# Voxel mesh
+# ---------------
# %%
-# Visualize meshes with different meshing algorithm
+vm["regions"] = [str(m) for m in vm["MaterialIDs"]]
+vm.scale([1, 1, 5]).plot(scalars="regions", show_edges=True)
-meshes = [sm, vm, pm, tm]
-names = [
- "to_region_simplified",
- "to_region_voxel",
- "to_region_prism",
- "to_region_tetraeder",
-]
+# %% [markdown]
+# Prism mesh
+# ---------------
+# %%
+pm["regions"] = [str(m) for m in pm["MaterialIDs"]]
+pm.scale([1, 1, 5]).plot(scalars="regions", show_edges=True)
-x_slices = np.reshape([mesh.slice("x") for mesh in meshes], (-1, 1))
-fig = mpl.plot(x_slices, "MaterialIDs")
-for ax, name in zip(fig.axes, names):
- ax.set_title(name)
+# %% [markdown]
+# Tetraeder mesh
+# ---------------
+# %%
+tm["regions"] = [str(m) for m in tm["MaterialIDs"]]
+tm.scale([1, 1, 5]).plot(scalars="regions", show_edges=True)
diff --git a/docs/examples/howto_meshplotlib/plot_animation.py b/docs/examples/howto_meshplotlib/plot_animation.py
index 8e5ab568dabd5612c9de9da98bc6e5580a8542d0..47ff7a8612abdb2ab61966b8e69fb3079d248649 100644
--- a/docs/examples/howto_meshplotlib/plot_animation.py
+++ b/docs/examples/howto_meshplotlib/plot_animation.py
@@ -41,7 +41,9 @@ timevalues = np.linspace(
# 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("Si", "", "%", "Saturation")
+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
diff --git a/docs/examples/howto_meshplotlib/plot_aspect_ratios.py b/docs/examples/howto_meshplotlib/plot_aspect_ratios.py
new file mode 100644
index 0000000000000000000000000000000000000000..015c7cfe4ba08d0f67b257407db23d4643c75547
--- /dev/null
+++ b/docs/examples/howto_meshplotlib/plot_aspect_ratios.py
@@ -0,0 +1,90 @@
+"""
+Aspect ratios
+=============
+
+.. sectionauthor:: Florian Zill (Helmholtz Centre for Environmental Research GmbH - UFZ)
+
+By default plots with meshplotlib try to retain the true mesh proportions.
+If the meshes aspect ratio lies outside of predefined limits
+(setup.min_ax_aspect, setup.max_ax_aspect) the axes get compressed to stay
+inside the given limits. The following examples shall illustrate this
+behaviour.
+"""
+
+# %%
+import numpy as np
+import pyvista as pv
+
+from ogstools.meshplotlib import plot, setup
+
+setup.reset()
+print(f"{setup.min_ax_aspect=}")
+print(f"{setup.max_ax_aspect=}")
+
+
+# sphinx_gallery_start_ignore
+def custom_mesh(dx: float, dy: float):
+ number_of_points = 50
+ x = np.linspace(0, dx, num=number_of_points)
+ y = np.linspace(0, dy, num=number_of_points)
+ xx, yy = np.meshgrid(x, y)
+ zz = np.zeros(xx.shape)
+ points = np.c_[xx.reshape(-1), yy.reshape(-1), zz.reshape(-1)]
+ mesh = pv.PolyData(points).delaunay_2d()
+ mesh.point_data["example"] = np.sin(mesh.points[:, 0]) * np.cos(
+ mesh.points[:, 1]
+ )
+ return mesh
+
+
+# sphinx_gallery_end_ignore
+
+# %% [markdown]
+# The following fits inside the defined limits and gets displayed with true
+# proportions.
+
+# %%
+fig = plot(custom_mesh(np.pi * 2, np.pi), "example")
+# %% [markdown]
+# This one would be too wide and thus and gets compressed to fit the maximum
+# aspect ratio.
+
+# %%
+fig = plot(custom_mesh(np.pi * 4, np.pi), "example")
+# %% [markdown]
+# When plotting multiple meshes together, this applies to each subplot.
+# So here each subplot has true proportions again since each one fits the limits.
+
+# %%
+fig = plot(
+ [custom_mesh(np.pi * 2, np.pi), custom_mesh(np.pi * 2, np.pi)], "example"
+)
+# %% [markdown]
+# The following figure would be to tall and is clipped to the minimum aspect
+# ratio.
+
+# %%
+fig = plot(custom_mesh(np.pi, np.pi * 2), "example")
+# %% [markdown]
+# The same is true here:
+
+# %%
+fig = plot(
+ [custom_mesh(np.pi, np.pi * 2), custom_mesh(np.pi, np.pi * 2)], "example"
+)
+
+# %% [markdown]
+# You can enforce true proportions regardless of the resulting figures
+# dimensions, by setting the limiting values to None. In this case we get a
+# very wide figure.
+
+# %%
+setup.min_ax_aspect = None
+setup.max_ax_aspect = None
+fig = plot(custom_mesh(np.pi * 3, np.pi), "example")
+
+# %% [markdown]
+# And in this case we get a very tall figure.
+
+# %%
+fig = plot(custom_mesh(np.pi, np.pi * 3), "example")
diff --git a/docs/examples/howto_meshplotlib/plot_meshplotlib_2d.py b/docs/examples/howto_meshplotlib/plot_meshplotlib_2d.py
index d1722ffda1fa9662016533c686fb1354d4a3281b..9be7ab03bbee74806a3a276ccb9f91e8a9ad28f3 100644
--- a/docs/examples/howto_meshplotlib/plot_meshplotlib_2d.py
+++ b/docs/examples/howto_meshplotlib/plot_meshplotlib_2d.py
@@ -37,11 +37,6 @@ fig = mpl.plot(mesh, presets.pressure.get_mask())
# Let's plot the fluid velocity field.
fig = mpl.plot(mesh, presets.velocity)
-# %%
-# It is also possible to plot a shape on top, e.g. to display an overburden.
-mpl.plot_on_top(fig.axes[0], mesh, lambda x: min(max(0, 0.1 * (x - 3)), 100))
-fig # noqa: B018
-
# %%
# We can also plot components of vector property:
diff --git a/docs/examples/howto_meshplotlib/plot_special.py b/docs/examples/howto_meshplotlib/plot_special.py
new file mode 100644
index 0000000000000000000000000000000000000000..77b0283435dbdc4785cb8a2364165326c50ce604
--- /dev/null
+++ b/docs/examples/howto_meshplotlib/plot_special.py
@@ -0,0 +1,46 @@
+"""
+How to plot limits and differences
+==================================
+
+.. 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
+a component transport example from the ogs benchmark gallery
+(https://www.opengeosys.org/docs/benchmarks/hydro-component/elder/).
+
+To see this benchmark results over all timesteps have a look at
+:ref:`sphx_glr_auto_examples_howto_meshplotlib_plot_animation.py`.
+"""
+
+# %%
+from ogstools.meshplotlib import examples, plot_diff, plot_limit, setup
+from ogstools.propertylib import Scalar
+
+setup.reset()
+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
+
+# %%
+fig = plot_limit(mesh_series, si, "max")
+
+# %% [markdown]
+# Minimum though all timesteps
+
+# %%
+fig = plot_limit(mesh_series, si, "min")
+
+
+# %% [markdown]
+# Difference between the last and he first timestep:
+
+# %%
+fig = plot_diff(mesh_series.read(-1), mesh_series.read(0), si)
diff --git a/ogstools/meshlib/mesh_series.py b/ogstools/meshlib/mesh_series.py
index 71ebe5e3f04d320603f84556ae969159b4f987cb..552e13f34fc3c8604bde450bf700b2e9d7f4b37e 100644
--- a/ogstools/meshlib/mesh_series.py
+++ b/ogstools/meshlib/mesh_series.py
@@ -11,6 +11,7 @@ from meshio.xdmf.time_series import (
cell_data_from_raw,
xdmf_to_numpy_type,
)
+from tqdm.auto import tqdm
class TimeSeriesReader(meshio.xdmf.TimeSeriesReader):
@@ -156,12 +157,22 @@ class MeshSeries:
self._pvd_reader = pv.PVDReader(filepath)
elif self._data_type == "xdmf":
self._xdmf_reader = TimeSeriesReader(filepath)
+ self._read_xdmf(0) # necessary to initialize hdf5_files
+ meshes = self.hdf5["meshes"]
+ self.hdf5_bulk_name = list(meshes.keys())[
+ np.argmax([meshes[m]["geometry"].shape[1] for m in meshes])
+ ]
elif self._data_type == "vtu":
self._vtu_reader = pv.XMLUnstructuredGridReader(filepath)
else:
msg = "Can only read 'pvd', 'xdmf' or 'vtu' files."
raise TypeError(msg)
+ @property
+ def hdf5(self):
+ # We assume there is only one h5 file
+ return next(iter(self._xdmf_reader.hdf5_files.values()))
+
def _read_pvd(self, timestep: int) -> pv.UnstructuredGrid:
self._pvd_reader.set_active_time_point(timestep)
return self._pvd_reader.read()[0]
@@ -251,3 +262,18 @@ class MeshSeries:
timevalue - t_vals[ts1]
)
return mesh
+
+ def values(self, data_name: str) -> np.ndarray:
+ """
+ Get the data in the MeshSeries for all timesteps.
+
+ :param data_name: Name of the data in the MeshSeries.
+
+ :returns: A numpy array of the requested data for all timesteps
+ """
+ mesh = self.read(0).copy()
+ if self._data_type == "xdmf":
+ return self.hdf5["meshes"][self.hdf5_bulk_name][data_name]
+ if self._data_type == "pvd":
+ return [self.read(t)[data_name] for t in tqdm(self.timesteps)]
+ return mesh[data_name]
diff --git a/ogstools/meshplotlib/__init__.py b/ogstools/meshplotlib/__init__.py
index acd2a332c58e484e05720eeb9f4bd0a6cacf6570..3702477db2f628fab93fe6f853e1748e22bb6355 100644
--- a/ogstools/meshplotlib/__init__.py
+++ b/ogstools/meshplotlib/__init__.py
@@ -3,13 +3,15 @@
from .plot_setup import _setup as setup # noqa: I001: noqa
-from .core import plot, subplot
+from .core import plot_diff, plot_limit, plot, subplot
from .plot_features import plot_contour, plot_on_top
__all__ = [
"setup",
"plot",
- "plot_on_top",
"plot_contour",
+ "plot_diff",
+ "plot_limit",
+ "plot_on_top",
"subplot",
]
diff --git a/ogstools/meshplotlib/animation.py b/ogstools/meshplotlib/animation.py
index 0df18f1cc812e25a312bc0e0e9fcbda54af6d99a..295d97e42b756873e41842e65a034a6c94406e30 100644
--- a/ogstools/meshplotlib/animation.py
+++ b/ogstools/meshplotlib/animation.py
@@ -1,8 +1,7 @@
import warnings
from collections.abc import Sequence
from functools import partial
-from typing import Optional as Opt
-from typing import Union
+from typing import Optional, Union
import matplotlib as mpl
import numpy as np
@@ -20,8 +19,8 @@ from .core import _plot_on_figure, plot
def animate(
mesh_series: MeshSeries,
property: Property,
- timesteps: Opt[Sequence] = None,
- titles: Opt[list[str]] = None,
+ timesteps: Optional[Sequence] = None,
+ titles: Optional[list[str]] = None,
) -> FuncAnimation:
"""
Create an animation for a property of a mesh series with given timesteps.
diff --git a/ogstools/meshplotlib/core.py b/ogstools/meshplotlib/core.py
index 0968c7eba11117d70192e3150d2a78de6ffad3a3..54f7755a1c9567b70dba0f98ca3c9382111db765 100644
--- a/ogstools/meshplotlib/core.py
+++ b/ogstools/meshplotlib/core.py
@@ -1,7 +1,8 @@
"""Meshplotlib core utilitites."""
+
import types
-from typing import Optional as Opt
-from typing import Union
+from copy import deepcopy
+from typing import Literal, Optional, Union
import numpy as np
import pyvista as pv
@@ -14,6 +15,7 @@ from matplotlib import ticker as mticker
from matplotlib import transforms as mtransforms
from matplotlib.patches import Rectangle as Rect
+from ogstools.meshlib import MeshSeries
from ogstools.propertylib import Property, Vector
from ogstools.propertylib.presets import _resolve_property
@@ -30,14 +32,29 @@ def _q_zero_line(property: Property, levels: np.ndarray):
)
-def get_data(mesh: pv.UnstructuredGrid, property: Property) -> np.ndarray:
+def has_masked_values(mesh: pv.UnstructuredGrid, property: Property) -> bool:
+ return (
+ not property.is_mask()
+ and property.mask in mesh.cell_data
+ and (len(mesh.cell_data[property.mask]) != 0)
+ )
+
+
+def get_data(
+ mesh: pv.UnstructuredGrid, property: Property, masked: bool = False
+) -> np.ndarray:
"""Get the data associated with a scalar or vector property from a mesh."""
- if property.data_name in mesh.point_data:
- return mesh.point_data[property.data_name]
- if property.data_name in mesh.cell_data:
- return mesh.cell_data[property.data_name]
- msg = f"Property not found in mesh {mesh}."
- raise IndexError(msg)
+ if (
+ property.data_name not in mesh.point_data
+ and property.data_name not in mesh.cell_data
+ ):
+ msg = f"Property not found in mesh {mesh}."
+ raise IndexError(msg)
+ if masked:
+ return mesh.ctp(True).threshold(value=[1, 1], scalars=property.mask)[
+ property.data_name
+ ]
+ return mesh[property.data_name]
def get_level_boundaries(levels: np.ndarray):
@@ -69,7 +86,6 @@ def get_cmap_norm(
conti_cmap = mcolors.LinearSegmentedColormap.from_list(
"temperature_cmap", np.vstack((cool_colors, warm_colors))
)
- bilinear = vmin < 0 < vmax
if bilinear:
vmin, vmax = np.max(np.abs([vmin, vmax])) * np.array([-1.0, 1.0])
if property.categoric:
@@ -113,6 +129,7 @@ def add_colorbars(
property: Property,
levels: np.ndarray,
pad: float = 0.05,
+ labelsize: Optional[float] = None,
) -> None:
"""Add a colorbar to the matplotlib figure."""
cmap, norm = get_cmap_norm(levels, property)
@@ -137,13 +154,13 @@ def add_colorbars(
unit_str = (
f" / {property.get_output_unit()}" if property.get_output_unit() else ""
)
- cb.set_label(
- property.output_name.replace("_", " ") + unit_str,
- size=setup.rcParams_scaled["font.size"],
+ labelsize = (
+ setup.rcParams_scaled["font.size"] if labelsize is None else labelsize
)
- cb.ax.tick_params(
- labelsize=setup.rcParams_scaled["font.size"], direction="out"
+ cb.set_label(
+ property.output_name.replace("_", " ") + unit_str, size=labelsize
)
+ cb.ax.tick_params(labelsize=labelsize, direction="out")
mf = mticker.ScalarFormatter(useMathText=True, useOffset=True)
mf.set_scientific(True)
mf.set_powerlimits([-3, 3])
@@ -175,7 +192,7 @@ def subplot(
mesh: pv.UnstructuredGrid,
property: Union[Property, str],
ax: plt.Axes,
- levels: Opt[np.ndarray] = None,
+ levels: Optional[np.ndarray] = None,
) -> None:
"""
Plot the property field of a mesh on a matplotlib.axis.
@@ -194,11 +211,7 @@ def subplot(
ax.axis("auto")
- if (
- not property.is_mask()
- and property.mask in mesh.cell_data
- and len(mesh.cell_data[property.mask])
- ):
+ if has_masked_values(mesh, property):
subplot(mesh, property.get_mask(), ax)
mesh = mesh.ctp(True).threshold(value=[1, 1], scalars=property.mask)
@@ -239,15 +252,15 @@ def subplot(
surf = mesh.extract_surface()
- if setup.show_region_bounds and "MaterialIDs" in mesh.cell_data:
- pf.plot_layer_boundaries(ax, surf, projection)
-
show_edges = setup.show_element_edges
if isinstance(setup.show_element_edges, str):
show_edges = setup.show_element_edges == property.data_name
if show_edges:
pf.plot_element_edges(ax, surf, projection)
+ if setup.show_region_bounds and "MaterialIDs" in mesh.cell_data:
+ pf.plot_layer_boundaries(ax, surf, projection)
+
if isinstance(property, Vector):
pf.plot_streamlines(ax, surf_tri, property, projection)
@@ -302,7 +315,7 @@ def _get_rows_cols(
# TODO: fixed_figure_size -> ax aspect automatic
-def _fig_init(rows: int, cols: int, ax_aspect: float = 1.0) -> mfigure.Figure:
+def _fig_init(rows: int, cols: int, aspect: float = 1.0) -> mfigure.Figure:
nx_cb = 1 if setup.combined_colorbar else cols
default_size = 8
cb_width = 4
@@ -310,14 +323,10 @@ def _fig_init(rows: int, cols: int, ax_aspect: float = 1.0) -> mfigure.Figure:
x_label_height = 1
figsize = setup.fig_scale * np.asarray(
[
- default_size * cols * ax_aspect + cb_width * nx_cb + y_label_width,
+ default_size * cols * aspect + cb_width * nx_cb + y_label_width,
default_size * rows + x_label_height,
]
)
- if figsize[0] / figsize[1] > setup.fig_aspect_limits[1]:
- figsize[0] = figsize[1] * setup.fig_aspect_limits[1]
- elif figsize[0] / figsize[1] < setup.fig_aspect_limits[0]:
- figsize[0] = figsize[1] * setup.fig_aspect_limits[0]
fig, _ = plt.subplots(
rows,
cols,
@@ -343,7 +352,9 @@ def get_combined_levels(
p_min, p_max = np.inf, -np.inf
unique_vals = np.array([])
for mesh in np.ravel(meshes):
- values = property.magnitude.strip_units(get_data(mesh, property))
+ values = property.magnitude.strip_units(
+ get_data(mesh, property, masked=has_masked_values(mesh, property))
+ )
if setup.log_scaled: # TODO: can be improved
values = np.log10(np.where(values > 1e-14, values, 1e-14))
p_min = min(p_min, np.nanmin(values)) if setup.p_min is None else p_min
@@ -445,23 +456,64 @@ def plot(
data_shape = _meshes[0][property].shape
property = _resolve_property(property, data_shape)
data_aspects = np.asarray([get_data_aspect(mesh) for mesh in _meshes])
- data_aspects[
- (data_aspects > setup.ax_aspect_limits[0])
- & (data_aspects < setup.ax_aspect_limits[1])
- ] = 1.0
- clamped_aspects = np.clip(data_aspects, *setup.ax_aspect_limits)
- ax_aspects = data_aspects / clamped_aspects
- _fig = _fig_init(
- rows=shape[0], cols=shape[1], ax_aspect=np.mean(ax_aspects)
- )
+ if setup.min_ax_aspect is None and setup.max_ax_aspect is None:
+ fig_aspect = np.mean(data_aspects)
+ else:
+ fig_aspect = np.mean(
+ np.clip(data_aspects, setup.min_ax_aspect, setup.max_ax_aspect)
+ )
+ ax_aspects = fig_aspect / data_aspects
+ _fig = _fig_init(rows=shape[0], cols=shape[1], aspect=fig_aspect)
n_axs = shape[0] * shape[1]
- # setting the aspect twice is intended
- # the first time results in properly spaced ticks
- # the second time fixes any deviations from the set aspect due to
- # additional colorbar(s) or secondary axes.
- for ax, aspect in zip(_fig.axes[: n_axs + 1], clamped_aspects):
- ax.set_aspect(aspect)
fig = _plot_on_figure(_fig, meshes, property)
- for ax, aspect in zip(fig.axes[: n_axs + 1], clamped_aspects):
- ax.set_aspect(aspect)
+ for ax, aspect in zip(fig.axes[: n_axs + 1], ax_aspects):
+ ax.set_aspect(1.0 / aspect)
return fig
+
+
+def plot_diff(
+ mesh1: pv.UnstructuredGrid,
+ mesh2: pv.UnstructuredGrid,
+ property: Union[Property, str],
+) -> mfigure.Figure:
+ if isinstance(property, str):
+ data_shape = mesh1[property].shape
+ property = _resolve_property(property, data_shape)
+ diff_mesh = deepcopy(mesh1)
+ diff_mesh[property.data_name] -= mesh2[property.data_name]
+ data_property = property.replace(output_unit=property.data_unit)
+ diff_unit = (data_property(1) - data_property(1)).units
+ diff_property = property.replace(
+ data_unit=diff_unit,
+ output_unit=diff_unit,
+ output_name=property.output_name + " difference",
+ bilinear_cmap=True,
+ )
+ return plot(diff_mesh, diff_property)
+
+
+def plot_limit(
+ mesh_series: MeshSeries,
+ property: Union[Property, str],
+ limit: Literal["min", "max"],
+) -> 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
+
+ :returns: A matplotlib Figure
+ """
+ mesh = mesh_series.read(0)
+ if isinstance(property, str):
+ data_shape = mesh[property].shape
+ property = _resolve_property(property, data_shape)
+ func = {"min": np.min, "max": np.max}[limit]
+ vals = mesh_series.values(property.data_name)
+ func(vals, out=mesh[property.data_name], axis=0)
+ limit_property = property.replace(
+ output_name=limit + " " + property.output_name
+ )
+ return plot(mesh, limit_property)
diff --git a/ogstools/meshplotlib/plot_setup.py b/ogstools/meshplotlib/plot_setup.py
index 9ca591098bf6014a5d6fb7c2cc3ec3346f97c9ab..c50115be9626a5c769a691fba4c80d1189432015 100644
--- a/ogstools/meshplotlib/plot_setup.py
+++ b/ogstools/meshplotlib/plot_setup.py
@@ -1,6 +1,5 @@
"""Plot configuration setup."""
-
from dataclasses import dataclass
from typing import Union
@@ -36,13 +35,10 @@ class PlotSetup:
"The resolution (dots per inch) for the figure."
fig_scale: float
"A scaling factor for the figure."
- ax_aspect_limits: tuple[float, float]
- "Lower and upper limit of the ax aspect ratio. For meshes with data ratios"
- "outside these bounds the aspect ratio gets clamped. Inside these bounds,"
- "the axes keep an aspect ratio of 1 (units are equally long on both axes)."
- fig_aspect_limits: tuple[float, float]
- "Lower and upper limit of the figure aspect ratio. If a figure would exceed"
- "them, they get clipped to prevent overly long or wide figures."
+ min_ax_aspect: float
+ "Minimum aspect ratio of subplots."
+ max_ax_aspect: float
+ "Maximum aspect ratio of subplots."
invert_colorbar: bool
"A boolean indicating whether to invert the colorbar."
layout: str
@@ -109,8 +105,8 @@ class PlotSetup:
"""Create a PlotSetup instance from a dictionary."""
return cls(
fig_scale=obj["fig_scale"],
- ax_aspect_limits=obj["ax_aspect_limits"],
- fig_aspect_limits=obj["fig_aspect_limits"],
+ min_ax_aspect=obj["min_ax_aspect"],
+ max_ax_aspect=obj["max_ax_aspect"],
invert_colorbar=obj["invert_colorbar"],
dpi=obj["dpi"],
num_levels=obj["num_levels"],
diff --git a/ogstools/meshplotlib/plot_setup_defaults.py b/ogstools/meshplotlib/plot_setup_defaults.py
index 42f14b713d6d02287d72fa71c7718996590865f2..29810fa350d783dfc11ac55f1a018040c8d492c1 100644
--- a/ogstools/meshplotlib/plot_setup_defaults.py
+++ b/ogstools/meshplotlib/plot_setup_defaults.py
@@ -10,10 +10,10 @@ setup_dict = {
"default_cmap": "RdBu_r",
"dpi": 120,
"fig_scale": 1.0,
- "ax_aspect_limits": (0.4, 2.2),
- "fig_aspect_limits": (0.2, 5.0),
+ "min_ax_aspect": 1.0,
+ "max_ax_aspect": 2.0,
"invert_colorbar": False,
- "layout": "constrained",
+ "layout": "compressed",
"length": ["m", "m"],
"material_names": None,
"num_levels": 11,
@@ -34,12 +34,15 @@ setup_dict = {
"custom_cmap": None,
"cmap_dict": {
"displacement": "Greens",
- "temperature": ["Blues", "plasma"],
+ "temperature": "plasma",
"pressure": "Blues",
"velocity": "coolwarm",
"MaterialIDs": "tab20",
},
- "cmap_dict_if_bilinear": {"displacement": "PRGn"},
+ "cmap_dict_if_bilinear": {
+ "displacement": "PRGn",
+ "temperature": ["Blues", "plasma"],
+ },
"cmap_if_mask": ["lightgrey", "green"],
"rcParams": {
"font.weight": "normal",
diff --git a/ogstools/physics/nuclearwasteheat/nuclearwaste.py b/ogstools/physics/nuclearwasteheat/nuclearwaste.py
index 3377556846152cb681efcc5303c82e8929718b6c..86f417d35ef41876fb6f258c3ba04f207dd9140c 100644
--- a/ogstools/physics/nuclearwasteheat/nuclearwaste.py
+++ b/ogstools/physics/nuclearwasteheat/nuclearwaste.py
@@ -1,6 +1,5 @@
from dataclasses import dataclass
-from typing import Optional as Opt
-from typing import Union
+from typing import Optional, Union
import numpy as np
from pint.facets.plain import PlainQuantity
@@ -31,7 +30,7 @@ class NuclearWaste:
self,
t: Union[PlainQuantity, float, np.ndarray],
baseline: bool = False,
- ncl_id: Opt[int] = None,
+ ncl_id: Optional[int] = None,
time_unit: str = "s",
power_unit: str = "W",
) -> Union[float, np.ndarray]:
@@ -84,7 +83,7 @@ class Repository:
self,
t: Union[PlainQuantity, float, np.ndarray],
baseline: bool = False,
- ncl_id: Opt[int] = None,
+ ncl_id: Optional[int] = None,
time_unit: str = "s",
power_unit: str = "W",
) -> Union[float, np.ndarray]:
diff --git a/ogstools/propertylib/presets.py b/ogstools/propertylib/presets.py
index 8712f01803d986453516f34c3626e025c0eeaa07..93bc44b6bdf9a0a6dc56e7fe0b968415237a83a3 100644
--- a/ogstools/propertylib/presets.py
+++ b/ogstools/propertylib/presets.py
@@ -22,7 +22,7 @@ hydraulic_height = Scalar("pressure", "m", "m", "hydraulic_height", H_mask)
qp_ratio = Scalar("sigma", "Pa", "percent", "QP_ratio", M_mask, v2s.qp_ratio)
strain = Matrix("epsilon", "", "percent", "strain", M_mask)
stress = Matrix("sigma", "Pa", "MPa", "stress", M_mask)
-temperature = Scalar("temperature", "K", "°C", mask=T_mask)
+temperature = Scalar("temperature", "K", "°C", mask=T_mask, bilinear_cmap=True)
velocity = Vector("velocity", "m/s", "m/s", "darcy_velocity", H_mask)
von_mises_stress = Scalar("sigma", "Pa", "MPa", "von_Mises_stress", M_mask, v2s.von_mises)
# fmt: on
diff --git a/pyproject.toml b/pyproject.toml
index c7d831d433f177bf4141a6c9b23a7354fc2255de..940c1ba39eca158119543c8053abb379308875f0 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -16,7 +16,7 @@ dependencies = [
"meshio==5.3.4",
"matplotlib>=3.7.1",
"nbconvert>=7.9.2",
- "pyvista[all,trame]>=0.39.1",
+ "pyvista[all,trame]>=0.40.1",
"rich>=13.4.2",
"scipy>=1.10.1",
"Pint>=0.22",
diff --git a/tests/test_meshplotlib.py b/tests/test_meshplotlib.py
index e51f85c7e5f2d8c6a998328552d4da4bb5449331..6113d32282e643101961abfcb29c03195e25fdb2 100644
--- a/tests/test_meshplotlib.py
+++ b/tests/test_meshplotlib.py
@@ -7,7 +7,7 @@ from tempfile import mkstemp
import numpy as np
from pyvista import examples as pv_examples
-from ogstools.meshplotlib import examples, plot, setup
+from ogstools.meshplotlib import examples, plot, plot_diff, plot_limit, setup
from ogstools.meshplotlib.animation import animate, save_animation
from ogstools.meshplotlib.levels import get_levels
from ogstools.meshplotlib.plot_features import plot_on_top
@@ -58,6 +58,17 @@ class MeshplotlibTest(unittest.TestCase):
fig.axes[0], mesh, lambda x: min(max(0, 0.1 * (x - 3)), 100)
)
+ def test_diff_plots(self):
+ """Test creation of difference plots."""
+ meshseries = examples.meshseries_CT_2D
+ plot_diff(meshseries.read(0), meshseries.read(1), "Si")
+
+ 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")
+
def test_animation(self):
"""Test creation of animation."""
meshseries = examples.meshseries_THM_2D
@@ -67,7 +78,7 @@ class MeshplotlibTest(unittest.TestCase):
anim.to_jshtml()
def test_save_animation(self):
- """Test creation of animation."""
+ """Test saving of an animation."""
meshseries = examples.meshseries_THM_2D
timevalues = np.linspace(0, meshseries.timevalues[-1], num=3)
anim = animate(meshseries, presets.temperature, timevalues)