diff --git a/ogstools/meshplotlib/core.py b/ogstools/meshplotlib/core.py
index 4e7284265b7f7036748c9075c173714bcab983d0..86c9475f658775bb814bed1932c17e2539d255aa 100644
--- a/ogstools/meshplotlib/core.py
+++ b/ogstools/meshplotlib/core.py
@@ -33,31 +33,6 @@ def _q_zero_line(property: Property, levels: 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 not in mesh.point_data
- and property.data_name not in mesh.cell_data
- ):
- msg = f"Property {property.data_name} not found in 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):
return np.array(
[
@@ -192,7 +167,7 @@ def subplot(
ax.axis("auto")
- if has_masked_values(mesh, property):
+ if property.mask_used(mesh):
subplot(mesh, property.get_mask(), ax)
mesh = mesh.ctp(True).threshold(value=[1, 1], scalars=property.mask)
@@ -205,9 +180,9 @@ def subplot(
# faces contains a padding indicating number of points per face which gets
# removed with this reshaping and slicing to get the array of tri's
- x, y = setup.length.strip_units(surf_tri.points.T[[x_id, y_id]])
+ x, y = setup.length(surf_tri.points.T[[x_id, y_id]])
tri = surf_tri.faces.reshape((-1, 4))[:, 1:]
- values = property.magnitude.strip_units(get_data(surf_tri, property))
+ values = property.magnitude(surf_tri)
if setup.log_scaled:
values_temp = np.where(values > 1e-14, values, 1e-14)
values = np.log10(values_temp)
@@ -333,9 +308,7 @@ 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, masked=has_masked_values(mesh, property))
- )
+ values = property.magnitude(mesh)
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
@@ -463,7 +436,9 @@ def plot_diff(
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_unit = (
+ data_property(1, strip_unit=False) - data_property(1, strip_unit=False)
+ ).units
diff_property = property.replace(
data_unit=diff_unit,
output_unit=diff_unit,
@@ -538,7 +513,7 @@ def plot_probe(
if isinstance(mesh_property, str):
data_shape = mesh_series.read(0)[mesh_property].shape
mesh_property = _resolve_property(mesh_property, data_shape)
- values = mesh_property.magnitude.strip_units(
+ values = mesh_property.magnitude(
mesh_series.probe(
points, mesh_property.data_name, interp_method, interp_backend_pvd
)
@@ -556,7 +531,7 @@ def plot_probe(
mesh_property_abscissa = _resolve_property(
mesh_property_abscissa, data_shape
)
- x_values = mesh_property_abscissa.magnitude.strip_units(
+ x_values = mesh_property_abscissa.magnitude(
mesh_series.probe(
points,
mesh_property_abscissa.data_name,
diff --git a/ogstools/meshplotlib/plot_features.py b/ogstools/meshplotlib/plot_features.py
index 9345093d32b524b053b799d2e884a1b91fb3d007..8d1baeca82420a538fafad6cbf432e62e67828e9 100644
--- a/ogstools/meshplotlib/plot_features.py
+++ b/ogstools/meshplotlib/plot_features.py
@@ -10,7 +10,7 @@ from matplotlib import patheffects
from matplotlib.collections import PolyCollection
from matplotlib.transforms import blended_transform_factory as btf
-from ogstools.propertylib.property import Vector
+from ogstools.propertylib import Vector
from . import setup
@@ -28,7 +28,7 @@ def plot_layer_boundaries(
for reg_id in np.unique(segments.cell_data["RegionId"]):
segment = segments.threshold((reg_id, reg_id), "RegionId")
edges = segment.extract_surface().strip(True, 10000)
- x_b, y_b = setup.length.strip_units(
+ x_b, y_b = setup.length(
edges.points[edges.lines % edges.n_points].T[[x_id, y_id]]
)
lw = setup.rcParams_scaled["lines.linewidth"]
@@ -72,7 +72,7 @@ def plot_element_edges(ax: plt.Axes, surf: pv.DataSet, projection: int) -> None:
cell_pts = [
cp for cp, ct in zip(cell_points, cell_types) if ct == cell_type
]
- verts = setup.length.strip_units(np.delete(cell_pts, projection, -1))
+ verts = setup.length(np.delete(cell_pts, projection, -1))
lw = 0.5 * setup.rcParams_scaled["lines.linewidth"]
pc = PolyCollection(verts, fc="None", ec="black", lw=lw)
ax.add_collection(pc)
@@ -103,7 +103,7 @@ def plot_streamlines(
grid.point_data[property.data_name], projection, 1
)
val = np.reshape(
- property.strip_units(grid.point_data[property.data_name]),
+ property(grid.point_data[property.data_name]),
(n_pts, n_pts, 2),
)
@@ -114,7 +114,7 @@ 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.strip_units(np.meshgrid(x, y))
+ x_g, y_g = setup.length(np.meshgrid(x, y))
ax.streamplot(x_g, y_g, val[..., 0], val[..., 1], color="k", linewidth=lw)
@@ -134,11 +134,11 @@ 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.strip_units(np.max(contour_vals)))
+ ax.set_ylim(top=setup.length(np.max(contour_vals)))
ax.fill_between(
- setup.length.strip_units(x_vals),
- setup.length.strip_units(y_vals),
- setup.length.strip_units(contour_vals),
+ setup.length(x_vals),
+ setup.length(y_vals),
+ setup.length(contour_vals),
facecolor="lightgrey",
)
diff --git a/ogstools/propertylib/__init__.py b/ogstools/propertylib/__init__.py
index 63a4235764384604326c80826d92cdaddcdd2110..6f954c171c6119757556472c84217279a8c4e7f0 100644
--- a/ogstools/propertylib/__init__.py
+++ b/ogstools/propertylib/__init__.py
@@ -1,10 +1,14 @@
# Author: Florian Zill (Helmholtz Centre for Environmental Research GmbH - UFZ)
"""Define easy-to-access Property classes and PropertyCollection instances."""
-from . import presets
-from .property import Matrix, Property, Scalar, Vector
+from . import mesh_dependent, presets, tensor_math
+from .matrix import Matrix
+from .property import Property, Scalar
+from .vector import Vector
__all__ = [
+ "tensor_math",
+ "mesh_dependent",
"presets",
"Property",
"Scalar",
diff --git a/ogstools/propertylib/custom_colormaps.py b/ogstools/propertylib/custom_colormaps.py
index fd101559fe1fed313f68af186a1384dfccafcf97..c91712f241097e9aad3d3a36279d55a044714b8b 100644
--- a/ogstools/propertylib/custom_colormaps.py
+++ b/ogstools/propertylib/custom_colormaps.py
@@ -8,7 +8,14 @@ mask_cmap = LSC.from_list("mask_cmap", ["lightgrey", "g"])
temperature_cmap = LSC.from_list(
"temperature_cmap",
np.r_[
- colormaps["Blues"](np.linspace(0, 0.75, 128, endpoint=True)),
- colormaps["plasma"](np.linspace(0, 1, 128, endpoint=True)),
+ colormaps["Blues"](np.linspace(0, 0.75, 128)),
+ colormaps["plasma"](np.linspace(0, 1, 128)),
+ ],
+)
+integrity_cmap = LSC.from_list(
+ "criterion_cmap",
+ np.r_[
+ colormaps["viridis"](np.linspace(0, 0.75, 128)),
+ colormaps["autumn_r"](np.linspace(0.0, 1.0, 128)),
],
)
diff --git a/ogstools/propertylib/matrix.py b/ogstools/propertylib/matrix.py
new file mode 100644
index 0000000000000000000000000000000000000000..d8ad499a181e7dfcc6c671f2e8c76e61e2767489
--- /dev/null
+++ b/ogstools/propertylib/matrix.py
@@ -0,0 +1,190 @@
+from dataclasses import dataclass
+from typing import Literal, Union
+
+from ogstools.propertylib import tensor_math
+from ogstools.propertylib.property import Property, Scalar
+from ogstools.propertylib.vector import Vector, VectorList
+
+
+@dataclass
+class Matrix(Property):
+ """Represent a matrix property of a dataset.
+
+ Matrix properties should contain either 4 (2D) or 6 (3D) components.
+ Matrix components can be accesses with brackets e.g. stress[0]
+ """
+
+ def __getitem__(
+ self, index: Union[int, Literal["xx", "yy", "zz", "xy", "yz", "xz"]]
+ ) -> Scalar:
+ "A scalar property as a matrix component."
+ int_index = (
+ index
+ if isinstance(index, int)
+ else ["xx", "yy", "zz", "xy", "yz", "xz"].index(index)
+ )
+ return Scalar.from_property(
+ self,
+ output_name=self.output_name + f"_{index}",
+ func=lambda x: self.func(x)[..., int_index],
+ bilinear_cmap=True,
+ )
+
+ @property
+ def magnitude(self) -> Scalar:
+ "A scalar property as the frobenius norm of the matrix."
+ return Scalar.from_property(
+ self,
+ output_name=self.output_name + "_magnitude",
+ func=lambda x: tensor_math.frobenius_norm(self.func(x)),
+ )
+
+ @property
+ def trace(self) -> Scalar:
+ "A scalar property as the trace of the matrix."
+ return Scalar.from_property(
+ self,
+ output_name=self.output_name + "_trace",
+ func=tensor_math.trace,
+ )
+
+ @property
+ def eigenvalues(self) -> Vector:
+ "A vector property as the eigenvalues of the matrix."
+ return Vector.from_property(
+ self,
+ output_name=self.output_name + "_eigenvalues",
+ func=lambda x: tensor_math.eigenvalues(self.func(x)),
+ )
+
+ @property
+ def eigenvectors(self) -> VectorList:
+ "A vector property as the eigenvectors of the matrix."
+ return VectorList.from_property(
+ self,
+ output_name=self.output_name + "_eigenvectors",
+ data_unit="",
+ output_unit="",
+ func=lambda x: tensor_math.eigenvectors(self.func(x)),
+ )
+
+ @property
+ def det(self) -> Scalar:
+ "A scalar property as the determinant of the matrix."
+ return Scalar.from_property(
+ self,
+ output_name=self.output_name + "_det",
+ func=lambda x: tensor_math.det(self.func(x)),
+ )
+
+ @property
+ def I1(self) -> Scalar: # pylint: disable=invalid-name
+ "A scalar property as the first invariant of the matrix."
+ return Scalar.from_property(
+ self,
+ output_name=self.output_name + "_I1",
+ func=lambda x: tensor_math.I1(self.func(x)),
+ )
+
+ @property
+ def I2(self) -> Scalar: # pylint: disable=invalid-name
+ "A scalar property as the second invariant of the matrix."
+ return Scalar.from_property(
+ self,
+ output_unit=self.output_unit + "^2",
+ output_name=self.output_name + "_I2",
+ func=lambda x: tensor_math.I2(self.func(x)),
+ process_with_units=True,
+ )
+
+ @property
+ def I3(self) -> Scalar: # pylint: disable=invalid-name
+ "A scalar property as the third invariant of the matrix."
+ return Scalar.from_property(
+ self,
+ output_name=self.output_name + "_I3",
+ func=lambda x: tensor_math.I3(self.func(x)),
+ )
+
+ @property
+ def mean(self) -> Scalar:
+ "A scalar property as the mean value of the matrix."
+ return Scalar.from_property(
+ self,
+ output_name="mean_" + self.output_name,
+ func=lambda x: tensor_math.mean(self.func(x)),
+ )
+
+ @property
+ def hydrostatic_component(self) -> "Matrix":
+ "A vector property as the effective pressure of the matrix."
+ return Matrix.from_property(
+ self,
+ output_name="hydrostatic_" + self.output_name + "_component",
+ func=lambda x: tensor_math.hydrostatic_component(self.func(x)),
+ )
+
+ @property
+ def deviator(self) -> "Matrix":
+ "A vector property as the deviator of the matrix."
+ return Matrix.from_property(
+ self,
+ output_name=self.output_name + "_deviator",
+ func=lambda x: tensor_math.deviator(self.func(x)),
+ )
+
+ @property
+ def J1(self) -> Scalar: # pylint: disable=invalid-name
+ "A scalar property as the first invariant of the matrix deviator."
+ return Scalar.from_property(
+ self,
+ output_name=self.output_name + "_J1",
+ func=lambda x: tensor_math.J1(self.func(x)),
+ )
+
+ @property
+ def J2(self) -> Scalar: # pylint: disable=invalid-name
+ "A scalar property as the second invariant of the matrix deviator."
+ return Scalar.from_property(
+ self,
+ output_name=self.output_name + "_J2",
+ func=lambda x: tensor_math.J2(self.func(x)),
+ )
+
+ @property
+ def J3(self) -> Scalar: # pylint: disable=invalid-name
+ "A scalar property as the third invariant of the matrix deviator."
+ return Scalar.from_property(
+ self,
+ output_name=self.output_name + "_J3",
+ func=lambda x: tensor_math.J3(self.func(x)),
+ )
+
+ @property
+ def octahedral_shear(self) -> Scalar:
+ "A scalar property as the octahedral shear component of the matrix."
+ return Scalar.from_property(
+ self,
+ output_name="octahedral_shear_" + self.output_name,
+ func=lambda x: tensor_math.octahedral_shear(self.func(x)),
+ )
+
+ @property
+ def von_Mises(self) -> Scalar:
+ "A scalar property as the von Mises stress."
+ return Scalar.from_property(
+ self,
+ output_name="von_Mises_" + self.output_name,
+ func=lambda x: tensor_math.von_mises(self.func(x)),
+ )
+
+ @property
+ def qp_ratio(self) -> Scalar:
+ "A scalar property as the qp stress ratio."
+ return Scalar.from_property(
+ self,
+ output_name="qp_ratio",
+ output_unit="percent",
+ func=lambda x: tensor_math.qp_ratio(self.func(x)),
+ process_with_units=True,
+ )
diff --git a/ogstools/propertylib/mesh_dependent.py b/ogstools/propertylib/mesh_dependent.py
new file mode 100644
index 0000000000000000000000000000000000000000..434a5a88e6842039a6cf002fb021f25845da62dc
--- /dev/null
+++ b/ogstools/propertylib/mesh_dependent.py
@@ -0,0 +1,104 @@
+"Functions related to stress analysis which can be only applied to a mesh."
+
+from functools import partial
+from typing import Union
+
+import numpy as np
+import pyvista as pv
+from pint.facets.plain import PlainQuantity
+
+from .property import Property
+from .tensor_math import eigenvalues, mean, octahedral_shear
+from .unit_registry import u_reg
+
+ValType = Union[PlainQuantity, np.ndarray]
+
+
+# TODO: find smart way to choose between y values and depth
+def depth(mesh: pv.UnstructuredGrid) -> np.ndarray:
+ # return -mesh.points[:, 1]
+ edges = mesh.extract_feature_edges()
+ cell_bounds = np.asarray([cell.bounds for cell in edges.cell])
+ cell_vec = np.abs(cell_bounds[..., 1:2] - cell_bounds[..., :2])
+ ids = np.argmax(cell_vec, axis=-1) == 0
+ ids2 = edges.cell_centers().points[:, 1] > edges.center[1]
+ top = edges.extract_cells(ids & ids2)
+
+ vert_diff = mesh.points[:, 1, None] - top.points[:, 1]
+ return np.abs(np.min(vert_diff, axis=-1))
+
+
+def p_fluid(mesh: pv.UnstructuredGrid) -> PlainQuantity:
+ """Return the fluid pressure in the mesh.
+
+ If not given, calculate by hypothetical water column."""
+ qty = u_reg.Quantity
+ if "pressure" in mesh.point_data:
+ return qty(mesh["pressure"], "Pa")
+ _depth = mesh["depth"] if "depth" in mesh.point_data else depth(mesh)
+ return qty(1000, "kg/m^3") * qty(9.81, "m/s^2") * qty(_depth, "m")
+
+
+def fluid_pressure_criterion(
+ mesh: pv.UnstructuredGrid, mesh_property: Property
+) -> PlainQuantity:
+ "Return the 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]
+
+
+def dilatancy_critescu(
+ mesh: pv.UnstructuredGrid,
+ mesh_property: Property,
+ a: float = -0.01697,
+ b: float = 0.8996,
+ effective: bool = False,
+) -> PlainQuantity:
+ """Return the dilatancy criterion defined as:
+
+ :math:`F_{dil} = \\frac{\\tau}{\\sigma_0} - a \\left( \\frac{\\sigma_m}{\\sigma_0} \\right)^2 - b \\frac{\\sigma_m}{\\sigma_0}`
+
+ <https://www.sciencedirect.com/science/article/pii/S0360544222000512?via%3Dihub>
+ """
+
+ qty = u_reg.Quantity
+ sigma = -qty(mesh[mesh_property.data_name], mesh_property.data_unit)
+ sigma_0 = qty(1, "MPa")
+ sigma_m = mean(sigma)
+ if effective:
+ sigma_m -= p_fluid(mesh)
+ tau_oct = octahedral_shear(sigma)
+ return (
+ tau_oct / sigma_0 - a * (sigma_m / sigma_0) ** 2 - b * sigma_m / sigma_0
+ )
+
+
+dilatancy_critescu_eff = partial(dilatancy_critescu, effective=True)
+
+
+def dilatancy_alkan(
+ mesh: pv.UnstructuredGrid,
+ mesh_property: Property,
+ b: float = 0.04,
+ effective: bool = False,
+) -> ValType:
+ """Return the dilatancy criterion defined as:
+
+ :math:`F_{dil} = \\tau - \\tau_{max} \\cdot b \\frac{\\sigma_m}{\\sigma_0 + b \\cdot \\sigma_m}`
+
+ <https://www.sciencedirect.com/science/article/pii/S1365160906000979>
+ """
+
+ qty = u_reg.Quantity
+ sigma = -qty(mesh[mesh_property.data_name], mesh_property.data_unit)
+ tau_max = qty(33, "MPa")
+ sigma_m = mean(sigma)
+ if effective:
+ sigma_m -= p_fluid(mesh)
+ tau = octahedral_shear(sigma)
+ return tau - tau_max * (b * sigma_m / (qty(1, "MPa") + b * sigma_m))
+
+
+dilatancy_alkan_eff = partial(dilatancy_alkan, effective=True)
diff --git a/ogstools/propertylib/presets.py b/ogstools/propertylib/presets.py
index 6bd7127284fc9317951a4b35ef1d58b2b562ba96..139faac875debf1d1ee142f1cd3480c04bdd578a 100644
--- a/ogstools/propertylib/presets.py
+++ b/ogstools/propertylib/presets.py
@@ -1,53 +1,168 @@
# flake8: noqa: E501
-"Predefined properties."
+"""Predefined properties.
+
+".. seealso:: :py:mod:`ogstools.propertylib.tensor_math`"
+"""
from typing import Optional
-from . import vector2scalar as v2s
-from .property import Matrix, Property, Scalar, Vector
-
-from .custom_colormaps import temperature_cmap
-
-T_mask = "temperature_active"
-H_mask = "pressure_active"
-M_mask = "displacement_active"
-
-# fmt: off
-material_id = Scalar("MaterialIDs", categoric=True, cmap="tab20")
-displacement = Vector("displacement", "m", "m", mask=M_mask, cmap="PRGn")
-effective_pressure = Scalar("sigma", "Pa", "MPa", "effective_pressure", M_mask, v2s.effective_pressure)
-heatflowrate = Scalar("HeatFlowRate", mask=T_mask)
-massflowrate = Scalar("MassFlowRate", mask=H_mask)
-nodal_forces = Vector("NodalForces", mask=M_mask)
-pressure = Scalar("pressure", "Pa", "MPa", "pore_pressure", H_mask, cmap="Blues")
-hydraulic_height = Scalar("pressure", "m", "m", "hydraulic_height", H_mask, cmap="Blues")
-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, cmap=temperature_cmap, 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
+from . import mesh_dependent, tensor_math
+from .custom_colormaps import integrity_cmap, temperature_cmap
+from .matrix import Matrix
+from .property import Property, Scalar
+from .vector import Vector
+
+T_MASK = "temperature_active"
+H_MASK = "pressure_active"
+M_MASK = "displacement_active"
+
+# general
+material_id = Scalar(data_name="MaterialIDs", categoric=True, cmap="tab20")
+
+# thermal
+temperature = Scalar(
+ data_name="temperature",
+ data_unit="K",
+ output_name="temperature",
+ output_unit="°C",
+ mask=T_MASK,
+ cmap=temperature_cmap,
+ bilinear_cmap=True,
+)
+heatflowrate = Scalar(data_name="HeatFlowRate", mask=T_MASK)
+
+# hydraulic
+pressure = Scalar(
+ data_name="pressure",
+ data_unit="Pa",
+ output_unit="MPa",
+ output_name="pore_pressure",
+ mask=H_MASK,
+ cmap="Blues",
+)
+hydraulic_height = Scalar(
+ data_name="pressure",
+ data_unit="m",
+ output_unit="m",
+ output_name="hydraulic_height",
+ mask=H_MASK,
+ cmap="Blues",
+)
+velocity = Vector(
+ data_name="velocity",
+ data_unit="m/s",
+ output_unit="m/s",
+ output_name="darcy_velocity",
+ mask=H_MASK,
+)
+massflowrate = Scalar(data_name="MassFlowRate", mask=H_MASK)
+
+# mechanical
+displacement = Vector(
+ data_name="displacement",
+ data_unit="m",
+ output_unit="m",
+ mask=M_MASK,
+ cmap="PRGn",
+ bilinear_cmap=True,
+)
+strain = Matrix(
+ data_name="epsilon",
+ data_unit="",
+ output_unit="percent",
+ output_name="strain",
+ mask=M_MASK,
+)
+stress = Matrix(
+ data_name="sigma",
+ data_unit="Pa",
+ output_unit="MPa",
+ output_name="stress",
+ mask=M_MASK,
+)
+effective_pressure = Scalar(
+ data_name="sigma",
+ data_unit="Pa",
+ output_unit="MPa",
+ output_name="effective_pressure",
+ mask=M_MASK,
+ func=tensor_math.effective_pressure,
+)
+dilatancy_critescu = Scalar(
+ data_name="sigma",
+ data_unit="Pa",
+ output_unit="",
+ output_name="dilatancy_criterion",
+ mask=M_MASK,
+ func=mesh_dependent.dilatancy_critescu,
+ mesh_dependent=True,
+ cmap=integrity_cmap,
+ bilinear_cmap=True,
+)
+dilatancy_critescu_eff = Scalar(
+ data_name="sigma",
+ data_unit="Pa",
+ output_unit="",
+ output_name="effective_dilatancy_criterion",
+ mask=M_MASK,
+ func=mesh_dependent.dilatancy_critescu_eff,
+ mesh_dependent=True,
+ cmap=integrity_cmap,
+ bilinear_cmap=True,
+)
+dilatancy_alkan = Scalar(
+ data_name="sigma",
+ data_unit="Pa",
+ output_unit="MPa",
+ output_name="dilatancy_criterion",
+ mask=M_MASK,
+ func=mesh_dependent.dilatancy_alkan,
+ mesh_dependent=True,
+ cmap=integrity_cmap,
+ bilinear_cmap=True,
+)
+dilatancy_alkan_eff = Scalar(
+ data_name="sigma",
+ data_unit="Pa",
+ output_unit="MPa",
+ output_name="effective_dilatancy_criterion",
+ mask=M_MASK,
+ func=mesh_dependent.dilatancy_alkan_eff,
+ mesh_dependent=True,
+ cmap=integrity_cmap,
+ bilinear_cmap=True,
+)
+fluid_pressure_crit = Scalar(
+ data_name="sigma",
+ data_unit="Pa",
+ output_unit="MPa",
+ output_name="fluid_pressure_criterion",
+ mask=M_MASK,
+ func=mesh_dependent.fluid_pressure_criterion,
+ mesh_dependent=True,
+ cmap=integrity_cmap,
+ bilinear_cmap=True,
+)
+nodal_forces = Vector(data_name="NodalForces", mask=M_MASK)
+
all_properties = [v for v in locals().values() if isinstance(v, Property)]
-def find_property_preset(property_name: str) -> Optional[Property]:
- """Return predefined property with given output_name."""
+def get_preset(property_name: str, shape: Optional[tuple] = None) -> Property:
+ """
+ Returns a Property preset or create one with correct type.
+
+ Searches for presets by data_name and output_name and returns if found.
+ Otherwise create Scalar, Vector or Matrix Property depending on shape.
+ """
for prop in all_properties:
if prop.output_name == property_name:
return prop
- # if not found by output name, find by data_name
for prop in all_properties:
if prop.data_name == property_name:
return prop
- return None
-
-
-def _resolve_property(property_name: str, shape: tuple) -> Property:
- if found_property := find_property_preset(property_name):
- return found_property
- if len(shape) == 1:
+ if shape is None or len(shape) == 1:
return Scalar(property_name)
if shape[1] in [2, 3]:
return Vector(property_name)
diff --git a/ogstools/propertylib/property.py b/ogstools/propertylib/property.py
index 711c06da2d11e07fc439e0ee8cb0845e8c182008..fe83d01dfc42c5ef893aad41b29f3fb77b5ea5b7 100644
--- a/ogstools/propertylib/property.py
+++ b/ogstools/propertylib/property.py
@@ -5,20 +5,19 @@ way (e.g. temperature, pressure, displacement, …). Unit conversion is handled
via pint.
"""
+from collections.abc import Sequence
from dataclasses import dataclass, replace
-from typing import Any, Callable, Union
+from typing import 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
from .unit_registry import u_reg
-from .utils import identity, sym_tensor_to_mat
-from .vector2scalar import trace
-# TODO: rename to BaseProperty?? / GenericProperty
@dataclass
class Property:
"""Represent a property of a dataset."""
@@ -33,14 +32,12 @@ class Property:
"""The output name of the property."""
mask: str = ""
"""The name of the mask data in the dataset."""
- func: Union[
- Callable[
- [Union[float, np.ndarray, PlainQuantity]],
- Union[float, np.ndarray, PlainQuantity],
- ],
- Callable[[Any], Any],
- ] = identity
+ func: Callable = identity
"""The function to be applied on the data."""
+ mesh_dependent: bool = False
+ """If the function to be applied is dependent on the mesh itself"""
+ process_with_units: bool = False
+ """If true, apply the function on values with units."""
cmap: Union[Colormap, str] = "coolwarm"
"""Colormap to use for plotting."""
bilinear_cmap: bool = False
@@ -69,33 +66,52 @@ class Property:
"""
return replace(self, **changes)
- def __call__(self, vals: np.ndarray) -> PlainQuantity:
- """
- Return transformed values with units.
-
- Apply property function and convert from data_unit to output_unit
-
- :param vals: The input values.
-
- :returns: The values with units.
- """
- Q_, _du, _ou = u_reg.Quantity, self.data_unit, self.output_unit
- if Q_(0, _du).dimensionality == Q_(0, _ou).dimensionality:
- return Q_(Q_(self.func(np.asarray(vals)), _du), _ou)
- return Q_(self.func(Q_(vals, _du)), _ou)
-
- def strip_units(self, vals: np.ndarray) -> np.ndarray:
- """
- Return transformed values without units.
+ @classmethod
+ def from_property(cls, new_property: "Property", **changes):
+ "Create a new Property object with modified attributes."
+ return cls(
+ data_name=new_property.data_name,
+ data_unit=new_property.data_unit,
+ output_unit=new_property.output_unit,
+ output_name=new_property.output_name,
+ mask=new_property.mask,
+ func=new_property.func,
+ mesh_dependent=new_property.mesh_dependent,
+ process_with_units=new_property.process_with_units,
+ cmap=new_property.cmap,
+ categoric=new_property.categoric,
+ ).replace(**changes)
+
+ def __call__(
+ self,
+ data: Union[int, float, np.ndarray, pv.DataSet, Sequence],
+ strip_unit: bool = True,
+ ) -> np.ndarray:
+ """
+ Return the transformed data values.
Apply property function, convert from data_unit to output_unit and
- strip the unit.
-
- :param vals: The input values.
-
- :returns: The values without units.
- """
- return self(vals).magnitude
+ return the values, optionally with the unit.
+
+ :param vals: The input data values.
+
+ :returns: The transformed data values.
+ """
+ qty, _du, _ou = u_reg.Quantity, self.data_unit, self.output_unit
+ if self.mesh_dependent:
+ if isinstance(data, pv.DataSet):
+ result = qty(self.func(data, self), _ou)
+ else:
+ msg = "This property can only be evaluated on a mesh."
+ raise TypeError(msg)
+ else:
+ if isinstance(data, pv.DataSet):
+ result = qty(self.func(qty(self.get_data(data), _du)), _ou)
+ elif self.process_with_units:
+ result = qty(self.func(qty(data, _du)), _ou)
+ else:
+ result = qty(qty(self.func(np.asarray(data)), _du), _ou)
+ return result.magnitude if strip_unit else result
def get_output_unit(self) -> str:
"""
@@ -125,114 +141,30 @@ class Property:
def magnitude(self) -> "Property":
return self
-
-@dataclass
-class Scalar(Property):
- "Represent a scalar property of a dataset."
-
-
-@dataclass
-class Vector(Property):
- """Represent a vector property of a dataset.
-
- Vector properties should contain either 2 (2D) or 3 (3D) components.
- Vector components can be accesses with brackets e.g. displacement[0]
- """
-
- def __getitem__(self, index: int) -> Scalar:
- """
- Get a scalar property as a specific component of the vector property.
-
- :param index: The index of the component.
-
- :returns: A scalar property as a vector component.
- """
- suffix = {False: index, True: ["x", "y", "z"][index]}
- return Scalar(
- data_name=self.data_name,
- data_unit=self.data_unit,
- output_unit=self.output_unit,
- output_name=self.output_name + f"_{suffix[0 <= index <= 2]}",
- mask=self.mask,
- func=lambda x: np.array(x)[..., index],
- cmap=self.cmap,
- bilinear_cmap=True,
- )
-
- @property
- def magnitude(self) -> Scalar:
- ":returns: A scalar property as the magnitude of the vector."
- return Scalar(
- data_name=self.data_name,
- data_unit=self.data_unit,
- output_unit=self.output_unit,
- output_name=self.output_name + "_magnitude",
- mask=self.mask,
- func=lambda x: np.linalg.norm(x, axis=-1),
- cmap=self.cmap,
+ def mask_used(self, mesh: pv.UnstructuredGrid) -> bool:
+ return (
+ not self.is_mask()
+ and self.mask in mesh.cell_data
+ and (len(mesh.cell_data[self.mask]) != 0)
)
- @property
- def log_magnitude(self) -> Scalar:
- ":returns: A scalar property as the log-magnitude of the vector."
- return Scalar(
- data_name=self.data_name,
- output_name=self.output_name + "_log10",
- mask=self.mask,
- func=lambda x: np.log10(np.linalg.norm(x, axis=-1)),
- cmap=self.cmap,
- )
+ def get_data(
+ self, mesh: pv.UnstructuredGrid, masked: bool = True
+ ) -> pv.UnstructuredGrid:
+ """Get the data associated with a scalar or vector property from a mesh."""
+ if (
+ self.data_name not in mesh.point_data
+ and self.data_name not in mesh.cell_data
+ ):
+ msg = f"Property {self.data_name} not found in mesh."
+ raise IndexError(msg)
+ if masked and self.mask_used(mesh):
+ return mesh.ctp(True).threshold(value=[1, 1], scalars=self.mask)[
+ self.data_name
+ ]
+ return mesh[self.data_name]
@dataclass
-class Matrix(Property):
- """Represent a matrix property of a dataset.
-
- Matrix properties should contain either 4 (2D) or 6 (3D) components.
- Matrix components can be accesses with brackets e.g. stress[0]
- """
-
- def __getitem__(self, index: int) -> Scalar:
- """
- Get a scalar property as a specific component of the matrix property.
-
- :param index: The index of the component.
-
- :returns: A scalar property as a matrix component.
- """
- suffix = {False: index, True: ["x", "y", "z", "xy", "yz", "xz"][index]}
- return Scalar(
- data_name=self.data_name,
- data_unit=self.data_unit,
- output_unit=self.output_unit,
- output_name=self.output_name + f"_{suffix[0 <= index <= 5]}",
- mask=self.mask,
- func=lambda x: np.array(x)[..., index],
- bilinear_cmap=True,
- cmap=self.cmap,
- )
-
- @property
- def magnitude(self) -> Scalar:
- ":returns: A scalar property as the frobenius norm of the matrix."
- return Scalar(
- data_name=self.data_name,
- data_unit=self.data_unit,
- output_unit=self.output_unit,
- output_name=self.output_name + "_magnitude",
- mask=self.mask,
- func=lambda x: np.linalg.norm(sym_tensor_to_mat(x), axis=(-2, -1)),
- cmap=self.cmap,
- )
-
- @property
- def trace(self) -> Scalar:
- ":returns: A scalar property as the trace of the matrix."
- return Scalar(
- data_name=self.data_name,
- data_unit=self.data_unit,
- output_unit=self.output_unit,
- output_name=self.output_name + "_trace",
- mask=self.mask,
- func=trace,
- )
+class Scalar(Property):
+ "Represent a scalar property of a dataset."
diff --git a/ogstools/propertylib/tensor_math.py b/ogstools/propertylib/tensor_math.py
new file mode 100644
index 0000000000000000000000000000000000000000..e58b6b789ef2007514d2738d23d18c0072cadf53
--- /dev/null
+++ b/ogstools/propertylib/tensor_math.py
@@ -0,0 +1,210 @@
+"""Common tensor transformation operations.
+
+They can be used as they are, but are also part of
+:py:obj:`ogstools.propertylib.matrix.Matrix` properties.
+All input arrays are expected to be in vector notation of symmetric tensors in
+the form of:
+
+[xx, yy, zz, xy] for 2D and
+
+[xx, yy, zz, xy, yz, xz] for 3D.
+
+This notation style is the default output of OGS:
+
+<https://www.opengeosys.org/docs/userguide/basics/conventions/#a-namesymmetric-tensorsa--symmetric-tensors-and-kelvin-mapping>
+
+A better overview for the theoretical background of the equations can be found
+here, for example:
+
+<https://en.wikipedia.org/wiki/Cauchy_stress_tensor#Cauchy's_stress_theorem%E2%80%94stress_tensor>
+"""
+
+from typing import TypeVar, Union
+
+import numpy as np
+from numpy import linalg as LA
+from pint.facets.plain import PlainQuantity
+
+from .unit_registry import u_reg
+
+ValType = Union[PlainQuantity, np.ndarray]
+
+
+T = TypeVar("T")
+
+
+def identity(vals: T) -> T:
+ ":returns: The input values."
+ return vals
+
+
+def sym_tensor_to_mat(vals: np.ndarray) -> np.ndarray:
+ "Convert an symmetric tensor to a 3x3 matrix."
+ assert np.shape(vals)[-1] in [4, 6]
+ shape = list(np.shape(vals))[:-1] + [3, 3]
+ mat = np.zeros(shape)
+ idx = {0: [0, 0], 1: [1, 1], 2: [2, 2], 3: [0, 1], 4: [1, 2], 5: [0, 2]}
+ for i in range(np.shape(vals)[-1]):
+ mat[..., idx[i][0], idx[i][1]] = vals[..., i]
+ return mat
+
+
+def trace(vals: ValType) -> ValType:
+ """Return the trace.
+
+ :math:`tr(\\mathbf{\\sigma}) = \\sum\\limits_{i=1}^3 \\sigma_{ii}`
+ """
+ return np.sum(vals[..., :3], axis=-1)
+
+
+def eigenvalues(vals: ValType) -> ValType:
+ "Return the eigenvalues."
+ if isinstance(vals, PlainQuantity):
+ unit = vals.units
+ vals = vals.magnitude
+ else:
+ unit = None
+ eigvals: np.ndarray = LA.eigvals(sym_tensor_to_mat(vals))
+ eigvals.sort(axis=-1)
+ assert np.all(eigvals[..., 0] <= eigvals[..., 1])
+ assert np.all(eigvals[..., 1] <= eigvals[..., 2])
+ return eigvals if unit is None else u_reg.Quantity(eigvals, unit)
+
+
+def eigenvectors(vals: ValType) -> ValType:
+ "Return the eigenvectors."
+ if isinstance(vals, PlainQuantity):
+ vals = vals.magnitude
+ eigvals, eigvecs = LA.eig(sym_tensor_to_mat(vals))
+ ids = eigvals.argsort(axis=-1)
+ eigvals = np.take_along_axis(eigvals, ids, axis=-1)
+ eigvecs = np.take_along_axis(eigvecs, ids[:, np.newaxis], axis=-1)
+ assert np.all(eigvals[..., 0] <= eigvals[..., 1])
+ assert np.all(eigvals[..., 1] <= eigvals[..., 2])
+ return eigvecs
+
+
+def det(vals: ValType) -> ValType:
+ "Return the determinants."
+ if isinstance(vals, PlainQuantity):
+ unit = vals.units
+ vals = vals.magnitude
+ else:
+ unit = None
+ result = np.linalg.det(sym_tensor_to_mat(vals))
+ return result if unit is None else u_reg.Quantity(result, unit)
+
+
+def frobenius_norm(val: ValType) -> ValType:
+ """Return the Frobenius norm.
+
+ :math:`||\\mathbf{\\sigma}||_F = \\sqrt{ \\sum\\limits_{i=1}^m \\sum\\limits_{j=1}^n |\\sigma_{ij}|^2 }`
+ """
+ return np.linalg.norm(sym_tensor_to_mat(val), axis=(-2, -1))
+
+
+def I1(vals: ValType) -> ValType: # pylint: disable=invalid-name
+ """Return the first invariant.
+
+ :math:`I1 = tr(\\mathbf{\\sigma})`
+ """
+ return trace(vals)
+
+
+def I2(vals: ValType) -> ValType: # pylint: disable=invalid-name
+ """Return the second invariant.
+
+ :math:`I2 = \\frac{1}{2} \\left[(tr(\\mathbf{\\sigma}))^2 - tr(\\mathbf{\\sigma}^2) \\right]`
+ """
+ return 0.5 * (trace(vals) ** 2 - trace(vals**2))
+
+
+def I3(vals: ValType) -> ValType: # pylint: disable=invalid-name
+ """Return the third invariant.
+
+ :math:`I3 = det(\\mathbf{\\sigma})`
+ """
+ return det(vals)
+
+
+def mean(vals: ValType) -> ValType:
+ """Return the mean value.
+ Also called hydrostatic component or octahedral normal component.
+
+ :math:`\\pi = \\frac{1}{3} I1`
+ """
+ return (1.0 / 3.0) * I1(vals)
+
+
+def effective_pressure(vals: ValType) -> ValType:
+ """Return the effective pressure.
+
+ :math:`\\pi = -\\frac{1}{3} I1`
+ """
+ return -mean(vals)
+
+
+def hydrostatic_component(vals: ValType) -> ValType:
+ """Return the hydrostatic component.
+
+ :math:`p_{ij} = \\pi \\delta_{ij}`
+ """
+ result = vals * 0.0
+ result[..., :3] = np.tile(mean(vals), (3, 1)).T
+ return result
+
+
+def deviator(vals: ValType) -> ValType:
+ """Return the deviator.
+
+ :math:`s_{ij} = \\sigma_{ij} - \\pi \\delta_{ij}`
+ """
+ return vals - hydrostatic_component(vals)
+
+
+def J1(vals: ValType) -> ValType: # pylint: disable=invalid-name
+ """Return the first invariant of the deviator.
+
+ :math:`J1 = 0`
+ """
+ return trace(deviator(vals))
+
+
+def J2(vals: ValType) -> ValType: # pylint: disable=invalid-name
+ """Return the second invariant of the deviator.
+
+ :math:`J2 = \\frac{1}{2} tr(\\mathbf{s}^2)`
+ """
+ return 0.5 * trace(deviator(vals) ** 2)
+
+
+def J3(vals: ValType) -> ValType: # pylint: disable=invalid-name
+ """Return the third invariant of the deviator.
+
+ :math:`J3 = \\frac{1}{3} tr(\\mathbf{s}^3)`
+ """
+ return (1.0 / 3.0) * trace(deviator(vals) ** 3)
+
+
+def octahedral_shear(vals: ValType) -> ValType:
+ """Return the octahedral shear value.
+
+ :math:`\\tau_{oct} = \\sqrt{\\frac{2}{3} J2}`
+ """
+ return np.sqrt((2.0 / 3.0) * J2(vals))
+
+
+def von_mises(vals: ValType) -> ValType:
+ """Return the von Mises stress.
+
+ :math:`\\sigma_{Mises} = \\sqrt{3 J2}`
+ """
+ return np.sqrt(3.0 * J2(vals))
+
+
+def qp_ratio(vals: ValType) -> ValType:
+ """Return the qp ratio (von Mises stress / effective pressure).
+
+ :math:`qp = \\sigma_{Mises} / \\pi`
+ """
+ return von_mises(vals) / effective_pressure(vals)
diff --git a/ogstools/propertylib/utils.py b/ogstools/propertylib/utils.py
deleted file mode 100644
index 3b2c1cc09cf8ae82d2fff9bc29451aaa92ec8ce8..0000000000000000000000000000000000000000
--- a/ogstools/propertylib/utils.py
+++ /dev/null
@@ -1,30 +0,0 @@
-from typing import TypeVar
-
-import numpy as np
-
-
-def dim_from_len(len: int):
- ":returns: The dimension corresponding to the length. (2|4 -> 2, 3|6 -> 3)"
- dim_map = {2: 2, 3: 3, 4: 2, 6: 3}
- if len in dim_map:
- return dim_map[len]
- raise ValueError("Can't determine dimension for length " + str(len))
-
-
-def sym_tensor_to_mat(vals: np.ndarray) -> np.ndarray:
- "Convert an symmetric tensor to a 3x3 matrix."
- assert np.shape(vals)[-1] in [4, 6]
- shape = list(np.shape(vals))[:-1] + [3, 3]
- mat = np.zeros(shape)
- idx = {0: [0, 0], 1: [1, 1], 2: [2, 2], 3: [0, 1], 4: [1, 2], 5: [0, 2]}
- for i in range(np.shape(vals)[-1]):
- mat[..., idx[i][0], idx[i][1]] = vals[..., i]
- return mat
-
-
-T = TypeVar("T")
-
-
-def identity(vals: T) -> T:
- ":returns: The input values."
- return vals
diff --git a/ogstools/propertylib/vector.py b/ogstools/propertylib/vector.py
new file mode 100644
index 0000000000000000000000000000000000000000..a5a2b73dc2c2fc327335955f36a54a4a2fba4ec3
--- /dev/null
+++ b/ogstools/propertylib/vector.py
@@ -0,0 +1,71 @@
+from dataclasses import dataclass
+from typing import Literal, Union
+
+import numpy as np
+from pint.facets.plain import PlainQuantity
+
+from ogstools.propertylib.property import Property, Scalar
+
+from .unit_registry import u_reg
+
+ValType = Union[PlainQuantity, np.ndarray]
+
+
+def vector_norm(vals: ValType) -> ValType:
+ "Return the magnitudes."
+ if isinstance(vals, PlainQuantity):
+ unit = vals.units
+ vals = vals.magnitude
+ else:
+ unit = None
+ result = np.linalg.norm(vals, axis=-1)
+ return result if unit is None else u_reg.Quantity(result, unit)
+
+
+@dataclass
+class Vector(Property):
+ """Represent a vector property of a dataset.
+
+ Vector properties should contain either 2 (2D) or 3 (3D) components.
+ Vector components can be accesses with brackets e.g. displacement[0]
+ """
+
+ def __getitem__(self, index: Union[int, Literal["x", "y", "z"]]) -> Scalar:
+ """
+ Get a scalar property as a specific component of the vector property.
+
+ :param index: The index of the component.
+
+ :returns: A scalar property as a vector component.
+ """
+ int_index = index if isinstance(index, int) else "xyz".index(index)
+ return Scalar.from_property(
+ self,
+ output_name=self.output_name + f"_{index}",
+ func=lambda x: self.func(x)[..., int_index],
+ bilinear_cmap=True,
+ )
+
+ @property
+ def magnitude(self) -> Scalar:
+ ":returns: A scalar property as the magnitude of the vector."
+ return Scalar.from_property(
+ self,
+ output_name=self.output_name + "_magnitude",
+ func=lambda x: vector_norm(self.func(x)),
+ )
+
+
+@dataclass
+class VectorList(Property):
+ """Represent a list of vector properties of a dataset."""
+
+ def __getitem__(self, index: int) -> Vector:
+ """
+ Get a vector property as a component of the vectorlist property.
+ """
+ return Vector.from_property(
+ self,
+ output_name=self.output_name + f"_{index}",
+ func=lambda x: np.take(self.func(x), index, axis=-1),
+ )
diff --git a/ogstools/propertylib/vector2scalar.py b/ogstools/propertylib/vector2scalar.py
deleted file mode 100644
index f441f5d7e39783a3d84e346fd8b5b526768f0ac8..0000000000000000000000000000000000000000
--- a/ogstools/propertylib/vector2scalar.py
+++ /dev/null
@@ -1,53 +0,0 @@
-"""Common engineering transformation functions."""
-
-import numpy as np
-
-from .utils import dim_from_len
-
-
-def trace(vals: np.ndarray) -> np.ndarray:
- """
- Calculate the trace of each vector in the input array.
-
- :param vals: The input array of vectors.
-
- :returns: The trace values of the vectors.
- """
- return np.sum(vals[..., : dim_from_len(vals.shape[-1])], axis=-1)
-
-
-def effective_pressure(vals: np.ndarray) -> np.ndarray:
- """
- Calculate the effective pressure based on the input array.
-
- :param vals: The input array.
-
- :returns: The effective pressure values.
- """
- return -(1.0 / 3.0) * trace(vals)
-
-
-def von_mises(vals: np.ndarray) -> np.ndarray:
- """
- Calculate the von Mises stress based on the input array.
-
- :param vals: The input array.
-
- :returns: The von Mises stress values.
- """
- return np.sqrt(
- 0.5
- * np.sum(np.square(np.diff(vals[..., :3], append=vals[..., :1])), -1)
- + 3 * np.sum(np.square(vals[..., 3:]), -1)
- )
-
-
-def qp_ratio(vals: np.ndarray) -> np.ndarray:
- """
- Calculate the QP ratio (von Mises stress / effective pressure).
-
- :param vals: The input array.
-
- :returns: The QP ratios.
- """
- return von_mises(vals) / effective_pressure(vals)
diff --git a/ogstools/studies/convergence/convergence.py b/ogstools/studies/convergence/convergence.py
index 5734f70efac3eaf83cc0ae77ca528e70e2ea7b5d..ba79673d6cb5e6d39410c18200488dc9f98abb77 100644
--- a/ogstools/studies/convergence/convergence.py
+++ b/ogstools/studies/convergence/convergence.py
@@ -147,7 +147,7 @@ def convergence_metrics(
"""
def _data(m: pv.DataSet):
- return property.magnitude.strip_units(m.point_data[property.data_name])
+ return property.magnitude(m.point_data[property.data_name])
grid_spacings = [
np.mean(add_grid_spacing(mesh)["grid_spacing"]) for mesh in meshes
diff --git a/ogstools/studies/templates/convergence_study.py b/ogstools/studies/templates/convergence_study.py
index bb9c0a5f6f2017a8b64b4a38e97dd967f96c7a3b..8afbefb91fb31b6d9f8dcdcfbf35c5dfc28ee773 100644
--- a/ogstools/studies/templates/convergence_study.py
+++ b/ogstools/studies/templates/convergence_study.py
@@ -46,7 +46,7 @@ data_shape = (
if property_name in meshes[0].point_data
else None
)
-mesh_property = propertylib.presets._resolve_property(property_name, data_shape)
+mesh_property = propertylib.presets.get_preset(property_name, data_shape)
richardson = studies.convergence.richardson_extrapolation(
meshes, mesh_property, topology, refinement_ratio
)
diff --git a/tests/test_propertylib.py b/tests/test_propertylib.py
index d173c4da9f5d719c7a77c37eeb14d616b542ddee..4568e46d20bf605c06801557bc6c1975c779242c 100644
--- a/tests/test_propertylib.py
+++ b/tests/test_propertylib.py
@@ -8,7 +8,7 @@ from pint.facets.plain import PlainQuantity
from ogstools.propertylib import presets as pp
from ogstools.propertylib.property import Scalar, u_reg
-Q_ = u_reg.Quantity
+qty = u_reg.Quantity
class PhysicalPropertyTest(unittest.TestCase):
@@ -25,7 +25,7 @@ class PhysicalPropertyTest(unittest.TestCase):
:param expected_result: The expected result.
"""
np.testing.assert_allclose(
- p(vals).magnitude,
+ p(vals),
res.to(p.output_unit).magnitude,
rtol=self.EPS,
verbose=True,
@@ -33,84 +33,91 @@ class PhysicalPropertyTest(unittest.TestCase):
def test_temperature(self):
"""Test temperature property."""
- self.equality(pp.temperature, 293.15, Q_(20, "°C"))
- self.equality(pp.temperature, [293.15, 303.15], Q_([20, 30], "°C"))
+ self.equality(pp.temperature, 293.15, qty(20, "°C"))
+ self.equality(pp.temperature, [293.15, 303.15], qty([20, 30], "°C"))
def test_pressure(self):
"""Test pressure property."""
- self.equality(pp.pressure, 1e6, Q_(1, "MPa"))
- self.equality(pp.pressure, [1e6, 2e6], Q_([1, 2], "MPa"))
+ self.equality(pp.pressure, 1e6, qty(1, "MPa"))
+ self.equality(pp.pressure, [1e6, 2e6], qty([1, 2], "MPa"))
def test_velocity(self):
"""Test velocity property."""
- self.equality(pp.velocity.magnitude, [3, 4], Q_(5, "m/s"))
+ self.equality(pp.velocity.magnitude, [3, 4], qty(5, "m/s"))
self.equality(
- pp.velocity.magnitude, [[3, 4], [1, 0]], Q_([5, 1], "m/s")
+ pp.velocity.magnitude, [[3, 4], [1, 0]], qty([5, 1], "m/s")
)
- self.equality(pp.velocity.log_magnitude, np.sqrt([50, 50]), Q_(1, ""))
+ self.equality(pp.velocity.log_magnitude, np.sqrt([50, 50]), qty(1, ""))
self.equality(
pp.velocity.log_magnitude,
[np.sqrt([50, 50]), [10, 0]],
- Q_([1, 1], ""),
+ qty([1, 1], ""),
)
def test_displacement(self):
"""Test displacement property."""
- self.equality(pp.displacement.magnitude, [3e-3, 4e-3], Q_(5, "mm"))
+ self.equality(pp.displacement.magnitude, [3e-3, 4e-3], qty(5, "mm"))
u = np.array([[2, 3, 6], [1, 4, 8]]) * 1e-3
- self.equality(pp.displacement.magnitude, u, Q_([7.0, 9.0], "mm"))
+ self.equality(pp.displacement.magnitude, u, qty([7.0, 9.0], "mm"))
def test_strain(self):
"""Test strain property."""
eps = np.array([1, 3, 9, 1, 2, 2]) * 1e-2
- self.equality(pp.strain.magnitude, eps, Q_(10, "%"))
- self.equality(pp.strain.magnitude, [eps, eps], Q_([10, 10], "%"))
- self.equality(pp.strain.trace, eps, Q_(13, "%"))
- self.equality(pp.strain.trace, [eps, eps], Q_([13, 13], "%"))
+ self.equality(pp.strain.magnitude, eps, qty(10, "%"))
+ self.equality(pp.strain.magnitude, [eps, eps], qty([10, 10], "%"))
+ self.equality(pp.strain.trace, eps, qty(13, "%"))
+ self.equality(pp.strain.trace, [eps, eps], qty([13, 13], "%"))
def test_components(self):
"""Test strain components."""
eps = np.array([0, 1, 2, 3, 4, 5]) * 1e-2
u = np.array([0, 1, 2]) * 1e-3
for i in range(len(eps)):
- self.equality(pp.strain[i], eps, Q_(i, "%"))
- self.equality(pp.strain[i], [eps, eps], Q_([i, i], "%"))
+ self.equality(pp.strain[i], eps, qty(i, "%"))
+ self.equality(pp.strain[i], [eps, eps], qty([i, i], "%"))
for i in range(len(u)):
- self.equality(pp.displacement[i], u, Q_(i, "mm"))
- self.equality(pp.displacement[i], [u, u], Q_([i, i], "mm"))
+ self.equality(pp.displacement[i], u, qty(i, "mm"))
+ self.equality(pp.displacement[i], [u, u], qty([i, i], "mm"))
assert pp.strain[0].bilinear_cmap is True
def test_von_mises(self):
"""Test von_mises_stress property."""
- sig_3D = np.array([4, 1, 2, 1, 1, 1]) * 1e6
- self.equality(pp.von_mises_stress, sig_3D, Q_(4, "MPa"))
- self.equality(pp.von_mises_stress, [sig_3D, sig_3D], Q_([4, 4], "MPa"))
- sig_2D = np.array([4, 1, 2, 3**0.5]) * 1e6
- self.equality(pp.von_mises_stress, sig_2D, Q_(4, "MPa"))
+ sig_3D = np.array([3, 1, 1, 1, 1, 1]) * 1e6
+ self.equality(pp.von_mises_stress, sig_3D, qty(2, "MPa"))
+ self.equality(pp.von_mises_stress, [sig_3D, sig_3D], qty([2, 2], "MPa"))
+ sig_2D = np.array([2, 1, 1, 1]) * 1e6
+ self.equality(pp.von_mises_stress, sig_2D, qty(1, "MPa"))
def test_eff_pressure(self):
"""Test effective_pressure property."""
sig = np.array([-1, -2, -3, 1, 1, 1]) * 1e6
- self.equality(pp.effective_pressure, sig, Q_(2, "MPa"))
- self.equality(pp.effective_pressure, [sig, sig], Q_([2, 2], "MPa"))
+ self.equality(pp.effective_pressure, sig, qty(2, "MPa"))
+ self.equality(pp.effective_pressure, [sig, sig], qty([2, 2], "MPa"))
def test_qp_ratio(self):
"""Test qp_ratio property."""
+ sig = np.array([4, 4, 1, 1, 1, 1]) * 1e6
+ self.equality(pp.qp_ratio, sig, qty(-100, "percent"))
+ self.equality(pp.qp_ratio, [sig] * 2, qty([-100] * 2, "percent"))
+
+ def test_dilatancy(self):
+ """Test dilatancy property."""
sig = np.array([4, 1, 2, 1, 1, 1]) * 1e6
- self.equality(pp.qp_ratio, sig, Q_(-100 * 12 / 7, "percent"))
- self.equality(
- pp.qp_ratio, [sig] * 2, Q_([-100 * 12 / 7] * 2, "percent")
- )
+ self.assertRaises(TypeError, pp.dilatancy_alkan_eff, sig)
def test_simple(self):
"""Test call functionality."""
- self.assertEqual(pp.temperature(273.15), Q_(0, "°C"))
- self.assertEqual(pp.displacement[0]([1, 2, 3]), Q_(1, "m"))
- self.assertEqual(pp.displacement([1, 2, 3])[1], Q_(2, "m"))
+ self.assertEqual(pp.temperature(273.15, strip_unit=False), qty(0, "°C"))
+ self.assertEqual(
+ pp.displacement[0]([1, 2, 3], strip_unit=False), qty(1, "m")
+ )
+ self.assertEqual(
+ pp.displacement([1, 2, 3], strip_unit=False)[1], qty(2, "m")
+ )
def test_values(self):
"""Test values functionality."""
- self.assertEqual(pp.temperature.strip_units(273.15), 0.0)
+ self.assertEqual(pp.temperature(273.15), 0.0)
def test_units(self):
"""Test get_output_unit functionality."""
@@ -122,9 +129,13 @@ class PhysicalPropertyTest(unittest.TestCase):
"""Test mask functionality."""
self.assertTrue(pp.temperature.get_mask().is_mask())
- def test_find_property(self):
+ def test_get_preset(self):
"""Test find property function."""
- self.assertEqual(pp.find_property_preset("pressure"), pp.pressure)
+ self.assertEqual(pp.get_preset("pressure"), pp.pressure)
+ self.assertEqual(pp.get_preset("pore_pressure"), pp.pressure)
+ self.assertEqual(pp.get_preset("test"), Scalar("test"))
+ self.assertEqual(pp.get_preset("test", shape=(100, 3)), Vector("test"))
+ self.assertEqual(pp.get_preset("test", shape=(100, 6)), Matrix("test"))
def test_copy_ctor(self):
"""Test replace constructor."""