[geom] Fix push(), Sphere.approximate_closest_surface()

phi/geom/_geom_ops.py

@@ -245,7 +245,7 @@ def expel(geometry: Geometry,
     signed_distance, vec_to_surface, *_ = geometry.approximate_closest_surface(location)
     if not invert:
         shift_amount = math.maximum(0, min_separation - signed_distance)  # expel
-        direction = vec_to_surface / -signed_distance  # this always points outward
+        direction = math.safe_div(vec_to_surface, -signed_distance)  # this always points outward
     else:
         raise NotImplementedError
     return location + direction * shift_amount

phi/geom/_sphere.py

@@ -1,7 +1,7 @@
 from typing import Union, Dict, Tuple
 
 from phi import math
-from phiml.math import Shape, dual, PI
+from phiml.math import Shape, dual, PI, non_channel
 from ._geom import Geometry, _keep_vector, NO_GEOMETRY
 from ..math import wrap, Tensor, expand
 from ..math.magic import slicing_dict
@@ -104,6 +104,17 @@ def approximate_signed_distance(self, location: Union[Tensor, tuple]):
         distance = math.sqrt(distance_squared)
         return math.min(distance - self.radius, self.shape.instance)  # union for instance dimensions
 
+    def approximate_closest_surface(self, location: Tensor) -> Tuple[Tensor, Tensor, Tensor, Tensor, Tensor]:
+        center_delta = location - self.center
+        center_dist = math.vec_length(center_delta)
+        sgn_dist = center_dist - self.radius
+        normal = center_delta / center_dist
+        surface_pos = self.center + self.radius * normal
+        delta = surface_pos - location
+        face_index = expand(0, non_channel(location))
+        offset = normal.vector @ surface_pos.vector
+        return sgn_dist, delta, normal, offset, face_index
+
     def sample_uniform(self, *shape: math.Shape):
         raise NotImplementedError('Not yet implemented')  # ToDo