From ab6b43c1b56af241aef9831a730a1ecbc1c2e90e Mon Sep 17 00:00:00 2001
From: FZill <florian.zill@ufz.de>
Date: Wed, 6 Mar 2024 13:32:31 +0100
Subject: [PATCH] [propertylib] improved depth calculation + tests
---
.../howto_meshplotlib/plot_solid_mechanics.py | 2 +-
ogstools/propertylib/mesh_dependent.py | 56 ++++++++++++++-----
tests/test_propertylib.py | 17 ++++++
3 files changed, 60 insertions(+), 15 deletions(-)
diff --git a/docs/examples/howto_meshplotlib/plot_solid_mechanics.py b/docs/examples/howto_meshplotlib/plot_solid_mechanics.py
index 0a717d109..22fd4b904 100644
--- a/docs/examples/howto_meshplotlib/plot_solid_mechanics.py
+++ b/docs/examples/howto_meshplotlib/plot_solid_mechanics.py
@@ -108,7 +108,7 @@ fig = plot(mesh, presets.pressure)
# %%
del mesh.point_data["pressure"]
-mesh["depth"] = mesh_dependent.depth(mesh) + 450 # in m
+mesh["depth"] = mesh_dependent.depth(mesh, use_coords=True)
fig = plot(mesh, "depth")
mesh["pressure"] = mesh_dependent.p_fluid(mesh)
fig = plot(mesh, presets.pressure)
diff --git a/ogstools/propertylib/mesh_dependent.py b/ogstools/propertylib/mesh_dependent.py
index d586feaf6..c26938ba0 100644
--- a/ogstools/propertylib/mesh_dependent.py
+++ b/ogstools/propertylib/mesh_dependent.py
@@ -14,23 +14,51 @@ 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:
+def depth(mesh: pv.UnstructuredGrid, use_coords: bool = False) -> np.ndarray:
"""Return the depth values of the mesh.
- This assumes the y-axes is the vertical axis. For each point, the shortest
- distance to the top boundary is returned.
+ For 2D, the last non-flat dimension is used as the vertical axis, 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.
"""
- # 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))
+ 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
+ if use_coords:
+ return -mesh.points[:, top_id]
+ point_upwards = edges.cell_normals[..., top_id] > 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]
+ if use_coords:
+ return -mesh.points[:, top_id]
+ # prevents inner edge (from holes) to be detected as a top boundary
+ edges = mesh.extract_feature_edges().connectivity(
+ "point_seed", point_ids=[0]
+ )
+ edge_horizontal_extent = [
+ np.diff(np.reshape(cell.bounds, (-1, 2))).ravel()[0]
+ for cell in edges.cell
+ ]
+ 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_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)),
+ 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])
def p_fluid(mesh: pv.UnstructuredGrid) -> PlainQuantity:
diff --git a/tests/test_propertylib.py b/tests/test_propertylib.py
index c94b2565f..2a203cf73 100644
--- a/tests/test_propertylib.py
+++ b/tests/test_propertylib.py
@@ -3,11 +3,13 @@
import unittest
import numpy as np
+import pyvista as pv
from pint.facets.plain import PlainQuantity
from ogstools.meshplotlib.examples import mesh_mechanics
from ogstools.propertylib import presets as pp
from ogstools.propertylib.matrix import Matrix
+from ogstools.propertylib.mesh_dependent import depth
from ogstools.propertylib.property import Scalar, u_reg
from ogstools.propertylib.vector import Vector
@@ -97,6 +99,21 @@ class PhysicalPropertyTest(unittest.TestCase):
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
+ mesh["depth"] = depth(mesh, use_coords=True)
+ # y Axis is vertical axis
+ self.assertTrue(np.all(mesh["depth"] == -mesh.points[..., 1]))
+ mesh["depth"] = depth(mesh)
+ self.assertTrue(np.all(mesh["depth"] < -mesh.points[..., 1]))
+
+ def test_depth_3D(self):
+ mesh = pv.SolidSphere(100, center=(0, 0, -101))
+ mesh["depth"] = depth(mesh, use_coords=True)
+ self.assertTrue(np.all(mesh["depth"] == -mesh.points[..., -1]))
+ mesh["depth"] = depth(mesh)
+ self.assertTrue(np.all(mesh["depth"] < -mesh.points[..., -1]))
+
def test_integrity_criteria(self):
"""Test integrity criteria."""
sig = np.array([4, 1, 2, 1, 1, 1]) * 1e6
--
GitLab