[geom] Add Sphere radius/volume conversions

phi/geom/_box.py

@@ -457,15 +457,15 @@ def shape(self):
 
     @property
     def size(self):
-        return 2 * self.half_size
+        return 2 * self._half_size
 
     @property
     def lower(self):
-        return self.center - self.half_size
+        return self._center - self._half_size
 
     @property
     def upper(self):
-        return self.center + self.half_size
+        return self._center + self._half_size
 
     @property
     def rotation_matrix(self) -> Optional[Tensor]:

phi/geom/_sphere.py

@@ -1,7 +1,7 @@
 from typing import Union, Dict, Tuple
 
 from phi import math
-from phiml.math import Shape, dual
+from phiml.math import Shape, dual, PI
 from ._geom import Geometry, _keep_vector, NO_GEOMETRY
 from ..math import wrap, Tensor, expand
 from ..math.magic import slicing_dict
@@ -16,6 +16,8 @@ class Sphere(Geometry):
     def __init__(self,
                  center: Tensor = None,
                  radius: Union[float, Tensor] = None,
+                 volume: Union[float, Tensor] = None,
+                 radius_variable=True,
                  **center_: Union[float, Tensor]):
         """
         Args:
@@ -31,8 +33,12 @@ def __init__(self,
             self._center = center
         else:
             self._center = wrap(tuple(center_.values()), math.channel(vector=tuple(center_.keys())))
-        assert radius is not None, "radius must be specified."
-        self._radius = wrap(radius)
+        if radius is None:
+            assert volume is not None, f"Either radius or volume must be specified but got neither."
+            self._radius = Sphere.radius_from_volume(volume, self._center.vector.size)
+        else:
+            self._radius = wrap(radius)
+        self._radius_variable = radius_variable
         assert 'vector' not in self._radius.shape, f"Sphere radius must not vary along vector but got {radius}"
 
     @property
@@ -51,17 +57,32 @@ def center(self):
 
     @property
     def volume(self) -> math.Tensor:
-        if self.spatial_rank == 1:
-            return 2 * self._radius
-        elif self.spatial_rank == 2:
-            return math.PI * self._radius ** 2
-        elif self.spatial_rank == 3:
-            return 4 / 3 * math.PI * self._radius ** 3
+        return Sphere.volume_from_radius(self._radius, self.spatial_rank)
+
+    @staticmethod
+    def volume_from_radius(radius: Union[float, Tensor], spatial_rank: int):
+        if spatial_rank == 1:
+            return 2 * radius
+        elif spatial_rank == 2:
+            return PI * radius ** 2
+        elif spatial_rank == 3:
+            return 4/3 * PI * radius ** 3
         else:
-            raise NotImplementedError()
+            raise NotImplementedError(f"spatial_rank>3 not supported, got {spatial_rank}")
             # n = self.spatial_rank
             # return math.pi ** (n // 2) / math.faculty(math.ceil(n / 2)) * self._radius ** n
 
+    @staticmethod
+    def radius_from_volume(volume: Union[float, Tensor], spatial_rank: int):
+        if spatial_rank == 1:
+            return volume / 2
+        elif spatial_rank == 2:
+            return math.sqrt(volume / PI)
+        elif spatial_rank == 3:
+            return (.75 / PI * volume) ** (1/3)
+        else:
+            raise NotImplementedError(f"spatial_rank>3 not supported, got {spatial_rank}")
+
     def lies_inside(self, location):
         distance_squared = math.sum((location - self.center) ** 2, dim='vector')
         return math.any(distance_squared <= self.radius ** 2, self.shape.instance)  # union for instance dimensions
@@ -93,23 +114,23 @@ def bounding_half_extent(self):
         return expand(self.radius, self._center.shape.only('vector'))
 
     def at(self, center: Tensor) -> 'Geometry':
-        return Sphere(center, self._radius)
+        return Sphere(center, self._radius, radius_variable=self._radius_variable)
 
     def rotated(self, angle):
         return self
 
     def scaled(self, factor: Union[float, Tensor]) -> 'Geometry':
-        return Sphere(self.center, self.radius * factor)
+        return Sphere(self.center, self.radius * factor, radius_variable=self._radius_variable)
 
     def __variable_attrs__(self):
-        return '_center', '_radius'
+        return ('_center', '_radius') if self._radius_variable else ('_center',)
 
     def __value_attrs__(self):
-        return '_center',
+        return ()
 
     def __getitem__(self, item):
         item = slicing_dict(self, item)
-        return Sphere(self._center[_keep_vector(item)], self._radius[item])
+        return Sphere(self._center[_keep_vector(item)], self._radius[item], radius_variable=self._radius_variable)
 
     def __hash__(self):
         return hash(self._center) + hash(self._radius)

phi/physics/sph.py

@@ -10,7 +10,7 @@
 from phi import math
 from phi.field import Field
 from phi.math import Tensor, pairwise_distances, vec_length, Shape, non_channel, dual, where, PI
-from phi.geom import Geometry, Graph, Box
+from phi.geom import Geometry, Graph, Box, Sphere
 from phiml.math import channel, stack, vec, concat, expand, clip
 
 _DEFAULT_DESIRED_NEIGHBORS = {
@@ -80,12 +80,7 @@ def _get_support_radius(volume: Tensor, desired_neighbors: float, spatial_rank:
     Returns:
         Support radius, i.e. neighbor search cutoff.
     """
-    if spatial_rank == 1:
-        return desired_neighbors * .5 * volume  # N(h) = 2 h / v
-    elif spatial_rank == 2:
-        return math.sqrt(desired_neighbors / math.PI * volume)  # N(h) = 𝜋 h² / v
-    else:
-        return (desired_neighbors / math.PI * .75 * volume) ** (1/3)  # N(h) = 4/3 𝜋 h³ / v
+    return Sphere.radius_from_volume(volume * desired_neighbors, spatial_rank)
 
 
 def expected_neighbors(volume: Tensor, support_radius, spatial_rank: int):
@@ -100,12 +95,7 @@ def expected_neighbors(volume: Tensor, support_radius, spatial_rank: int):
     Returns:
         Number of expected neighbors.
     """
-    if spatial_rank == 1:
-        return 2 * support_radius / volume
-    elif spatial_rank == 2:
-        return PI * support_radius**2 / volume
-    else:
-        return 4/3 * PI * support_radius**3 / volume
+    return Sphere.volume_from_radius(support_radius, spatial_rank) / volume
 
 
 def evaluate_kernel(delta, distance, h, spatial_rank: int, kernel: str, types: Sequence[str] = ('kernel',)) -> Dict[str, Tensor]: