Merge pull request #948 from tanocalogero/main

Fixed bugs in EllipticalCylinder and indices_in_cylinder.

Author: zerothi

changes/947.fix.rst

@@ -0,0 +1 @@
+Fixed finding coordinates within elliptical cylinders

src/sisl/_indices.pyx

@@ -146,16 +146,17 @@ def indices(ints_st[::1] element, ints_st[::1] test_element, ints_st offset=0,
 @cython.wraparound(False)
 @cython.initializedcheck(False)
 def indices_in_cylinder(floats_st[:, ::1] dxyz, const floats_st R, const floats_st h):
-    """ Indices for all coordinates that are within a cylinde radius `R` and height `h`
+    """ Indices for all coordinates that are within a cylinder radius `R` and height `h`
 
     Parameters
     ----------
     dxyz :
-       coordinates centered around the cylinder
+       coordinates centered around the cylinder.
+       The last axis is the cylinder height.
     R :
-       radius of cylinder to check
+       radius of cylinder to check.
     h :
-       height of cylinder to check
+       height of cylinder to check.
 
     Returns
     -------
@@ -168,25 +169,23 @@ def indices_in_cylinder(floats_st[:, ::1] dxyz, const floats_st R, const floats_
     cdef ndarray[int32_t] IDX = np.empty([n], dtype=np.int32)
     cdef int[::1] idx = IDX
 
-    cdef floats_st R2 = R * R
-    cdef floats_st L2
+    cdef floats_st R2
+    cdef floats_st hhalve, L2
     cdef Py_ssize_t i, j, m
-    cdef bint skip
+
+    # Handle radius input
+    R2 = R * R
+    hhalve = h / 2
 
     # Reset number of elements
     m = 0
 
     with nogil:
         for i in range(n):
-            skip = 0
-            for j in range(nxyz):
-                skip |= dxyz[i, j] > R
-            if skip or dxyz[i, nxyz] > h: continue
-
-            L2 = 0.
-            for j in range(nxyz):
-                L2 += dxyz[i, j] * dxyz[i, j]
-            if L2 > R2: continue
+            if dxyz[i,nxyz] > hhalve or dxyz[i,nxyz] < -hhalve: continue
+            # Calculate the distance of the circle
+            L2 = (dxyz[i, 0]*dxyz[i, 0])/R2 + (dxyz[i, 1]*dxyz[i, 1])/R2
+            if L2 > 1.0: continue
             idx[m] = <int> i
             m += 1
 

src/sisl/shape/_cylinder.py

@@ -97,7 +97,7 @@ def __init__(
         self._h = h
 
     def copy(self) -> Self:
-        return self.__class__(self.radial_vector, self.height, self.center)
+        return self.__class__(self.radial_vector, self.height, center=self.center)
 
     @property
     def volume(self) -> float:
@@ -140,7 +140,7 @@ def scale(self, scale: SeqOrScalarFloat) -> Self:
         else:
             v = self._v * scale
             h = self._h * scale
-        return self.__class__(v, h, self.center)
+        return self.__class__(v, h, center=self.center)
 
     def expand(self, radius: SeqOrScalarFloat) -> Self:
         """Expand elliptical cylinder by a constant value along each vector and height
@@ -163,7 +163,7 @@ def expand(self, radius: SeqOrScalarFloat) -> Self:
             raise ValueError(
                 f"{self.__class__.__name__}.expand requires the radius to be either (1,) or (3,)"
             )
-        return self.__class__([v0, v1], h, self.center)
+        return self.__class__([v0, v1], h, center=self.center)
 
     @deprecate_argument(
         "tol",
@@ -191,7 +191,7 @@ def within_index(self, other, rtol: float = 1.0e-8):
         # Get indices where we should do the more
         # expensive exact check of being inside shape
         # I.e. this reduces the search space to the box
-        return indices_in_cylinder(tmp, 1.0 + rtol, 1.0 + rtol)
+        return indices_in_cylinder(tmp, 1.0 + rtol)
 
     @deprecation(
         "toSphere is deprecated, use shape.to.Sphere(...) instead.", "0.15", "0.17"

src/sisl/shape/tests/test_cylinder.py

@@ -39,10 +39,24 @@ def test_create_ellipticalcylinder():
 
 def test_ellipticalcylinder_within():
     el = EllipticalCylinder(1.0, 1.0)
-    # center of cylinder
-    assert el.within_index([0, 0, 0])[0] == 0
+    # points in an ellipsis
+    points = [
+        [0, 0, 0],
+        [0, 0, 0.5],
+        [0, 0, -0.5],
+        [1, 0, -0.5],
+        [0, 1, -0.5],
+    ]
+    assert len(el.within_index(points)) == len(points)
+
     # should not be in a circle
-    assert el.within_index([0.2, 0.2, 0.9])[0] == 0
+    points = [
+        [0, 0, 0.6],
+        [0, 0, -0.6],
+        [0.2, 0.2, 0.9],
+        [0.2, 0.2, -0.9],
+    ]
+    assert len(el.within_index(points)) == 0
 
 
 def test_tosphere():