diff --git a/docs/conf.py b/docs/conf.py
index dc7d748bd74d37f351df59e8b30c057d164734d0..df9585dcb7ddf1c62073540a33f79481b225a7a3 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -117,6 +117,8 @@ autodoc_member_order = "bysource"
autodoc_preserve_defaults = True
autodoc_typehints = "description"
+autodoc_default_options = {"special-members": "__call__"}
+
### sphinx-gallery setup ###
# necessary when building the sphinx gallery
pyvista.BUILDING_GALLERY = True
diff --git a/docs/examples/howto_meshplotlib/plot_meshplotlib_2d.py b/docs/examples/howto_meshplotlib/plot_meshplotlib_2d.py
index 006e86e3c2e3f449b88db95cf15ff74b6b0be742..363064503006f47f1256d4fce93449625443648b 100644
--- a/docs/examples/howto_meshplotlib/plot_meshplotlib_2d.py
+++ b/docs/examples/howto_meshplotlib/plot_meshplotlib_2d.py
@@ -28,6 +28,14 @@ fig = mpl.plot(mesh, presets.material_id)
fig = mpl.plot(mesh, presets.temperature)
+# %%
+# We can also plot components of vector properties:
+
+fig = mpl.plot(mesh, presets.displacement[0])
+
+# %%
+fig = mpl.plot(mesh, presets.displacement[1])
+
# %%
# This example has hydraulically deactivated subdomains:
@@ -37,19 +45,8 @@ fig = mpl.plot(mesh, presets.pressure.get_mask())
# Let's plot the fluid velocity field.
fig = mpl.plot(mesh, presets.velocity)
-
# %%
# 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)
-mpl.setup.p_min = None
-mpl.setup.log_scaled = False
-
-# %%
-# We can also plot components of vector property:
-
-fig = mpl.plot(mesh, presets.displacement[0])
-
-# %%
-fig = mpl.plot(mesh, presets.displacement[1])
diff --git a/docs/examples/howto_meshplotlib/plot_solid_mechanics.py b/docs/examples/howto_meshplotlib/plot_solid_mechanics.py
index 22fd4b9044d42c0088634ec6698e2e809e1d67db..0cbfd245f7b0c3618fde696c8876d78832ef3e0c 100644
--- a/docs/examples/howto_meshplotlib/plot_solid_mechanics.py
+++ b/docs/examples/howto_meshplotlib/plot_solid_mechanics.py
@@ -42,23 +42,29 @@ fig = plot(mesh, presets.stress["xy"])
# %% [markdown]
# Principal stresses
# ------------------
-# Let's calculate the the principal stress components and overlay the direcetion
-# of the corresponding eigenvector. Note: the eigenvalues are sorted by
-# increasing order, i.e. eigenvalue[0] is the most negative / largest
+# Let's plot the the principal stress components and also overlay the direction
+# of the corresponding eigenvector in the plot. Note: the eigenvalues are sorted
+# by increasing order, i.e. eigenvalue[0] is the most negative / largest
# compressive principal stress.
# %%
eigvecs = presets.stress.eigenvectors
fig = plot(mesh, mesh_property=presets.stress.eigenvalues[0])
-plot_streamlines(ax=fig.axes[0], mesh=mesh, property=eigvecs[0], arrows=True)
+plot_streamlines(
+ ax=fig.axes[0], mesh=mesh, mesh_property=eigvecs[0], plot_type="lines"
+)
# %%
fig = plot(mesh, mesh_property=presets.stress.eigenvalues[1])
-plot_streamlines(ax=fig.axes[0], mesh=mesh, property=eigvecs[1], arrows=True)
+plot_streamlines(
+ ax=fig.axes[0], mesh=mesh, mesh_property=eigvecs[1], plot_type="lines"
+)
# %%
fig = plot(mesh, mesh_property=presets.stress.eigenvalues[2])
-plot_streamlines(ax=fig.axes[0], mesh=mesh, property=eigvecs[2], arrows=True)
+plot_streamlines(
+ ax=fig.axes[0], mesh=mesh, mesh_property=eigvecs[2], plot_type="lines"
+)
# %% [markdown]
# We can also plot the mean of the principal stress, i.e. the magnitude of the
diff --git a/docs/examples/howto_propertylib/plot_propertylib.py b/docs/examples/howto_propertylib/plot_propertylib.py
index 709895eb1581bd83ac854c3838cfce8c2a9cf0bf..1d536b115223a126d85a33b33be545b49549f2e2 100644
--- a/docs/examples/howto_propertylib/plot_propertylib.py
+++ b/docs/examples/howto_propertylib/plot_propertylib.py
@@ -11,35 +11,14 @@ mesh data. There are several predefined properties stored under the module
"""
# %%
-
-# sphinx_gallery_start_ignore
-
-import pandas as pd
-
from ogstools.meshplotlib import examples, plot
from ogstools.propertylib import Scalar, presets
-data = ["preset,data_name,data_unit,output_unit,output_name,type".split(",")]
-for preset_name in dir(presets):
- if isinstance(preset := presets.__dict__[preset_name], presets.Property):
- data += [
- [
- preset_name,
- preset.data_name,
- preset.data_unit,
- preset.output_unit,
- preset.output_name,
- preset.type_name,
- ]
- ]
-
-pd.DataFrame(data[1:], columns=data[0]).sort_values(
- ["data_name", "preset"]
-).set_index("preset")
-
-# sphinx_gallery_end_ignore
+presets.get_dataframe()
# %% [markdown]
+# Scalar, Vector and Matrix are inherit from the class Property with its
+# :meth:`~ogstools.propertylib.Property.__call__` operator.
# Calling a property converts the argument from data_unit to output_unit and
# applies a function if specified. In this case we convert from K to °C:
@@ -57,7 +36,7 @@ custom_temperature = Scalar(
custom_temperature(273.15)
# %% [markdown]
-# Or us existing presets as a template and replace some parameters:
+# Or use existing presets as a template and replace some parameters:
custom_temperature = presets.temperature.replace(output_unit="°F")
custom_temperature(273.15)
@@ -81,8 +60,8 @@ presets.strain["xx"]([0.01, 0.02, 0.03, 0.04, 0.05, 0.06])
presets.displacement.magnitude([0.03, 0.04])
# %% [markdown]
-# The main benefit of having specified how these properties should be
-# transformed is when making use of these in post processing. When plotting
+# We suggest specifying the properties and their transformations once.
+# These can be reused in different kind of post processing. When plotting
# with :py:mod:`ogstools.meshplotlib` we can use these presets to simplify the
# task of processing the data (e.g. calculate the von Mises stress):
diff --git a/ogstools/meshplotlib/core.py b/ogstools/meshplotlib/core.py
index a5e0444c16ded9378e569ac33cf0b8e3bb2761c6..87ea86b24bec855fe75b5843fef15fad11d710a2 100644
--- a/ogstools/meshplotlib/core.py
+++ b/ogstools/meshplotlib/core.py
@@ -21,7 +21,7 @@ from ogstools.propertylib.unit_registry import u_reg
from . import plot_features as pf
from . import setup
-from .levels import get_levels, median_exponent
+from .levels import compute_levels, median_exponent
from .utils import get_style_cycler
# TODO: define default data_name for regions in setup
@@ -236,7 +236,7 @@ def subplot(
if levels is None:
num_levels = min(setup.num_levels, len(np.unique(values)))
- levels = get_levels(p_min, p_max, num_levels)
+ levels = compute_levels(p_min, p_max, num_levels)
cmap, norm = get_cmap_norm(levels, mesh_property)
if (
@@ -372,7 +372,7 @@ def get_combined_levels(
and setup.p_max is None
):
return unique_vals[(p_min <= unique_vals) & (unique_vals <= p_max)]
- return get_levels(p_min, p_max, setup.num_levels)
+ return compute_levels(p_min, p_max, setup.num_levels)
def _plot_on_figure(
diff --git a/ogstools/meshplotlib/levels.py b/ogstools/meshplotlib/levels.py
index 55cfde5e589153c3d5a8423f6c8115cdccf33831..997048377ef9f59134303a0cef79816253295ec1 100644
--- a/ogstools/meshplotlib/levels.py
+++ b/ogstools/meshplotlib/levels.py
@@ -47,7 +47,7 @@ def adaptive_rounding(vals: np.ndarray, precision: int) -> np.ndarray:
return np.stack([np.round(v, 12 - median_exp) for v in vals])
-def get_levels(lower: float, upper: float, n_ticks: int) -> np.ndarray:
+def compute_levels(lower: float, upper: float, n_ticks: int) -> np.ndarray:
"""
Return an array in the interval [lower, upper] with terminating decimals.
@@ -56,10 +56,10 @@ def get_levels(lower: float, upper: float, n_ticks: int) -> np.ndarray:
"""
if lower == upper:
return np.asarray([lower, nextafter(lower, np.inf)])
- levels = nice_range(lower, upper, n_ticks)
+ result = nice_range(lower, upper, n_ticks)
return np.unique(
adaptive_rounding(
- np.append(np.append(lower, levels), upper), precision=3
+ np.append(np.append(lower, result), upper), precision=3
)
)
diff --git a/ogstools/meshplotlib/plot_features.py b/ogstools/meshplotlib/plot_features.py
index 41a8d6b9ee21acaf74d0188fce040e7e03e6a72c..7a9a3e1d4f16dda90da0352162488921f4ada4ad 100644
--- a/ogstools/meshplotlib/plot_features.py
+++ b/ogstools/meshplotlib/plot_features.py
@@ -1,6 +1,6 @@
"""Specialized plot features."""
-from typing import Callable, Optional
+from typing import Callable, Literal, Optional
import matplotlib.pyplot as plt
import numpy as np
@@ -81,14 +81,23 @@ def plot_element_edges(ax: plt.Axes, surf: pv.DataSet, projection: int) -> None:
def plot_streamlines(
ax: plt.Axes,
mesh: pv.DataSet,
- property: Vector,
+ mesh_property: Vector,
projection: Optional[int] = None,
- arrows: bool = False,
+ plot_type: Literal["streamlines", "arrows", "lines"] = "streamlines",
) -> None:
- """Plot vector streamlines on a matplotlib axis."""
+ """
+ Plot the vector streamlines or arrows on a matplotlib axis.
+
+ :param ax: Matplotlib axis to plot onto
+ :param mesh: Mesh containing the vector property
+ :param mesh_property: Vector property to visualize
+ :param projection: Index of flat dimension (e.g. 2 for z axis),
+ gets automatically determined if not given
+ :param plot_type: Whether to plot streamlines, arrows or lines.
+ """
if (n_pts := setup.num_streamline_interp_pts) is None:
return
- if arrows:
+ if plot_type != "streamlines":
n_pts = 50
if projection is None:
mean_normal = np.abs(
@@ -103,47 +112,37 @@ def plot_streamlines(
_mesh = mesh.copy()
for key in _mesh.point_data:
- if key not in [property.data_name, property.mask]:
+ if key not in [mesh_property.data_name, mesh_property.mask]:
del _mesh.point_data[key]
grid = pv.RectilinearGrid(
[x, y, z][x_id], [x, y, z][y_id], [x, y, z][projection]
)
grid = grid.sample(_mesh, pass_cell_data=False)
- values = property(grid.point_data[property.data_name])
+ values = mesh_property(grid.point_data[mesh_property.data_name])
values[np.argwhere(grid["vtkValidPointMask"] == 0), :] = np.nan
if np.shape(values)[-1] == 3:
values = np.delete(values, projection, 1)
val = np.reshape(values, (n_pts, n_pts, 2))
- if property.mask in grid.point_data:
- mask = np.reshape(grid.point_data[property.mask], (n_pts, n_pts))
+ if mesh_property.mask in grid.point_data:
+ mask = np.reshape(grid.point_data[mesh_property.mask], (n_pts, n_pts))
val[mask == 0, :] = 0
val_norm = np.linalg.norm(np.nan_to_num(val), axis=-1)
lw = 2.5 * val_norm / max(1e-16, np.max(val_norm))
lw *= setup.rcParams_scaled["lines.linewidth"]
x_g, y_g = setup.length(np.meshgrid(x, y))
- if arrows:
- ax.quiver(
- x_g,
- y_g,
- val[..., 0],
- val[..., 1],
- scale=1 / 0.03,
- headlength=0,
- headaxislength=0,
- headwidth=1,
- pivot="mid",
- )
+ 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)
else:
- ax.streamplot(
- x_g,
- y_g,
- val[..., 0],
- val[..., 1],
- color="k",
- linewidth=lw,
- density=1.5,
+ line_args = (
+ dict( # noqa: C408
+ headlength=0, headaxislength=0, headwidth=1, pivot="mid"
+ )
+ if plot_type == "lines"
+ else {}
)
+ ax.quiver(*plot_args, **line_args, scale=1 / 0.03)
def plot_on_top(
diff --git a/ogstools/propertylib/matrix.py b/ogstools/propertylib/matrix.py
index d8ad499a181e7dfcc6c671f2e8c76e61e2767489..e37c88117119b9fc5730d910943239323244fa44 100644
--- a/ogstools/propertylib/matrix.py
+++ b/ogstools/propertylib/matrix.py
@@ -78,32 +78,32 @@ class Matrix(Property):
)
@property
- def I1(self) -> Scalar: # pylint: disable=invalid-name
+ def invariant_1(self) -> Scalar:
"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)),
+ func=lambda x: tensor_math.invariant_1(self.func(x)),
)
@property
- def I2(self) -> Scalar: # pylint: disable=invalid-name
+ def invariant_2(self) -> Scalar:
"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)),
+ func=lambda x: tensor_math.invariant_2(self.func(x)),
process_with_units=True,
)
@property
- def I3(self) -> Scalar: # pylint: disable=invalid-name
+ def invariant_3(self) -> Scalar:
"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)),
+ func=lambda x: tensor_math.invariant_3(self.func(x)),
)
@property
@@ -134,30 +134,30 @@ class Matrix(Property):
)
@property
- def J1(self) -> Scalar: # pylint: disable=invalid-name
+ def deviator_invariant_1(self) -> Scalar:
"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)),
+ func=lambda x: tensor_math.deviator_invariant_1(self.func(x)),
)
@property
- def J2(self) -> Scalar: # pylint: disable=invalid-name
+ def deviator_invariant_2(self) -> Scalar:
"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)),
+ func=lambda x: tensor_math.deviator_invariant_2(self.func(x)),
)
@property
- def J3(self) -> Scalar: # pylint: disable=invalid-name
+ def deviator_invariant_3(self) -> Scalar:
"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)),
+ func=lambda x: tensor_math.deviator_invariant_3(self.func(x)),
)
@property
diff --git a/ogstools/propertylib/mesh_dependent.py b/ogstools/propertylib/mesh_dependent.py
index c26938ba0956036d77e17c3ebb48310efe2fb778..6f9682442eb8792cfe03cda95661e0256afbbb14 100644
--- a/ogstools/propertylib/mesh_dependent.py
+++ b/ogstools/propertylib/mesh_dependent.py
@@ -1,6 +1,5 @@
"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
@@ -17,8 +16,9 @@ ValType = Union[PlainQuantity, np.ndarray]
def depth(mesh: pv.UnstructuredGrid, use_coords: bool = False) -> np.ndarray:
"""Return the depth values of the mesh.
- For 2D, the last non-flat dimension is used as the vertical axis, for 3D,
- the z-axes is used.
+ For 2D, the last axis of the plane wherein the mesh is lying is used as the
+ vertical axis (i.e. y if the mesh is in the xy-plane, z if it is in the
+ xz-plane), for 3D, the z-axes is used.
If `use_coords` is True, returns the negative coordinate value of the
vertical axis. Otherwise, the vertical distance to the top facing edges /
surfaces are returned.
@@ -26,18 +26,18 @@ def depth(mesh: pv.UnstructuredGrid, use_coords: bool = False) -> np.ndarray:
if mesh.volume > 0:
# prevents inner edge (from holes) to be detected as a top boundary
edges = mesh.extract_surface().connectivity("point_seed", point_ids=[0])
- top_id = 2
+ vertical_dim = 2
if use_coords:
- return -mesh.points[:, top_id]
- point_upwards = edges.cell_normals[..., top_id] > 0
+ return -mesh.points[:, vertical_dim]
+ point_upwards = edges.cell_normals[..., vertical_dim] > 0
top_cells = point_upwards
else:
mean_normal = np.abs(
np.mean(mesh.extract_surface().cell_normals, axis=0)
)
- top_id = np.delete([0, 1, 2], int(np.argmax(mean_normal)))[-1]
+ vertical_dim = np.delete([0, 1, 2], int(np.argmax(mean_normal)))[-1]
if use_coords:
- return -mesh.points[:, top_id]
+ return -mesh.points[:, vertical_dim]
# prevents inner edge (from holes) to be detected as a top boundary
edges = mesh.extract_feature_edges().connectivity(
"point_seed", point_ids=[0]
@@ -49,16 +49,23 @@ 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[..., top_id] > adj_cell_centers[..., top_id]
+ are_above = (
+ edge_centers[..., vertical_dim]
+ > adj_cell_centers[..., vertical_dim]
+ )
are_non_vertical = np.asarray(edge_horizontal_extent) > 1e-12
top_cells = are_above & are_non_vertical
top = edges.extract_cells(top_cells)
eucl_distance_projected_top_points = np.sum(
- np.abs(np.delete(mesh.points[:, None] - top.points, top_id, axis=-1)),
+ np.abs(
+ np.delete(mesh.points[:, None] - top.points, vertical_dim, axis=-1)
+ ),
axis=-1,
)
matching_top = np.argmin(eucl_distance_projected_top_points, axis=-1)
- return np.abs(mesh.points[..., top_id] - top.points[matching_top, top_id])
+ return np.abs(
+ mesh.points[..., vertical_dim] - top.points[matching_top, vertical_dim]
+ )
def p_fluid(mesh: pv.UnstructuredGrid) -> PlainQuantity:
@@ -74,11 +81,11 @@ def p_fluid(mesh: pv.UnstructuredGrid) -> PlainQuantity:
where `h` is the depth below surface. If "depth" is not given in the mesh,
it is calculated via :py:func:`ogstools.propertylib.mesh_dependent.depth`.
"""
- qty = u_reg.Quantity
+ Qty = u_reg.Quantity
if "pressure" in mesh.point_data:
- return qty(mesh["pressure"], "Pa")
+ 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")
+ return Qty(1000, "kg/m^3") * Qty(9.81, "m/s^2") * Qty(_depth, "m")
def fluid_pressure_criterion(
@@ -97,8 +104,8 @@ def fluid_pressure_criterion(
:py:func:`ogstools.propertylib.mesh_dependent.p_fluid`.
"""
- qty = u_reg.Quantity
- sigma = -qty(mesh[mesh_property.data_name], mesh_property.data_unit)
+ Qty = u_reg.Quantity
+ sigma = -Qty(mesh[mesh_property.data_name], mesh_property.data_unit)
return p_fluid(mesh) - eigenvalues(sigma)[..., 0]
@@ -115,12 +122,21 @@ def dilatancy_critescu(
F_{dil} = \\frac{\\tau_{oct}}{\\sigma_0} - a \\left( \\frac{\\sigma_m}{\\sigma_0} \\right)^2 - b \\frac{\\sigma_m}{\\sigma_0}
+ for total stresses and defined as:
+
+ .. math::
+
+ F'_{dil} = \\frac{\\tau_{oct}}{\\sigma_0} - a \\left( \\frac{\\sigma'_m}{\\sigma_0} \\right)^2 - b \\frac{\\sigma'_m}{\\sigma_0}
+
+ for effective stresses
+ (uses :func:`~ogstools.propertylib.mesh_dependent.p_fluid`).
+
<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")
+ 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)
@@ -130,17 +146,6 @@ def dilatancy_critescu(
)
-dilatancy_critescu_eff = partial(dilatancy_critescu, effective=True)
-"""Return the dilatancy criterion defined as:
-
-.. math::
-
- F'_{dil} = \\frac{\\tau_{oct}}{\\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>
-"""
-
-
def dilatancy_alkan(
mesh: pv.UnstructuredGrid,
mesh_property: Property,
@@ -153,25 +158,23 @@ def dilatancy_alkan(
F_{dil} = \\tau_{oct} - \\tau_{max} \\cdot b \\frac{\\sigma'_m}{\\sigma_0 + b \\cdot \\sigma'_m}
+ for total stresses and defined as:
+
+ .. math::
+
+ F_{dil} = \\tau_{oct} - \\tau_{max} \\cdot b \\frac{\\sigma'_m}{\\sigma_0 + b \\cdot \\sigma'_m}
+
+ for effective stresses
+ (uses :func:`~ogstools.propertylib.mesh_dependent.p_fluid`).
+
<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")
+ 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)
-"""Return the dilatancy criterion defined as:
-
-.. math::
-
- F_{dil} = \\tau_{oct} - \\tau_{max} \\cdot b \\frac{\\sigma'_m}{\\sigma_0 + b \\cdot \\sigma'_m}
-
-<https://www.sciencedirect.com/science/article/pii/S1365160906000979>
-"""
+ return tau - tau_max * (b * sigma_m / (Qty(1, "MPa") + b * sigma_m))
diff --git a/ogstools/propertylib/presets.py b/ogstools/propertylib/presets.py
index 139faac875debf1d1ee142f1cd3480c04bdd578a..f53b6c17b2813509e2f29c6303336b980f5a9270 100644
--- a/ogstools/propertylib/presets.py
+++ b/ogstools/propertylib/presets.py
@@ -4,8 +4,11 @@
".. seealso:: :py:mod:`ogstools.propertylib.tensor_math`"
"""
+from functools import partial
from typing import Optional
+import pandas as pd
+
from . import mesh_dependent, tensor_math
from .custom_colormaps import integrity_cmap, temperature_cmap
from .matrix import Matrix
@@ -99,17 +102,11 @@ dilatancy_critescu = Scalar(
cmap=integrity_cmap,
bilinear_cmap=True,
)
-dilatancy_critescu_eff = Scalar(
- data_name="sigma",
- data_unit="Pa",
- output_unit="",
+dilatancy_critescu_eff = dilatancy_critescu.replace(
output_name="effective_dilatancy_criterion",
- mask=M_MASK,
- func=mesh_dependent.dilatancy_critescu_eff,
- mesh_dependent=True,
- cmap=integrity_cmap,
- bilinear_cmap=True,
+ func=partial(mesh_dependent.dilatancy_critescu, effective=True),
)
+
dilatancy_alkan = Scalar(
data_name="sigma",
data_unit="Pa",
@@ -121,17 +118,11 @@ dilatancy_alkan = Scalar(
cmap=integrity_cmap,
bilinear_cmap=True,
)
-dilatancy_alkan_eff = Scalar(
- data_name="sigma",
- data_unit="Pa",
- output_unit="MPa",
+dilatancy_alkan_eff = dilatancy_alkan.replace(
output_name="effective_dilatancy_criterion",
- mask=M_MASK,
- func=mesh_dependent.dilatancy_alkan_eff,
- mesh_dependent=True,
- cmap=integrity_cmap,
- bilinear_cmap=True,
+ func=partial(mesh_dependent.dilatancy_alkan, effective=True),
)
+
fluid_pressure_crit = Scalar(
data_name="sigma",
data_unit="Pa",
@@ -167,3 +158,27 @@ def get_preset(property_name: str, shape: Optional[tuple] = None) -> Property:
if shape[1] in [2, 3]:
return Vector(property_name)
return Matrix(property_name)
+
+
+def get_dataframe() -> pd.DataFrame:
+ data = [
+ "preset,data_name,data_unit,output_unit,output_name,type".split(",")
+ ]
+ for preset_name, preset_value in globals().items():
+ if isinstance(preset := preset_value, Property):
+ data += [
+ [
+ preset_name,
+ preset.data_name,
+ preset.data_unit,
+ preset.output_unit,
+ preset.output_name,
+ preset.type_name,
+ ]
+ ]
+
+ return (
+ pd.DataFrame(data[1:], columns=data[0])
+ .sort_values(["data_name", "preset"])
+ .set_index("preset")
+ )
diff --git a/ogstools/propertylib/property.py b/ogstools/propertylib/property.py
index fe83d01dfc42c5ef893aad41b29f3fb77b5ea5b7..c61a6c8120e8fac14d5013ae35e789c88e6cc6e5 100644
--- a/ogstools/propertylib/property.py
+++ b/ogstools/propertylib/property.py
@@ -33,7 +33,8 @@ class Property:
mask: str = ""
"""The name of the mask data in the dataset."""
func: Callable = identity
- """The function to be applied on the data."""
+ """The function to be applied on the data.
+ .. seealso:: :meth:`~ogstools.propertylib.Property.__call__`"""
mesh_dependent: bool = False
"""If the function to be applied is dependent on the mesh itself"""
process_with_units: bool = False
@@ -90,27 +91,30 @@ class Property:
"""
Return the transformed data values.
- Apply property function, convert from data_unit to output_unit and
- return the values, optionally with the unit.
+ Converts the data from data_unit to output_unit and apply a function
+ the transformation function of this property. The result is returned by
+ default without units. if `strip_unit` is False, a quantity is returned.
- :param vals: The input data values.
-
- :returns: The transformed data values.
+ 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
+ data is passed to the function with units (i.e. as a pint quantity) or
+ without.
"""
- qty, _du, _ou = u_reg.Quantity, self.data_unit, self.output_unit
+ Qty, d_u, o_u = 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)
+ 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):
- result = qty(self.func(qty(self.get_data(data), _du)), _ou)
+ result = Qty(self.func(Qty(self._get_data(data), d_u)), o_u)
elif self.process_with_units:
- result = qty(self.func(qty(data, _du)), _ou)
+ result = Qty(self.func(Qty(data, d_u)), o_u)
else:
- result = qty(qty(self.func(np.asarray(data)), _du), _ou)
+ result = Qty(Qty(self.func(np.asarray(data)), d_u), o_u)
return result.magnitude if strip_unit else result
def get_output_unit(self) -> str:
@@ -142,13 +146,14 @@ class Property:
return self
def mask_used(self, mesh: pv.UnstructuredGrid) -> bool:
+ "Check whether the mesh contains the mask of this property."
return (
not self.is_mask()
and self.mask in mesh.cell_data
and (len(mesh.cell_data[self.mask]) != 0)
)
- def get_data(
+ 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."""
diff --git a/ogstools/propertylib/tensor_math.py b/ogstools/propertylib/tensor_math.py
index 78a4d8cea05b20b7492e3c3b9fe2d8d8ec857d9c..599a5a78a991fe6fb19ec86f17271d97e702e589 100644
--- a/ogstools/propertylib/tensor_math.py
+++ b/ogstools/propertylib/tensor_math.py
@@ -103,7 +103,7 @@ def frobenius_norm(val: ValType) -> ValType:
return np.linalg.norm(sym_tensor_to_mat(val), axis=(-2, -1))
-def I1(vals: ValType) -> ValType: # pylint: disable=invalid-name
+def invariant_1(vals: ValType) -> ValType:
"""Return the first invariant.
:math:`I1 = tr(\\mathbf{\\sigma})`
@@ -111,7 +111,7 @@ def I1(vals: ValType) -> ValType: # pylint: disable=invalid-name
return trace(vals)
-def I2(vals: ValType) -> ValType: # pylint: disable=invalid-name
+def invariant_2(vals: ValType) -> ValType:
"""Return the second invariant.
:math:`I2 = \\frac{1}{2} \\left[(tr(\\mathbf{\\sigma}))^2 - tr(\\mathbf{\\sigma}^2) \\right]`
@@ -119,7 +119,7 @@ def I2(vals: ValType) -> ValType: # pylint: disable=invalid-name
return 0.5 * (trace(vals) ** 2 - trace(vals**2))
-def I3(vals: ValType) -> ValType: # pylint: disable=invalid-name
+def invariant_3(vals: ValType) -> ValType:
"""Return the third invariant.
:math:`I3 = det(\\mathbf{\\sigma})`
@@ -133,7 +133,7 @@ def mean(vals: ValType) -> ValType:
:math:`\\pi = \\frac{1}{3} I1`
"""
- return (1.0 / 3.0) * I1(vals)
+ return (1.0 / 3.0) * invariant_1(vals)
def effective_pressure(vals: ValType) -> ValType:
@@ -167,7 +167,7 @@ def deviator(vals: ValType) -> ValType:
return vals - hydrostatic_component(vals)
-def J1(vals: ValType) -> ValType: # pylint: disable=invalid-name
+def deviator_invariant_1(vals: ValType) -> ValType:
"""Return the first invariant of the deviator.
:math:`J1 = 0`
@@ -175,7 +175,7 @@ def J1(vals: ValType) -> ValType: # pylint: disable=invalid-name
return trace(deviator(vals))
-def J2(vals: ValType) -> ValType: # pylint: disable=invalid-name
+def deviator_invariant_2(vals: ValType) -> ValType:
"""Return the second invariant of the deviator.
:math:`J2 = \\frac{1}{2} tr(\\mathbf{s}^2)`
@@ -183,7 +183,7 @@ def J2(vals: ValType) -> ValType: # pylint: disable=invalid-name
return 0.5 * trace(deviator(vals) ** 2)
-def J3(vals: ValType) -> ValType: # pylint: disable=invalid-name
+def deviator_invariant_3(vals: ValType) -> ValType:
"""Return the third invariant of the deviator.
:math:`J3 = \\frac{1}{3} tr(\\mathbf{s}^3)`
@@ -196,7 +196,7 @@ def octahedral_shear(vals: ValType) -> ValType:
:math:`\\tau_{oct} = \\sqrt{\\frac{2}{3} J2}`
"""
- return np.sqrt((2.0 / 3.0) * J2(vals))
+ return np.sqrt((2.0 / 3.0) * deviator_invariant_2(vals))
def von_mises(vals: ValType) -> ValType:
@@ -204,7 +204,7 @@ def von_mises(vals: ValType) -> ValType:
:math:`\\sigma_{Mises} = \\sqrt{3 J2}`
"""
- return np.sqrt(3.0 * J2(vals))
+ return np.sqrt(3.0 * deviator_invariant_2(vals))
def qp_ratio(vals: ValType) -> ValType:
diff --git a/ogstools/propertylib/vector.py b/ogstools/propertylib/vector.py
index a5a2b73dc2c2fc327335955f36a54a4a2fba4ec3..25abd336fa89bfe8f211626a6b98a9f9e39cc0f6 100644
--- a/ogstools/propertylib/vector.py
+++ b/ogstools/propertylib/vector.py
@@ -12,7 +12,7 @@ ValType = Union[PlainQuantity, np.ndarray]
def vector_norm(vals: ValType) -> ValType:
- "Return the magnitudes."
+ ":returns: The norm of the vector."
if isinstance(vals, PlainQuantity):
unit = vals.units
vals = vals.magnitude
@@ -61,9 +61,7 @@ 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.
- """
+ ":returns: A vector property as a component of the vectorlist property."
return Vector.from_property(
self,
output_name=self.output_name + f"_{index}",
diff --git a/pyproject.toml b/pyproject.toml
index 292e85a3735f5e3cc4550fa1d86dac51cf2ea9a2..d4253813d11c0f99e7cfbaae4734cd50dc9da821 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -49,7 +49,7 @@ feflow2ogs = 'ogstools.feflowlib._cli:cli'
[project.optional-dependencies]
dev = ["pre-commit>=2.20", "build", "black"]
-test = ["pytest", "coverage", "gmsh", "ogs6py"]
+test = ["pytest", "coverage", "gmsh", "ogs6py", "parameterized"]
docs = [
"sphinx",
"sphinx-argparse",
diff --git a/tests/test_mechanics.py b/tests/test_mechanics.py
index 611eaba7fd242d8f79c29da4e260650753e63b30..d32072841d0be1f8e7ec45aa43f144f66b8bd5b6 100644
--- a/tests/test_mechanics.py
+++ b/tests/test_mechanics.py
@@ -3,127 +3,127 @@
import unittest
import numpy as np
+from parameterized import parameterized
from ogstools.propertylib import tensor_math
from ogstools.propertylib.tensor_math import sym_tensor_to_mat
+def assert_allclose(vals1: np.ndarray, vals2: np.ndarray, rtol=1e-7, atol=1e-9):
+ """Assert the equality of two arrays."""
+ np.testing.assert_allclose(vals1, vals2, verbose=True, rtol=rtol, atol=atol)
+
+
class MechanicsTest(unittest.TestCase):
"""Test case for physical properties."""
rng = np.random.default_rng()
- N_SAMPLES = 100000
-
- def equality(
- self, vals1: np.ndarray, vals2: np.ndarray, rtol=1e-7, atol=1e-9
- ):
- """Assert the equality of two arrays."""
- np.testing.assert_allclose(
- vals1, vals2, verbose=True, rtol=rtol, atol=atol
- )
- def test_frobenius_norm(self):
+ def generate_random_sig(self, symten_len: int):
+ """Generate random stress matrix."""
+ return self.rng.random((100000, symten_len)) * 1e6
+
+ @parameterized.expand([(4,), (6,)])
+ def test_frobenius_norm(self, symten_len: int):
"""Test Frobenius norm."""
- for symten_len in [4, 6]:
- sig = self.rng.random((self.N_SAMPLES, symten_len)) * 1e6
- sig_mat = sym_tensor_to_mat(sig)
- frob2 = np.sqrt(
- tensor_math.trace(
- np.diagonal(
- np.transpose(sig_mat, (0, 2, 1)) @ sig_mat, 0, 2
- )
- )
+ sig = self.generate_random_sig(symten_len)
+ sig_mat = sym_tensor_to_mat(sig)
+ frob2 = np.sqrt(
+ tensor_math.trace(
+ np.diagonal(np.transpose(sig_mat, (0, 2, 1)) @ sig_mat, 0, 2)
)
- self.equality(tensor_math.frobenius_norm(sig), frob2)
+ )
+ assert_allclose(tensor_math.frobenius_norm(sig), frob2)
- def test_I1(self):
+ @parameterized.expand([(4,), (6,)])
+ def test_invariant_1(self, symten_len: int):
"""Test first invariant."""
- for symten_len in [4, 6]:
- sig = self.rng.random((self.N_SAMPLES, symten_len)) * 1e6
- self.equality(
- tensor_math.I1(sig),
- tensor_math.trace(tensor_math.eigenvalues(sig)),
- )
+ sig = self.generate_random_sig(symten_len)
+ assert_allclose(
+ tensor_math.invariant_1(sig),
+ tensor_math.trace(tensor_math.eigenvalues(sig)),
+ )
- def test_I2(self):
+ @parameterized.expand([(4,), (6,)])
+ def test_invariant_2(self, symten_len: int):
"""Test second invariant."""
- for symten_len in [4, 6]:
- sig = self.rng.random((self.N_SAMPLES, symten_len)) * 1e6
- eig_vals = tensor_math.eigenvalues(sig)
- self.equality(
- tensor_math.I2(sig),
- np.sum(eig_vals * np.roll(eig_vals, 1, axis=-1), axis=-1),
- )
+ sig = self.generate_random_sig(symten_len)
+ eig_vals = tensor_math.eigenvalues(sig)
+ assert_allclose(
+ tensor_math.invariant_2(sig),
+ np.sum(eig_vals * np.roll(eig_vals, 1, axis=-1), axis=-1),
+ )
- def test_I3(self):
+ @parameterized.expand([(4,), (6,)])
+ def test_invariant_3(self, symten_len: int):
"""Test third invariant."""
- for symten_len in [4, 6]:
- sig = self.rng.random((self.N_SAMPLES, symten_len)) * 1e6
- eig_vals = tensor_math.eigenvalues(sig)
- self.equality(tensor_math.I3(sig), np.prod(eig_vals, axis=-1))
+ sig = self.generate_random_sig(symten_len)
+ eig_vals = tensor_math.eigenvalues(sig)
+ assert_allclose(
+ tensor_math.invariant_3(sig), np.prod(eig_vals, axis=-1)
+ )
- def test_J1(self):
+ @parameterized.expand([(4,), (6,)])
+ def test_deviator_invariant_1(self, symten_len: int):
"""Test first deviator invariant."""
- for symten_len in [4, 6]:
- sig = self.rng.random((self.N_SAMPLES, symten_len)) * 1e6
- j1 = tensor_math.J1(sig)
- self.equality(j1, np.zeros(j1.shape))
+ sig = self.generate_random_sig(symten_len)
+ j1 = tensor_math.deviator_invariant_1(sig)
+ assert_allclose(j1, np.zeros(j1.shape))
- def test_J2(self):
+ @parameterized.expand([(4,), (6,)])
+ def test_deviator_invariant_2(self, symten_len: int):
"""Test second deviator invariant."""
- for symten_len in [4, 6]:
- sig = self.rng.random((self.N_SAMPLES, symten_len)) * 1e6
- self.equality(
- tensor_math.J2(sig),
- 0.5
- * (
- tensor_math.trace(sig**2)
- - (1.0 / 3.0) * tensor_math.trace(sig) ** 2
- ),
- )
+ sig = self.generate_random_sig(symten_len)
+ assert_allclose(
+ tensor_math.deviator_invariant_2(sig),
+ 0.5
+ * (
+ tensor_math.trace(sig**2)
+ - (1.0 / 3.0) * tensor_math.trace(sig) ** 2
+ ),
+ )
- def test_J3(self):
+ @parameterized.expand([(4,), (6,)])
+ def test_deviator_invariant_3(self, symten_len: int):
"""Test third deviator invariant."""
- for symten_len in [4, 6]:
- sig = self.rng.random((self.N_SAMPLES, symten_len)) * 1e6
- # not exactly sure why, but for stresses where the below condition
- # is very close to zero, the error of this test shoots up.
- # Probably some floating point precision issue.
- mask = (
- np.abs(
- tensor_math.trace(sig - tensor_math.mean(sig)[..., None])
- )
- > 1e-9
- )
- sig = sig[mask]
- self.equality(
- tensor_math.J3(sig),
- (1.0 / 3.0)
- * (
- tensor_math.trace(sig**3)
- - tensor_math.trace(sig**2) * tensor_math.trace(sig)
- + (2.0 / 9.0) * tensor_math.trace(sig) ** 3.0
- ),
- )
+ sig = self.generate_random_sig(symten_len)
+ # not exactly sure why, but for stresses where the below condition
+ # is very close to zero, the error of this test shoots up.
+ # Probably some floating point precision issue.
+ mask = (
+ np.abs(tensor_math.trace(sig - tensor_math.mean(sig)[..., None]))
+ > 1e-9
+ )
+ sig = sig[mask]
+ assert_allclose(
+ tensor_math.deviator_invariant_3(sig),
+ (1.0 / 3.0)
+ * (
+ tensor_math.trace(sig**3)
+ - tensor_math.trace(sig**2) * tensor_math.trace(sig)
+ + (2.0 / 9.0) * tensor_math.trace(sig) ** 3.0
+ ),
+ )
- def test_von_mises(self):
+ @parameterized.expand([(4,), (6,)])
+ def test_von_mises(self, symten_len: int):
"""Test von Mises invariant."""
- for symten_len in [4, 6]:
- sig = self.rng.random((self.N_SAMPLES, symten_len)) * 1e6
- ev = tensor_math.eigenvalues(sig)
- self.equality(
- tensor_math.von_mises(sig),
- np.sqrt(0.5 * np.sum(np.square(ev - np.roll(ev, 1, -1)), -1)),
- )
+ sig = self.generate_random_sig(symten_len)
+ ev = tensor_math.eigenvalues(sig)
+ assert_allclose(
+ tensor_math.von_mises(sig),
+ np.sqrt(0.5 * np.sum(np.square(ev - np.roll(ev, 1, -1)), -1)),
+ )
- def test_octahedral_shear_stress(self):
+ @parameterized.expand([(4,), (6,)])
+ def test_octahedral_shear_stress(self, symten_len: int):
"""Test octahedral shear stress invariant."""
- for symten_len in [4, 6]:
- sig = self.rng.random((self.N_SAMPLES, symten_len)) * 1e6
- self.equality(
- tensor_math.octahedral_shear(sig),
- (1.0 / 3.0)
- * np.sqrt(
- 2 * tensor_math.I1(sig) ** 2 - 6 * tensor_math.I2(sig)
- ),
- )
+ sig = self.generate_random_sig(symten_len)
+ assert_allclose(
+ tensor_math.octahedral_shear(sig),
+ (1.0 / 3.0)
+ * np.sqrt(
+ 2 * tensor_math.invariant_1(sig) ** 2
+ - 6 * tensor_math.invariant_2(sig)
+ ),
+ )
diff --git a/tests/test_meshplotlib.py b/tests/test_meshplotlib.py
index 6e7911ef83d9bfbe9f63e9363b16e288c472dbc7..8df8486d581abf2a5fca1d49134637137e843d35 100644
--- a/tests/test_meshplotlib.py
+++ b/tests/test_meshplotlib.py
@@ -17,12 +17,12 @@ from ogstools.meshplotlib import (
)
from ogstools.meshplotlib.animation import animate, save_animation
from ogstools.meshplotlib.core import get_ticklabels
-from ogstools.meshplotlib.levels import get_levels
+from ogstools.meshplotlib.levels import compute_levels
from ogstools.meshplotlib.plot_features import plot_on_top
from ogstools.meshplotlib.utils import justified_labels
from ogstools.propertylib import Scalar, presets
-equality = partial(
+assert_allclose = partial(
np.testing.assert_allclose, rtol=1e-7, atol=1e-100, verbose=True
)
@@ -40,18 +40,26 @@ class MeshplotlibTest(unittest.TestCase):
def test_levels(self):
"""Test levels calculation property."""
- equality(get_levels(0.5, 10.1, 10), [0.5, *range(1, 11), 10.1])
- equality(get_levels(293, 350, 10), [293, *range(295, 355, 5)])
- equality(get_levels(1e-3, 1.2, 5), [1e-3, *np.arange(0.2, 1.4, 0.2)])
- equality(get_levels(1e5, 9e6, 20), [1e5, *np.arange(5e5, 9.5e6, 5e5)])
- equality(get_levels(1, 40, 20), [1, *range(2, 42, 2)])
- equality(get_levels(0.0, 0.0, 10), [0.0, 0.0])
- equality(get_levels(1e9, 1e9, 10), [1e9, 1e9])
+ assert_allclose(
+ compute_levels(0.5, 10.1, 10), [0.5, *range(1, 11), 10.1]
+ )
+ assert_allclose(
+ compute_levels(293, 350, 10), [293, *range(295, 355, 5)]
+ )
+ assert_allclose(
+ compute_levels(1e-3, 1.2, 5), [1e-3, *np.arange(0.2, 1.4, 0.2)]
+ )
+ assert_allclose(
+ compute_levels(1e5, 9e6, 20), [1e5, *np.arange(5e5, 9.5e6, 5e5)]
+ )
+ assert_allclose(compute_levels(1, 40, 20), [1, *range(2, 42, 2)])
+ assert_allclose(compute_levels(0.0, 0.0, 10), [0.0, 0.0])
+ assert_allclose(compute_levels(1e9, 1e9, 10), [1e9, 1e9])
def test_ticklabels(self):
def compare(lower, upper, precision, ref_labels, ref_offset=None):
labels, offset = get_ticklabels(
- np.asarray(get_levels(lower, upper, n_ticks=precision))
+ np.asarray(compute_levels(lower, upper, n_ticks=precision))
)
self.assertTrue(np.all(labels == ref_labels))
self.assertEqual(offset, ref_offset)
diff --git a/tests/test_propertylib.py b/tests/test_propertylib.py
index 2a203cf732b5185e52c79c117e1b26f46a41cfbc..373c39af3601a5db6895fe1484a1a7fb7d01ec44 100644
--- a/tests/test_propertylib.py
+++ b/tests/test_propertylib.py
@@ -13,7 +13,7 @@ from ogstools.propertylib.mesh_dependent import depth
from ogstools.propertylib.property import Scalar, u_reg
from ogstools.propertylib.vector import Vector
-qty = u_reg.Quantity
+Qty = u_reg.Quantity
class PhysicalPropertyTest(unittest.TestCase):
@@ -38,66 +38,66 @@ class PhysicalPropertyTest(unittest.TestCase):
def test_temperature(self):
"""Test temperature property."""
- self.equality(pp.temperature, 293.15, qty(20, "°C"))
- self.equality(pp.temperature, [293.15, 303.15], qty([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, qty(1, "MPa"))
- self.equality(pp.pressure, [1e6, 2e6], qty([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], qty(5, "m/s"))
+ self.equality(pp.velocity.magnitude, [3, 4], Qty(5, "m/s"))
self.equality(
- pp.velocity.magnitude, [[3, 4], [1, 0]], qty([5, 1], "m/s")
+ pp.velocity.magnitude, [[3, 4], [1, 0]], Qty([5, 1], "m/s")
)
def test_displacement(self):
"""Test displacement property."""
- self.equality(pp.displacement.magnitude, [3e-3, 4e-3], qty(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, qty([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, 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], "%"))
+ 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, qty(i, "%"))
- self.equality(pp.strain[i], [eps, eps], qty([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, qty(i, "mm"))
- self.equality(pp.displacement[i], [u, u], qty([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([3, 1, 1, 1, 1, 1]) * 1e6
- self.equality(pp.stress.von_Mises, sig_3D, qty(2, "MPa"))
- self.equality(pp.stress.von_Mises, [sig_3D, sig_3D], qty([2, 2], "MPa"))
+ self.equality(pp.stress.von_Mises, sig_3D, Qty(2, "MPa"))
+ self.equality(pp.stress.von_Mises, [sig_3D, sig_3D], Qty([2, 2], "MPa"))
sig_2D = np.array([2, 1, 1, 1]) * 1e6
- self.equality(pp.stress.von_Mises, sig_2D, qty(1, "MPa"))
+ self.equality(pp.stress.von_Mises, 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, qty(2, "MPa"))
- self.equality(pp.effective_pressure, [sig, sig], qty([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.stress.qp_ratio, sig, qty(-100, "percent"))
- self.equality(pp.stress.qp_ratio, [sig] * 2, qty([-100] * 2, "percent"))
+ self.equality(pp.stress.qp_ratio, sig, Qty(-100, "percent"))
+ self.equality(pp.stress.qp_ratio, [sig] * 2, Qty([-100] * 2, "percent"))
def test_depth_2D(self):
mesh = mesh_mechanics
@@ -135,24 +135,24 @@ class PhysicalPropertyTest(unittest.TestCase):
pp.stress.eigenvectors[i].magnitude(sig), 1.0
)
self.assertGreater(pp.stress.det(sig), 0)
- self.assertGreater(pp.stress.I1(sig), 0)
- self.assertGreater(pp.stress.I2(sig), 0)
- self.assertGreater(pp.stress.I3(sig), 0)
+ self.assertGreater(pp.stress.invariant_1(sig), 0)
+ self.assertGreater(pp.stress.invariant_2(sig), 0)
+ self.assertGreater(pp.stress.invariant_3(sig), 0)
self.assertGreater(pp.stress.mean(sig), 0)
self.assertGreater(pp.stress.deviator.magnitude(sig), 0)
- self.assertGreater(pp.stress.J1(sig), 0)
- self.assertGreater(pp.stress.J2(sig), 0)
- self.assertGreater(pp.stress.J3(sig), 0)
+ self.assertGreater(pp.stress.deviator_invariant_1(sig), 0)
+ self.assertGreater(pp.stress.deviator_invariant_2(sig), 0)
+ self.assertGreater(pp.stress.deviator_invariant_3(sig), 0)
self.assertGreater(pp.stress.octahedral_shear(sig), 0)
def test_simple(self):
"""Test call functionality."""
- self.assertEqual(pp.temperature(273.15, strip_unit=False), qty(0, "°C"))
+ 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")
+ 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")
+ pp.displacement([1, 2, 3], strip_unit=False)[1], Qty(2, "m")
)
def test_values(self):