(WIP) Fully anisotropic PR comments 2

CHANGELOG.md

@@ -13,6 +13,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Adjoint processing is done server side by default, to avoid unnecessary downloading of data.
 - `run_local` and `run_async_local` options in `tidy3d.plugins.adjoint.web` to provide way to run adjoint processing locally.
 - `JaxPolySlab` in `adjoint` plugin, which can track derivatives through its `.vertices`.
+- Fully anisotropic medium class (`FullyAnisotropicMedium`) that allows to simulate materials with permittivity and conductivity tensors oriented arbitrary with respect to simulation grid.
 
 ### Changed
 - Perfect electric conductors (PECs) are now modeled as high-conductivity media in both the frontend and backend mode solvers, and their presence triggers the use of a preconditioner to improve numerical stability and robustness. Consequently, the mode solver provides more accurate eigenvectors and field distributions when PEC structures are present.

tests/test_components/test_medium.py

@@ -361,7 +361,7 @@ def test_fully_anisotropic_media():
     assert np.allclose(m.permittivity, perm)
     assert np.allclose(m.conductivity, cond)
 
-    perm_d, cond_d = m.eps_sigma_diag
+    perm_d, cond_d, _ = m.eps_sigma_diag
 
     assert all(np.isin(np.round(perm_d), np.round(np.diag(perm_diag))))
     assert all(np.isin(np.round(cond_d), np.round(np.diag(cond_diag))))

tests/test_components/test_meshgenerate.py

@@ -665,8 +665,8 @@ def test_shapely_strtree_warnings():
 
 def test_anisotropic_material_meshing():
     """Make sure the largest propagation index defines refinement in all directions."""
-    perm_d = [[3, 0, 0], [0, 2, 0], [0, 0, 1]]
-    cond_d = [[0.2, 0, 0], [0, 0.15, 0], [0, 0, 0.1]]
+    perm_diag = [3, 2, 1]
+    cond_diag = [0.2, 0.15, 0.1]
 
     box = td.Box(
         center=(0, 0, 0),
@@ -675,25 +675,23 @@ def test_anisotropic_material_meshing():
 
     box_iso = td.Structure(
         geometry=box,
-        medium=td.Medium(permittivity=perm_d[0][0], conductivity=cond_d[0][0]),
+        medium=td.Medium(permittivity=perm_diag[0], conductivity=cond_diag[0]),
     )
 
     box_diag = td.Structure(
         geometry=box,
         medium=td.AnisotropicMedium(
-            xx=td.Medium(permittivity=perm_d[0][0], conductivity=cond_d[0][0]),
-            yy=td.Medium(permittivity=perm_d[1][1], conductivity=cond_d[1][1]),
-            zz=td.Medium(permittivity=perm_d[2][2], conductivity=cond_d[2][2]),
+            xx=td.Medium(permittivity=perm_diag[0], conductivity=cond_diag[0]),
+            yy=td.Medium(permittivity=perm_diag[1], conductivity=cond_diag[1]),
+            zz=td.Medium(permittivity=perm_diag[2], conductivity=cond_diag[2]),
         ),
     )
 
     box_full = td.Structure(
         geometry=box,
-        medium=td.FullyAnisotropicMedium.from_diagonal(
-            xx=td.Medium(permittivity=perm_d[0][0], conductivity=cond_d[0][0]),
-            yy=td.Medium(permittivity=perm_d[1][1], conductivity=cond_d[1][1]),
-            zz=td.Medium(permittivity=perm_d[2][2], conductivity=cond_d[2][2]),
-            rotation=td.RotationAroundAxis(axis=(1, 2, 3), angle=1.23),
+        medium=td.FullyAnisotropicMedium(
+            permittivity=np.diag(perm_diag),
+            conductivity=np.diag(cond_diag),
         ),
     )
 

tests/test_components/test_types.py

@@ -1,7 +1,7 @@
 """Tests type definitions."""
 import pytest
 import tidy3d as td
-from tidy3d.components.types import ArrayLike, Complex, constrained_array
+from tidy3d.components.types import ArrayLike, Complex, constrained_array, Tuple
 from tidy3d.components.base import Tidy3dBaseModel
 from tidy3d.exceptions import ValidationError
 import numpy as np
@@ -15,6 +15,15 @@ class S(Tidy3dBaseModel):
     with pytest.raises(pydantic.ValidationError):
         s = S(f=np.array([1.0, 2.0, 3.0]))
 
+    class MyClass(Tidy3dBaseModel):
+        f: constrained_array(ndim=3, shape=(1, 2, 3))
+
+    with pytest.raises(pydantic.ValidationError):
+        _ = MyClass(f=np.ones((2, 2, 3)))
+
+    with pytest.raises(pydantic.ValidationError):
+        _ = MyClass(f=np.ones((1, 2, 3, 4)))
+
 
 def test_schemas():
     class S(Tidy3dBaseModel):
@@ -35,13 +44,15 @@ class MyClass(Tidy3dBaseModel):
         c: constrained_array(dtype=float) = None  # must be float-like
         d: constrained_array(ndim=1, dtype=complex) = None  # 1D complex
         e: ArrayLike
+        f: constrained_array(ndim=3, shape=(1, 2, 3)) = None  # must have certain shape
 
     my_obj = MyClass(
         a=1.0 + 2j,
         b=np.array([[1.0, 2.0]]),
         c=[1, 3.0],
         d=[1.0],
         e=[[[[1.0]]]],
+        f=np.ones((1, 2, 3)),
     )
 
     assert np.all(my_obj.a == [1.0 + 2j])  # scalars converted to list of len 1
@@ -57,6 +68,7 @@ class MyClass(Tidy3dBaseModel):
 
         a: ArrayLike
         b: constrained_array(ndim=1)
+        c: Tuple[ArrayLike, ...]
 
-    c = MyClass(a=[1.0], b=[2.0, 1.0])
+    c = MyClass(a=[1.0], b=[2.0, 1.0], c=([2.0, 1.0]))
     hash(c.json())

tidy3d/__init__.py

@@ -12,7 +12,7 @@
 from .components.medium import Medium, PoleResidue, AnisotropicMedium, PEC, PECMedium, Medium2D
 from .components.medium import Sellmeier, Debye, Drude, Lorentz
 from .components.medium import CustomMedium, FullyAnisotropicMedium
-from .components.medium import RotationAroundAxis
+from .components.transformation import RotationAroundAxis
 
 # structures
 from .components.structure import Structure, MeshOverrideStructure

tidy3d/components/base.py

@@ -66,6 +66,13 @@ class Tidy3dBaseModel(pydantic.BaseModel):
     `Pydantic Models <https://pydantic-docs.helpmanual.io/usage/models/>`_
     """
 
+    def __hash__(self) -> int:
+        """Hash method."""
+        try:
+            return super().__hash__(self)
+        except TypeError:
+            return hash(self.json())
+
     def __init__(self, **kwargs):
         """Init method, includes post-init validators."""
         super().__init__(**kwargs)

tidy3d/components/geometry.py

@@ -24,6 +24,7 @@
 from ..constants import MICROMETER, LARGE_NUMBER, RADIAN, fp_eps, inf
 from .data.dataset import TriangleMeshDataset
 from .data.data_array import TriangleMeshDataArray, DATA_ARRAY_MAP
+from .transformation import RotationAroundAxis
 
 try:
     import trimesh
@@ -644,31 +645,8 @@ def rotate_points(points: Points, axis: Coordinate, angle: float) -> Points:
         angle : float
             Angle of rotation counter-clockwise around the axis (rad).
         """
-
-        if isclose(angle % (2 * np.pi), 0):
-            return points
-
-        # Normalized axis vector components
-        (ux, uy, uz) = axis / np.linalg.norm(axis)
-
-        # General rotation matrix
-        rot_mat = np.zeros((3, 3))
-        cos = np.cos(angle)
-        sin = np.sin(angle)
-        rot_mat[0, 0] = cos + ux**2 * (1 - cos)
-        rot_mat[0, 1] = ux * uy * (1 - cos) - uz * sin
-        rot_mat[0, 2] = ux * uz * (1 - cos) + uy * sin
-        rot_mat[1, 0] = uy * ux * (1 - cos) + uz * sin
-        rot_mat[1, 1] = cos + uy**2 * (1 - cos)
-        rot_mat[1, 2] = uy * uz * (1 - cos) - ux * sin
-        rot_mat[2, 0] = uz * ux * (1 - cos) - uy * sin
-        rot_mat[2, 1] = uz * uy * (1 - cos) + ux * sin
-        rot_mat[2, 2] = cos + uz**2 * (1 - cos)
-
-        if len(points.shape) == 1:
-            return rot_mat @ points
-
-        return np.tensordot(rot_mat, points, axes=1)
+        rotation = RotationAroundAxis(axis=axis, angle=angle)
+        return rotation.rotate_vector(points)
 
     def reflect_points(
         self,