[geom] Add line_trace(), Mesh distance

phi/geom/__init__.py

@@ -23,5 +23,6 @@
 from ._sdf import SDF, numpy_sdf
 from ._geom_ops import union, intersection
 from ._convert import surface_mesh, as_sdf
+from ._geom_functions import line_trace
 
 __all__ = [key for key in globals().keys() if not key.startswith('_')]

phi/geom/_convert.py

@@ -103,7 +103,7 @@ def surface_mesh(geo: Geometry,
     if isinstance(geo, NoGeometry):
         return mesh_from_numpy([], [], build_faces=False, element_rank=2, build_normals=False)
     if method == 'auto' and isinstance(geo, BaseBox):
-        assert rel_dx is None and abs_dx is None, f"When method='auto', boxes will always use their corners as vertices. rel_dx and abs_dx must be left unspecified."
+        assert rel_dx is None and abs_dx is None, f"When method='auto', boxes will always use their corners as vertices. Leave rel_dx,abs_dx unspecified or pass 'lewiner' or 'lorensen' as method"
         vertices = pack_dims(geo.corners, dual, instance('vertices'))
         corner_count = vertices.vertices.size
         vertices = pack_dims(vertices, instance(geo) + instance('vertices'), instance('vertices'))

phi/geom/_geom.py

@@ -333,12 +333,8 @@ def bounding_radius(self) -> Tensor:
         Returns the radius of a Sphere object that fully encloses this geometry.
         The sphere is centered at the center of this geometry.
 
-        :return: radius of type float
-
-        Args:
-
-        Returns:
-
+        If this geometry consists of multiple parts listed along instance/spatial dims, these dims are retained, giving the bounds of each part.
+        If these dims are not present on the result, all parts are assumed to have the same bounds.
         """
         raise NotImplementedError(self.__class__)
 
@@ -350,8 +346,8 @@ def bounding_half_extent(self) -> Tensor:
         Let the bounding half-extent have value `e` in dimension `d` (`extent[...,d] = e`).
         Then, no point of the geometry lies further away from its center point than `e` along `d` (in both axis directions).
 
-        When this geometry consists of multiple parts listed along instance/spatial dims, these dims are retained, giving the bounds of each part.
-        If these dims are not present, all parts are assumed to have the same bounds.
+        If this geometry consists of multiple parts listed along instance/spatial dims, these dims are retained, giving the bounds of each part.
+        If these dims are not present on the result, all parts are assumed to have the same bounds.
         """
         raise NotImplementedError(self.__class__)
 

phi/geom/_geom_functions.py

@@ -0,0 +1,48 @@
+import warnings
+from typing import Tuple, Optional
+
+from phiml import math
+from phiml.math import Tensor
+
+from ._geom import Geometry
+
+
+def line_trace(geo: Geometry, origin: Tensor, direction: Tensor, max_iter=64, epsilon=None) -> Tuple[Tensor, Tensor, Tensor, Tensor, Optional[Tensor]]:
+    """
+    Trace a line until it hits the surface of `geo`.
+    The surface can be hit either from the outside or the inside.
+
+    Args:
+        geo: `Geometry` that implements `approximate_closest_surface`.
+        origin: Line start location.
+        direction: Unit vector pointing in the line direction.
+        max_iter: Maximum number of steps per line.
+        epsilon: Surface distance tolerance.
+
+    Returns:
+        hit: Whether a surface intersection was found for the line.
+        distance: Distance between the line and the surface.
+        position: Hit location or point until which the line was traced.
+        normal: Surface normal at hit location
+        hit_index: Geometry face index at hit location
+    """
+    if epsilon is None:
+        epsilon = 1e-4 * geo.bounding_box().size.min
+    walked = 0
+    has_hit = False
+    initial_sdf = None
+    for i in range(max_iter):
+        sgn_dist, delta, normal, _, face_index = geo.approximate_closest_surface(origin + walked * direction)
+        initial_sdf = sgn_dist if initial_sdf is None else initial_sdf
+        has_crossed = sgn_dist * initial_sdf < 0
+        intersection = (normal.vector @ delta.vector) / (normal.vector @ direction.vector)
+        intersection = math.where(math.is_nan(intersection), math.INF, intersection)
+        max_walk = math.minimum(abs(sgn_dist), math.where(intersection < -epsilon, math.INF, intersection))
+        max_walk = math.where(has_hit, 0, max_walk)
+        walked += max_walk * .9
+        has_hit = (abs(sgn_dist) <= epsilon) | has_crossed
+        if math.all(has_hit):
+            break
+    else:
+        warnings.warn(f"thickness reached maximum iterations {max_iter}", RuntimeWarning, stacklevel=2)
+    return has_hit, walked, origin + walked * direction, normal, face_index

phi/geom/_mesh.py

@@ -202,8 +202,10 @@ def boundary_connectivity(self) -> Tensor:
 
     @property
     def connectivity(self) -> Tensor:
+        if self._element_connectivity is not None:
+            return self._element_connectivity
         if self._face_areas is None and self._face_normals is None and self._face_centers is None:
-            return math.sparse_tensor(None, None, self.face_shape)
+            return None
         if is_sparse(self._face_areas):
             return tensor_like(self._face_areas, True)
         else:
@@ -328,21 +330,10 @@ def lies_inside(self, location: Tensor) -> Tensor:
 
     def approximate_signed_distance(self, location: Union[Tensor, tuple]) -> Tensor:
         if self.element_rank == 2 and self.spatial_rank == 3:
-            warnings.warn("Use geom.as_sdf(mesh, method='closest-face') for improved performance", RuntimeWarning)
             closest_elem = math.find_closest(self._center, location)
             center = self._center[closest_elem]
-            vertices = self._elements[closest_elem]
-            assert dual(vertices._indices).size == 3, f"signed distance currently only supports triangles"
-            corners = self._vertices.center[vertices._indices]
-            c1, c2, c3 = unstack(corners, dual)
-            v1, v2 = c2 - c1, c3 - c1
-            normal = math.vec_normalize(math.cross_product(v1, v2))
+            normal = self._normals[closest_elem]
             return plane_sgn_dist(center, normal, location)
-            # closest_vertex = math.find_closest(self._vertices.center, location, index_dim=None)
-            # faces = self._elements[{dual: closest_vertex}]   # sparse (elements, locations)
-            # # ToDo write this first for constant number of elements per vertex, then use sparse
-            # tri_vert_ids = self._elements._indices[{instance: faces._indices}]  # ToDo this loses the sparse structure from faces -> Matrix-matrix outer product?
-            # A, B, C = self._vertices.center[tri_vert_ids].vertices.dual
         if self._center is None:
             raise NotImplementedError("Mesh.approximate_signed_distance only available when faces are built.")
         idx = math.find_closest(self._center, location)
@@ -351,6 +342,16 @@ def approximate_signed_distance(self, location: Union[Tensor, tuple]) -> Tensor:
         return math.max(distances, dual)
 
     def approximate_closest_surface(self, location: Tensor) -> Tuple[Tensor, Tensor, Tensor, Tensor, Tensor]:
+        if self.element_rank == 2 and self.spatial_rank == 3:
+            closest_elem = math.find_closest(self._center, location)
+            center = self._center[closest_elem]
+            normal = self._normals[closest_elem]
+            face_size = math.sqrt(self._volume) * 4
+            size = face_size[closest_elem]
+            sgn_dist = plane_sgn_dist(center, normal, location)
+            delta = center - location  # this is not accurate...
+            outward = math.where(abs(sgn_dist) < size, normal, math.normalize(delta))
+            return sgn_dist, delta, outward, None, closest_elem
         # idx = math.find_closest(self._center, location)
         # for i in range(self._max_cell_walk):
         #     idx, leaves_mesh, is_outside, distances, nb_idx = self.cell_walk_towards(location, idx, allow_exit=False)