Added PerturbedDroplet3DAxisSym (#45)

Added special axisymmetric, perturbed droplet in 3D

droplets/droplets.py

@@ -12,6 +12,7 @@
    DiffuseDroplet
    PerturbedDroplet2D
    PerturbedDroplet3D
+   PerturbedDroplet3DAxisSym
 
 
 Inheritance structure of the classes:
@@ -854,6 +855,16 @@ def _get_phase_field(self, grid: GridBase, dtype=np.double) -> np.ndarray:
 
         return result.astype(dtype)
 
+    def _get_mpl_patch(self, dim=None, *, color=None, **kwargs):
+        """return the patch representing the droplet for plotting
+
+        Args:
+            dim (int, optional):
+                The dimension in which the data is plotted. If omitted, the actual
+                physical dimension is assumed
+        """
+        raise NotImplementedError(f"Plotting {self.__class__.__name__} not implemented")
+
 
 class PerturbedDroplet2D(PerturbedDropletBase):
     r"""Represents a single droplet in two dimensions with a perturbed shape
@@ -1193,15 +1204,79 @@ def volume_approx(self) -> float:
             volume += self.amplitudes[0] * 2 * np.sqrt(np.pi) * self.radius**2
         return volume
 
-    def _get_mpl_patch(self, dim=None, *, color=None, **kwargs):
-        """return the patch representing the droplet for plotting
+
+class PerturbedDroplet3DAxisSym(PerturbedDropletBase):
+    r"""Represents a droplet axisymmetrically perturbed shape in three dimensions
+
+    The shape is described using the distance :math:`R(\theta)` of the interface from
+    the origin as a function of the azimuthal angle :math:`\theta`, while polar symmetry
+    is assumed. This function is developed as a truncated series of spherical harmonics
+    :math:`Y_{l,m}(\theta, 0)`:
+
+    .. math::
+        R(\theta) = R_0 \left[1 + \sum_{l=1}^{N_l}
+                                                \epsilon_{l} Y_{l,0}(\theta, 0) \right]
+
+    where :math:`N_l` is the number of perturbation modes considered, which is deduced
+    from the length of the `amplitudes` array.
+    """
+
+    dim = 3
+
+    __slots__ = ["data"]
+
+    def check_data(self):
+        """method that checks the validity and consistency of self.data"""
+        super().check_data()
+        if not np.allclose(self.position[:2], 0):
+            raise ValueError("Droplet must lie on z-axis")
+
+    @preserve_scalars
+    def interface_distance(self, θ: np.ndarray) -> np.ndarray:  # type: ignore
+        r"""calculates the distance of the droplet interface to the origin
 
         Args:
-            dim (int, optional):
-                The dimension in which the data is plotted. If omitted, the actual
-                physical dimension is assumed
+            θ (float or :class:`~np.ndarray`):
+                Azimuthal angle (in :math:`[0, \pi]`)
+
+        Returns:
+            Array with distances of the interfacial points associated with the angles
+        """
+        dist = np.ones(θ.shape, dtype=np.double)
+        for order, a in enumerate(self.amplitudes, 1):  # skip zero-th mode!
+            if a != 0:
+                dist += a * spherical.spherical_harmonic_symmetric(order, θ)  # type: ignore
+        return self.radius * dist
+
+    @preserve_scalars
+    def interface_curvature(self, θ: np.ndarray) -> np.ndarray:  # type: ignore
+        r"""calculates the mean curvature of the interface of the droplet
+
+        For simplicity, the effect of the perturbations are only included to
+        linear order in the perturbation amplitudes :math:`\epsilon_{l,m}`.
+
+        Args:
+            θ (float or :class:`~np.ndarray`):
+                Azimuthal angle (in :math:`[0, \pi]`)
+
+        Returns:
+            Array with curvature at the interfacial points associated with the angles
         """
-        raise NotImplementedError("Plotting PerturbedDroplet3D is not implemented")
+        Yl = spherical.spherical_harmonic_symmetric
+        correction = 0
+        for order, a in enumerate(self.amplitudes, 1):  # skip zero-th mode!
+            if a != 0:
+                hl = (order**2 + order - 2) / 2
+                correction = a * hl * Yl(order, θ)  # type: ignore
+        return 1 / self.radius + correction / self.radius**2  # type: ignore
+
+    @property
+    def volume_approx(self) -> float:
+        """float: approximate volume to linear order in the perturbation"""
+        volume = spherical.volume_from_radius(self.radius, 3)
+        if len(self.amplitudes) > 0:
+            volume += self.amplitudes[0] * 2 * np.sqrt(np.pi) * self.radius**2
+        return volume
 
 
 def droplet_from_data(droplet_class: str, data: np.ndarray) -> DropletBase:
@@ -1267,4 +1342,5 @@ def get_triangulation(self, num_est: int = 1) -> Dict[str, Any]:
     "DiffuseDroplet",
     "PerturbedDroplet2D",
     "PerturbedDroplet3D",
+    "PerturbedDroplet3DAxisSym",
 ]

droplets/image_analysis.py

@@ -41,6 +41,7 @@
     DiffuseDroplet,
     PerturbedDroplet2D,
     PerturbedDroplet3D,
+    PerturbedDroplet3DAxisSym,
     SphericalDroplet,
 )
 from .emulsions import Emulsion
@@ -417,7 +418,10 @@ def locate_droplets(
             if dim == 2:
                 droplet_class = PerturbedDroplet2D
             elif dim == 3:
-                droplet_class = PerturbedDroplet3D
+                if isinstance(phase_field.grid, CylindricalSymGrid):
+                    droplet_class = PerturbedDroplet3DAxisSym
+                else:
+                    droplet_class = PerturbedDroplet3D
             else:
                 raise NotImplementedError(f"Dimension {dim} is not supported")
             args["amplitudes"] = np.zeros(modes)

droplets/tools/spherical.py

@@ -90,8 +90,10 @@ def radius_from_volume(volume: TNumArr, dim: int) -> TNumArr:
     """Return the radius of a sphere with a given volume
 
     Args:
-        volume (float or :class:`~numpy.ndarray`): Volume of the sphere
-        dim (int): Dimension of the space
+        volume (float or :class:`~numpy.ndarray`):
+            Volume of the sphere
+        dim (int):
+            Dimension of the space
 
     Returns:
         float or :class:`~numpy.ndarray`: Radius of the sphere
@@ -110,7 +112,8 @@ def make_radius_from_volume_compiled(dim: int) -> Callable[[TNumArr], TNumArr]:
     """Return a function calculating the radius of a sphere with a given volume
 
     Args:
-        dim (int): Dimension of the space. If omitted, a general function is returned
+        dim (int):
+            Dimension of the space
 
     Returns:
         function: A function that takes a volume and returns the radius
@@ -159,7 +162,8 @@ def make_volume_from_radius_compiled(dim: int) -> Callable[[TNumArr], TNumArr]:
     """Return a function calculating the volume of a sphere with a given radius
 
     Args:
-        dim (int): Dimension of the space
+        dim (int):
+            Dimension of the space
 
     Returns:
         function: A function that takes a radius and returns the volume
@@ -208,8 +212,10 @@ def surface_from_radius(radius: TNumArr, dim: int) -> TNumArr:
     """Return the surface area of a sphere with a given radius
 
     Args:
-        radius (float or :class:`~numpy.ndarray`): Radius of the sphere
-        dim (int): Dimension of the space
+        radius (float or :class:`~numpy.ndarray`):
+            Radius of the sphere
+        dim (int):
+            Dimension of the space
 
     Returns:
         float or :class:`~numpy.ndarray`: Surface area of the sphere
@@ -233,8 +239,10 @@ def radius_from_surface(surface: TNumArr, dim: int) -> TNumArr:
     """Return the radius of a sphere with a given surface area
 
     Args:
-        surface (float or :class:`~numpy.ndarray`): Surface area of the sphere
-        dim (int): Dimension of the space
+        surface (float or :class:`~numpy.ndarray`):
+            Surface area of the sphere
+        dim (int):
+            Dimension of the space
 
     Returns:
         float or :class:`~numpy.ndarray`: Radius of the sphere
@@ -345,8 +353,10 @@ def spherical_index_k(degree: int, order: int = 0) -> int:
     """returns the mode `k` from the degree `degree` and order `order`
 
     Args:
-        degree (int): Degree of the spherical harmonics
-        order (int): Order of the spherical harmonics
+        degree (int):
+            Degree :math:`l` of the spherical harmonics
+        order (int):
+            Order :math:`m` of the spherical harmonics
 
     Raises:
         ValueError: if `order < -degree` or `order > degree`
@@ -363,7 +373,8 @@ def spherical_index_lm(k: int) -> Tuple[int, int]:
     """returns the degree `l` and the order `m` from the mode `k`
 
     Args:
-        k (int): The combined index for the spherical harmonics
+        k (int):
+            The combined index for the spherical harmonics
 
     Returns:
         tuple: The degree `l` and order `m` of the spherical harmonics
@@ -379,7 +390,8 @@ def spherical_index_count(l: int) -> int:
     The returned value is one less than the maximal mode `k` required.
 
     Args:
-        l (int): Maximal degree of the spherical harmonics
+        l (int):
+            Maximal degree of the spherical harmonics
 
     Returns:
         int: The number of modes
@@ -391,7 +403,11 @@ def spherical_index_count_optimal(k_count: int) -> bool:
     """checks whether the modes captures all orders for maximal degree
 
     Args:
-        k_count (int): The number of modes considered
+        k_count (int):
+            The number of modes considered
+
+    Returns:
+        bool: indiciates whether `k_count` is optimally chosen.
     """
     is_square = bool(int(np.sqrt(k_count) + 0.5) ** 2 == k_count)
     return is_square
@@ -401,9 +417,10 @@ def spherical_harmonic_symmetric(degree: int, θ: float) -> float:
     r"""axisymmetric spherical harmonics with degree `degree`, so `m=0`.
 
     Args:
-        degree (int): Degree of the spherical harmonics
-        θ (float): Azimuthal angle at which the spherical harmonics is
-            evaluated (in :math:`[0, \pi]`)
+        degree (int):
+            Degree of the spherical harmonics
+        θ (float):
+            Azimuthal angle at which function is evaluated (in :math:`[0, \pi]`)
 
     Returns:
         float: The value of the spherical harmonics
@@ -417,12 +434,14 @@ def spherical_harmonic_real(degree: int, order: int, θ: float, φ: float) -> fl
     r"""real spherical harmonics of degree l and order m
 
     Args:
-        degree (int): Degree :math:`l` of the spherical harmonics
-        order (int): Order :math:`m` of the spherical harmonics
-        θ (float): Azimuthal angle (in :math:`[0, \pi]`) at which the
-            spherical harmonics is evaluated.
-        φ (float): Polar angle (in :math:`[0, 2\pi]`) at which the spherical
-            harmonics is evaluated.
+        degree (int):
+            Degree :math:`l` of the spherical harmonics
+        order (int):
+            Order :math:`m` of the spherical harmonics
+        θ (float):
+            Azimuthal angle (in :math:`[0, \pi]`) at which fucntion is evaluated.
+        φ (float):
+            Polar angle (in :math:`[0, 2\pi]`) at which function is evaluated.
 
     Returns:
         float: The value of the spherical harmonics
@@ -448,12 +467,12 @@ def spherical_harmonic_real_k(k: int, θ: float, φ: float) -> float:
     r"""real spherical harmonics described by mode k
 
     Args:
-        k (int): Combined index determining the degree and order of the
-            spherical harmonics
-        θ (float): Azimuthal angle (in :math:`[0, \pi]`) at which the
-            spherical harmonics is evaluated.
-        φ (float): Polar angle (in :math:`[0, 2\pi]`) at which the spherical
-            harmonics is evaluated.
+        k (int):
+            Combined index determining the degree and order of the spherical harmonics
+        θ (float):
+            Azimuthal angle (in :math:`[0, \pi]`) at which fucntion is evaluated.
+        φ (float):
+            Polar angle (in :math:`[0, 2\pi]`) at which function is evaluated.
 
     Returns:
         float: The value of the spherical harmonics

tests/test_droplets.py

@@ -64,14 +64,25 @@ def test_perturbed_droplet_2d():
 
 
 def test_perturbed_droplet_3d():
-    """test methods of perturbed droplets in 2d"""
+    """test methods of perturbed droplets in 3d"""
     d = droplets.PerturbedDroplet3D([0, 1, 2], 1, 0.1, [0.0, 0.1, 0.2, 0.3])
     d.volume_approx
     d.interface_distance(0.1, 0.2)
     d.interface_position(0.1, 0.2)
     d.interface_curvature(0.1, 0.2)
 
 
+def test_perturbed_droplet_3d_axis_sym():
+    """test methods of axisymmetrically perturbed droplets in 3d"""
+    d = droplets.PerturbedDroplet3DAxisSym([0, 0, 0], 1, 0.1, [0.0, 0.1])
+    d.volume_approx
+    d.interface_distance(0.1)
+    d.interface_curvature(0.1)
+
+    with pytest.raises(ValueError):
+        droplets.PerturbedDroplet3DAxisSym([0, 1, 0], 1)
+
+
 def test_perturbed_volume():
     """test volume calculation of perturbed droplets"""
     pos = np.random.randn(2)