diff --git a/CHANGELOG.md b/CHANGELOG.md index 1ff0e9455..36a415cd5 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +### Added +- The method `Geometry.reflected` can be used to create a reflected copy of any geometry off a plane. As for other transformations, for efficiency, `reflected` `PolySlab` directly returns an updated `PolySlab` object rather than a `Transformed` object, except when the normal of the plane of reflection has a non-zero component along the slab axis, in which case `Transformed` is still returned. + ## [2.8.3] - 2025-04-24 ### Added diff --git a/tests/test_components/test_geometry.py b/tests/test_components/test_geometry.py index 64ed89575..d2fd82ead 100644 --- a/tests/test_components/test_geometry.py +++ b/tests/test_components/test_geometry.py @@ -439,6 +439,7 @@ def test_planar_transform(axis): geo = ( td.Box(size=(3 * axis, 2 * abs(axis - 1), 4 * (2 - axis))) .rotated(2.0, axis) + .reflected((axis, 2 * (axis - 1), 3 * (axis - 2))) .translated(-1, 2, 3) .scaled(1.4, -1.2, 1.3) ) @@ -479,6 +480,30 @@ def test_transforms(): assert len(geo.intersections_plane(x=0)) == 1 assert len(geo.intersections_plane(z=0)) == 3 + # Test reflection of a Box across the XY plane and verify point inclusion. + xyz = (np.array([1, 1, 1, 3]), np.array([1, 1, 1, 3]), np.array([1, -1, -1.5, 3])) + geo = td.Box(center=(1, 1, 1), size=(2, 2, 2)) + assert (geo.inside(*xyz) == (True, False, False, False)).all() + geo = geo.reflected((0, 0, 1)) + assert (geo.inside(*xyz) == (False, True, True, False)).all() + + # Test Sphere multiple reflections not influencing point inclusion. + xyz = (np.array([1, 2, 2, 1]), np.array([2, 1, 2, 3]), np.array([2, -2, -0.5, -2])) + geo = td.Sphere(radius=3.5) + assert (geo.inside(*xyz) == (True, True, True, False)).all() + geo = geo.reflected((2, 3, 1)).reflected((1, 2, 3)) + assert (geo.inside(*xyz) == (True, True, True, False)).all() + + # Test PolySlab reflection across non-axis plane and verify point inclusion. + xyz = (np.array([0, 1.5, -1.5, -1.5]), np.array([0, 1.5, -1.5, -2.5]), np.array([0, 0, 0, 0])) + geo = td.PolySlab( + vertices=[(1, 0), (3, 2), (2, 2), (0, 0), (2, -2), (3, -2)], + slab_bounds=(-1, 1), + ) + assert (geo.inside(*xyz) == (True, True, False, False)).all() + geo = geo.reflected((1, 1, 0)) + assert (geo.inside(*xyz) == (True, False, True, True)).all() + def test_polyslab_transforms(): # More tests on PolySlab tranforms matching direct Transformed @@ -501,6 +526,10 @@ def test_polyslab_transforms(): geo = geo.rotated(0.3, (0, -0.2, 0)) assert geo.type != geo_trans.type assert np.allclose(geo.inside(*xyz), geo_trans.inside(*xyz)) + geo_trans = td.Transformed(geometry=geo, transform=td.Transformed.reflection((1, 0, 2))) + geo = geo.reflected((1, 0, 2)) + assert geo.type != geo_trans.type + assert np.allclose(geo.inside(*xyz), geo_trans.inside(*xyz)) def test_general_rotation(): @@ -514,6 +543,21 @@ def test_general_rotation(): assert np.allclose(td.Transformed.rotation(0.3, 2), td.Transformed.rotation(-0.3, [0, 0, -4])) +def test_general_reflection(): + # Magnitude of normal direction does not affect the transformation. + assert np.allclose(td.Transformed.reflection((1, 1, 0)), td.Transformed.reflection((5, 5, 0))) + assert np.allclose(td.Transformed.reflection((1, 0, 1)), td.Transformed.reflection((5, 0, 5))) + assert np.allclose(td.Transformed.reflection((0, 1, 1)), td.Transformed.reflection((0, 5, 5))) + # Negative normal direction means the same transformation. + assert np.allclose(td.Transformed.reflection((1, 1, 0)), td.Transformed.reflection((-1, -1, 0))) + assert np.allclose(td.Transformed.reflection((1, 0, 1)), td.Transformed.reflection((-1, 0, -1))) + assert np.allclose(td.Transformed.reflection((0, 1, 1)), td.Transformed.reflection((0, -1, -1))) + # Magnitude and sign of normal direction does not affect the transformation. + assert np.allclose(td.Transformed.reflection((1, 1, 0)), td.Transformed.reflection((-5, -5, 0))) + assert np.allclose(td.Transformed.reflection((1, 0, 1)), td.Transformed.reflection((-5, 0, -5))) + assert np.allclose(td.Transformed.reflection((0, 1, 1)), td.Transformed.reflection((0, -5, -5))) + + def test_flattening(): flat = list( flatten_groups( diff --git a/tests/test_components/test_structure.py b/tests/test_components/test_structure.py index 151a8634b..5deea90b4 100644 --- a/tests/test_components/test_structure.py +++ b/tests/test_components/test_structure.py @@ -138,6 +138,31 @@ def test_invalid_polyslab(axis): geo7 = td.GeometryGroup(geometries=[(ps - box).rotated(np.pi / 4, j)]).rotated(-np.pi / 4, j) _ = td.Structure(geometry=geo7, medium=medium) + # normal with zero component along the slab axis + n1 = (axis, 2 * (axis - 1), 3 * (axis - 2)) + # normal with non-zero component along the slab axis + n2 = (2 * (axis - 1), 3 * (axis - 2), axis) + + geo8 = ps.reflected(n1) + _ = td.Structure(geometry=geo8, medium=medium) + + geo9 = ps.reflected(n2).scaled(2, 2, 2).translated(-1, 0.5, 2).reflected(n2) + _ = td.Structure(geometry=geo9, medium=medium) + + geo10 = ps.reflected(n1) + _ = td.Structure(geometry=geo10, medium=medium) + + geo11 = ps.reflected(n2) + with pytest.raises(pd.ValidationError): + _ = td.Structure(geometry=geo11, medium=medium) + + geo12 = td.GeometryGroup(geometries=[ps]).reflected(n2) + with pytest.raises(pd.ValidationError): + _ = td.Structure(geometry=geo12, medium=medium) + + geo13 = td.GeometryGroup(geometries=[(ps - box).reflected(n2)]).reflected(n2) + _ = td.Structure(geometry=geo13, medium=medium) + def test_validation_of_structures_with_2d_materials(): med2d = td.Medium2D(ss=td.PEC, tt=td.PEC) diff --git a/tidy3d/components/geometry/base.py b/tidy3d/components/geometry/base.py index 6edbfc8e6..ddc0955b7 100644 --- a/tidy3d/components/geometry/base.py +++ b/tidy3d/components/geometry/base.py @@ -30,7 +30,7 @@ from ..autograd import AutogradFieldMap, TracedCoordinate, TracedSize, get_static from ..autograd.derivative_utils import DerivativeInfo, integrate_within_bounds from ..base import Tidy3dBaseModel, cached_property -from ..transformation import RotationAroundAxis +from ..transformation import ReflectionFromPlane, RotationAroundAxis from ..types import ( ArrayFloat2D, ArrayFloat3D, @@ -975,6 +975,22 @@ def rotated(self, angle: float, axis: Union[Axis, Coordinate]) -> Geometry: """ return Transformed(geometry=self, transform=Transformed.rotation(angle, axis)) + def reflected(self, normal: Coordinate) -> Geometry: + """Return a reflected copy of this geometry. + + Parameters + ---------- + normal : Tuple[float, float, float] + The 3D normal vector of the plane of reflection. The plane is assumed + to pass through the origin (0,0,0). + + Returns + ------- + :class:`Geometry` + Reflected copy of this geometry. + """ + return Transformed(geometry=self, transform=Transformed.reflection(normal)) + """ Field and coordinate transformations """ @staticmethod @@ -2877,6 +2893,25 @@ def rotation(angle: float, axis: Union[Axis, Coordinate]) -> MatrixReal4x4: transform[:3, :3] = RotationAroundAxis(angle=angle, axis=axis).matrix return transform + @staticmethod + def reflection(normal: Coordinate) -> MatrixReal4x4: + """Return a reflection matrix. + + Parameters + ---------- + normal : Tuple[float, float, float] + Normal of the plane of reflection. + + Returns + ------- + numpy.ndarray + Transform matrix with shape (4, 4). + """ + + transform = np.eye(4) + transform[:3, :3] = ReflectionFromPlane(normal=normal).matrix + return transform + @staticmethod def preserves_axis(transform: MatrixReal4x4, axis: Axis) -> bool: """Indicate if the transform preserves the orientation of a given axis. diff --git a/tidy3d/components/geometry/polyslab.py b/tidy3d/components/geometry/polyslab.py index 5fb532361..f9f7e840e 100644 --- a/tidy3d/components/geometry/polyslab.py +++ b/tidy3d/components/geometry/polyslab.py @@ -19,7 +19,7 @@ from ..autograd.derivative_utils import DerivativeInfo, DerivativeSurfaceMesh from ..autograd.types import TracedFloat from ..base import cached_property, skip_if_fields_missing -from ..transformation import RotationAroundAxis +from ..transformation import ReflectionFromPlane, RotationAroundAxis from ..types import ( ArrayFloat1D, ArrayFloat2D, @@ -1726,6 +1726,33 @@ def rotated(self, angle: float, axis: Union[Axis, Coordinate]) -> PolySlab: return super().rotated(angle=angle, axis=axis) + def reflected(self, normal: Coordinate) -> PolySlab: + """Return a reflected copy of this geometry. + + Parameters + ---------- + normal : Tuple[float, float, float] + The 3D normal vector of the plane of reflection. The plane is assumed + to pass through the origin (0,0,0). + + Returns + ------- + ------- + :class:`PolySlab` + Reflected copy of this ``PolySlab``. + """ + if math.isclose(normal[self.axis], 0): + _, plane_axs = self.pop_axis((0, 1, 2), self.axis) + verts_3d = np.zeros((3, self.vertices.shape[0])) + verts_3d[plane_axs[0], :] = self.vertices[:, 0] + verts_3d[plane_axs[1], :] = self.vertices[:, 1] + reflection = ReflectionFromPlane(normal=normal) + reflected_vertices = reflection.reflect_vector(verts_3d) + reflected_vertices = reflected_vertices[plane_axs, :].T + return self.updated_copy(vertices=reflected_vertices) + + return super().reflected(normal=normal) + class ComplexPolySlabBase(PolySlab): """Interface for dividing a complex polyslab where self-intersecting polygon can diff --git a/tidy3d/components/transformation.py b/tidy3d/components/transformation.py index 48de4621f..ea7c1ee49 100644 --- a/tidy3d/components/transformation.py +++ b/tidy3d/components/transformation.py @@ -127,4 +127,78 @@ def matrix(self) -> TensorReal: return R +class AbstractReflection(ABC, Tidy3dBaseModel): + """Abstract reflection of vectors and tensors.""" + + @cached_property + @abstractmethod + def matrix(self) -> TensorReal: + """Reflection matrix.""" + + def reflect_vector(self, vector: ArrayFloat2D) -> ArrayFloat2D: + """Reflect a vector/point or a list of vectors/points. + + Parameters + ---------- + vector : ArrayLike[float] + Array of shape ``(3, ...)``. + + Returns + ------- + Coordinate + Reflected vector. + """ + + if len(vector.shape) == 1: + return self.matrix @ vector + + return np.tensordot(self.matrix, vector, axes=1) + + def reflect_tensor(self, tensor: TensorReal) -> TensorReal: + """Reflect a tensor. + + Parameters + ---------- + tensor : ArrayLike[float] + Array of shape ``(3, 3)``. + + Returns + ------- + TensorReal + Reflected tensor. + """ + + return np.matmul(self.matrix, np.matmul(tensor, self.matrix.T)) + + +class ReflectionFromPlane(AbstractReflection): + """Reflection of vectors and tensors around a given vector.""" + + normal: Coordinate = pd.Field( + (1, 0, 0), + title="Normal of the reflecting plane", + description="A vector that specifies the normal of the plane of reflection", + ) + + @pd.validator("normal") + def _guarantee_nonzero_normal(cls, val): + norm = np.linalg.norm(val) + if np.isclose(norm, 0): + raise ValidationError( + "The norm of vector 'normal' cannot be zero. Please provide a proper normal vector." + ) + return val + + @cached_property + def matrix(self) -> TensorReal: + """Reflection matrix.""" + + norm = np.linalg.norm(self.normal) + n = self.normal / norm + R = np.eye(3) - 2 * np.outer(n, n) + + return R + + RotationType = Union[RotationAroundAxis] +ReflectionType = Union[ReflectionFromPlane]