Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Rotating geometries does not work #1538

Closed
samuelkim16 opened this issue Mar 12, 2024 · 7 comments
Closed

Rotating geometries does not work #1538

samuelkim16 opened this issue Mar 12, 2024 · 7 comments
Assignees

Comments

@samuelkim16
Copy link

samuelkim16 commented Mar 12, 2024

I'm unable to rotate geometries - or rather, the simulation seems to build just fine but I'm unable to plot simulations with rotated geometries. I have not tried running the simulation yet. It seems that the objects are somehow being transformed to on or outside the simulation boundaries, which may be messing up the simulation?

import numpy as np
import tidy3d as td
import matplotlib.pyplot as plt

# Central frequency and bandwidth of pulsed excitation, in Hz
freq0 = 2e14
fwidth = 1e13

# Total time to run in seconds
run_time = 2/fwidth

psource = td.PlaneWave(
    center=(0,0,1.5),
    direction='-',
    size=(td.inf, td.inf, 0),
    source_time = td.GaussianPulse(
        freq0=freq0,
        fwidth=fwidth),
    pol_angle=np.pi/2)

material1 = td.Medium(permittivity=6.)
material2 = td.Medium(permittivity=2.)

box_geom = td.Box(center=[0, 0, 0], size=[2, td.inf, 1])
# box_geom = box_geom.rotated(axis=[0, 0, 1], angle=np.pi/4)

box = td.Structure(
    geometry=box_geom,
    medium=material1
)

# Triangle in the xy-plane with a finite extent in z
equi_tri_verts = [[-1/2, -1/4],
                  [1/2, -1/4],
                  [0, np.sqrt(3)/2 - 1/4]]

poly = td.PolySlab(
    vertices=(2*np.array(equi_tri_verts)).tolist(),
    slab_bounds=(.5, 1.0),
    axis=2
)

poly = poly.rotated(axis=2, angle=np.pi/4)

poly = td.Structure(
    geometry=poly,
    medium=material2)

# Initialize simulation
sim = td.Simulation(size=[4, 4, 4],
                    structures=[box, poly],
                    sources=[psource],
                    run_time=run_time)

# Plot the geometry in the Tidy3D Simulation object
fig, (ax1,ax2,ax3) = plt.subplots(1, 3, tight_layout=True, figsize=(12, 4))

# Waveguide cross-section
ax1 = sim.plot(y=0, ax=ax1)
# Top view
ax2 = sim.plot(z=0, ax=ax2)
ax3 = sim.plot(z=0.8, ax=ax3)

which gives the below error message

16:58:27 Eastern Daylight Time WARNING: Structure at structures[1] was detected 
                               as being less than half of a central wavelength  
                               from a PML on side x-min. To avoid inaccurate    
                               results or divergence, please increase gap       
                               between any structures and PML or fully extend   
                               structure through the pml.                       
                               WARNING: Suppressed 1 WARNING message.           
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[20], line 36
     33 fig, (ax1,ax2) = plt.subplots(1, 2, tight_layout=True, figsize=(10, 4))
     35 # Waveguide cross-section
---> 36 ax1 = sim.plot(y=0, ax=ax1)
     38 # Top view
     39 ax2 = sim.plot(z=0, ax=ax2)

File ~\AppData\Local\anaconda3\envs\asdf\lib\site-packages\tidy3d\components\viz.py:66, in equal_aspect.<locals>._plot(*args, **kwargs)
     63 @wraps(plot)
     64 def _plot(*args, **kwargs) -> Ax:
     65     """New plot function with equal aspect ratio axes returned."""
---> 66     ax = plot(*args, **kwargs)
     67     ax.set_aspect("equal")
     68     return ax

File ~\AppData\Local\anaconda3\envs\asdf\lib\site-packages\tidy3d\components\viz.py:52, in add_ax_if_none.<locals>._plot(*args, **kwargs)
     50     ax = make_ax()
     51     kwargs["ax"] = ax
---> 52 return plot(*args, **kwargs)

File ~\AppData\Local\anaconda3\envs\asdf\lib\site-packages\tidy3d\components\simulation.py:1791, in Simulation.plot(self, x, y, z, ax, source_alpha, monitor_alpha, hlim, vlim, **patch_kwargs)
   1761 """Plot each of simulation's components on a plane defined by one nonzero x,y,z coordinate.
   1762 
   1763 Parameters
   (...)
   1785     The supplied or created matplotlib axes.
   1786 """
   1787 hlim, vlim = Scene._get_plot_lims(
   1788     bounds=self.simulation_bounds, x=x, y=y, z=z, hlim=hlim, vlim=vlim
   1789 )
-> 1791 ax = self.plot_structures(ax=ax, x=x, y=y, z=z, hlim=hlim, vlim=vlim)
   1792 ax = self.plot_sources(ax=ax, x=x, y=y, z=z, hlim=hlim, vlim=vlim, alpha=source_alpha)
   1793 ax = self.plot_monitors(ax=ax, x=x, y=y, z=z, hlim=hlim, vlim=vlim, alpha=monitor_alpha)

File ~\AppData\Local\anaconda3\envs\asdf\lib\site-packages\tidy3d\components\viz.py:66, in equal_aspect.<locals>._plot(*args, **kwargs)
     63 @wraps(plot)
     64 def _plot(*args, **kwargs) -> Ax:
     65     """New plot function with equal aspect ratio axes returned."""
---> 66     ax = plot(*args, **kwargs)
     67     ax.set_aspect("equal")
     68     return ax

File ~\AppData\Local\anaconda3\envs\asdf\lib\site-packages\tidy3d\components\viz.py:52, in add_ax_if_none.<locals>._plot(*args, **kwargs)
     50     ax = make_ax()
     51     kwargs["ax"] = ax
---> 52 return plot(*args, **kwargs)

File ~\AppData\Local\anaconda3\envs\asdf\lib\site-packages\tidy3d\components\base_sim\simulation.py:457, in AbstractSimulation.plot_structures(self, x, y, z, ax, hlim, vlim)
    430 """Plot each of simulation's structures on a plane defined by one nonzero x,y,z coordinate.
    431 
    432 Parameters
   (...)
    450     The supplied or created matplotlib axes.
    451 """
    453 hlim_new, vlim_new = Scene._get_plot_lims(
    454     bounds=self.simulation_bounds, x=x, y=y, z=z, hlim=hlim, vlim=vlim
    455 )
--> 457 return self.scene.plot_structures(x=x, y=y, z=z, ax=ax, hlim=hlim_new, vlim=vlim_new)

File ~\AppData\Local\anaconda3\envs\asdf\lib\site-packages\tidy3d\components\viz.py:66, in equal_aspect.<locals>._plot(*args, **kwargs)
     63 @wraps(plot)
     64 def _plot(*args, **kwargs) -> Ax:
     65     """New plot function with equal aspect ratio axes returned."""
---> 66     ax = plot(*args, **kwargs)
     67     ax.set_aspect("equal")
     68     return ax

File ~\AppData\Local\anaconda3\envs\asdf\lib\site-packages\tidy3d\components\viz.py:52, in add_ax_if_none.<locals>._plot(*args, **kwargs)
     50     ax = make_ax()
     51     kwargs["ax"] = ax
---> 52 return plot(*args, **kwargs)

File ~\AppData\Local\anaconda3\envs\asdf\lib\site-packages\tidy3d\components\scene.py:385, in Scene.plot_structures(self, x, y, z, ax, hlim, vlim)
    351 @equal_aspect
    352 @add_ax_if_none
    353 def plot_structures(
   (...)
    360     vlim: Tuple[float, float] = None,
    361 ) -> Ax:
    362     """Plot each of scene's structures on a plane defined by one nonzero x,y,z coordinate.
    363 
    364     Parameters
   (...)
    382         The supplied or created matplotlib axes.
    383     """
--> 385     medium_shapes = self._get_structures_2dbox(
    386         structures=self.structures, x=x, y=y, z=z, hlim=hlim, vlim=vlim
    387     )
    388     medium_map = self.medium_map
    389     for medium, shape in medium_shapes:

File ~\AppData\Local\anaconda3\envs\asdf\lib\site-packages\tidy3d\components\scene.py:535, in Scene._get_structures_2dbox(self, structures, x, y, z, hlim, vlim)
    533 medium_shapes = []
    534 for structure in structures:
--> 535     intersections = plane.intersections_with(structure.geometry)
    536     for shape in intersections:
    537         if not shape.is_empty:

File ~\AppData\Local\anaconda3\envs\asdf\lib\site-packages\tidy3d\components\geometry\base.py:1956, in Box.intersections_with(self, other)
   1954 pos = self.center[normal_ind]
   1955 xyz_kwargs = {dim: pos}
-> 1956 shapes_plane = other.intersections_plane(**xyz_kwargs)
   1958 # intersect all shapes with the input self
   1959 bs_min, bs_max = (self.pop_axis(bounds, axis=normal_ind)[1] for bounds in self.bounds)

File ~\AppData\Local\anaconda3\envs\asdf\lib\site-packages\tidy3d\components\geometry\base.py:238, in Geometry.intersections_plane(self, x, y, z)
    236     last, indices = self.pop_axis((0, 1, 2), axis)
    237     to_2D = to_2D[list(indices) + [last, 3]]
--> 238 return self.intersections_tilted_plane(normal, origin, to_2D)

File ~\AppData\Local\anaconda3\envs\asdf\lib\site-packages\tidy3d\components\geometry\base.py:2275, in Transformed.intersections_tilted_plane(self, normal, origin, to_2D)
   2254 def intersections_tilted_plane(
   2255     self, normal: Coordinate, origin: Coordinate, to_2D: MatrixReal4x4
   2256 ) -> List[Shapely]:
   2257     """Return a list of shapely geometries at the plane specified by normal and origin.
   2258 
   2259     Parameters
   (...)
   2273         `Shapely's Documentaton <[https://shapely.readthedocs.io/en/stable/project.html>`_](https://shapely.readthedocs.io/en/stable/project.html%3E%60_).
   2274     """
-> 2275     return self.geometry.intersections_tilted_plane(
   2276         tuple(np.dot((normal[0], normal[1], normal[2], 0.0), self.transform)[:3]),
   2277         tuple(np.dot(self.inverse, (origin[0], origin[1], origin[2], 1.0))[:3]),
   2278         np.dot(to_2D, self.transform),
   2279     )

File ~\AppData\Local\anaconda3\envs\asdf\lib\site-packages\tidy3d\components\geometry\base.py:55, in requires_trimesh.<locals>._fn(*args, **kwargs)
     49 if trimesh is None:
     50     raise Tidy3dImportError(
     51         "The package 'trimesh' is required for this operation, but it was not found. "
     52         "Please install the 'trimesh' dependencies using, for example, "
     53         "'pip install -r requirements/trimesh.txt'."
     54     )
---> 55 return fn(*args, **kwargs)

File ~\AppData\Local\anaconda3\envs\asdf\lib\site-packages\tidy3d\components\geometry\polyslab.py:579, in PolySlab.intersections_tilted_plane(self, normal, origin, to_2D)
    577     return []
    578 path, _ = section.to_planar(to_2D=to_2D)
--> 579 return path.polygons_full.tolist()

AttributeError: 'list' object has no attribute 'tolist'

If I uncomment the line to rotate the box rather than the poly, it gives a different error message:

17:07:45 Eastern Daylight Time WARNING:                                         
                               'geometry=Transformed(type='Transformed',        
                               geometry=Box(type='Box', center=(0.0, 0.0, 0.0), 
                               size=(2.0, inf, 1.0)), transform=array([[        
                               0.70710678, -0.70710678,  0.        ,  0.        
                               ],                                               
                                      [ 0.70710678,  0.70710678,  0.        ,   
                               0.        ],                                     
                                      [ 0.        ,  0.        ,  1.        ,   
                               0.        ],                                     
                                      [ 0.        ,  0.        ,  0.        ,   
                               1.        ]])) name=None type='Structure'        
                               medium=Medium(name=None, frequency_range=None,   
                               allow_gain=False, nonlinear_spec=None,           
                               modulation_spec=None, heat_spec=None,            
                               type='Medium', permittivity=6.0,                 
                               conductivity=0.0)' (at                           
                               `simulation.structures[0]`) is completely outside
                               of simulation domain.                            
---------------------------------------------------------------------------
ValidationError                           Traceback (most recent call last)
Cell In[21], line 36
     33 fig, (ax1,ax2) = plt.subplots(1, 2, tight_layout=True, figsize=(10, 4))
     35 # Waveguide cross-section
---> 36 ax1 = sim.plot(y=0, ax=ax1)
     38 # Top view
     39 ax2 = sim.plot(z=0, ax=ax2)

File ~\AppData\Local\anaconda3\envs\asdf\lib\site-packages\tidy3d\components\viz.py:66, in equal_aspect.<locals>._plot(*args, **kwargs)
     63 @wraps(plot)
     64 def _plot(*args, **kwargs) -> Ax:
     65     """New plot function with equal aspect ratio axes returned."""
---> 66     ax = plot(*args, **kwargs)
     67     ax.set_aspect("equal")
     68     return ax

File ~\AppData\Local\anaconda3\envs\asdf\lib\site-packages\tidy3d\components\viz.py:52, in add_ax_if_none.<locals>._plot(*args, **kwargs)
     50     ax = make_ax()
     51     kwargs["ax"] = ax
---> 52 return plot(*args, **kwargs)

File ~\AppData\Local\anaconda3\envs\asdf\lib\site-packages\tidy3d\components\simulation.py:1791, in Simulation.plot(self, x, y, z, ax, source_alpha, monitor_alpha, hlim, vlim, **patch_kwargs)
   1761 """Plot each of simulation's components on a plane defined by one nonzero x,y,z coordinate.
   1762 
   1763 Parameters
   (...)
   1785     The supplied or created matplotlib axes.
   1786 """
   1787 hlim, vlim = Scene._get_plot_lims(
   1788     bounds=self.simulation_bounds, x=x, y=y, z=z, hlim=hlim, vlim=vlim
   1789 )
-> 1791 ax = self.plot_structures(ax=ax, x=x, y=y, z=z, hlim=hlim, vlim=vlim)
   1792 ax = self.plot_sources(ax=ax, x=x, y=y, z=z, hlim=hlim, vlim=vlim, alpha=source_alpha)
   1793 ax = self.plot_monitors(ax=ax, x=x, y=y, z=z, hlim=hlim, vlim=vlim, alpha=monitor_alpha)

File ~\AppData\Local\anaconda3\envs\asdf\lib\site-packages\tidy3d\components\viz.py:66, in equal_aspect.<locals>._plot(*args, **kwargs)
     63 @wraps(plot)
     64 def _plot(*args, **kwargs) -> Ax:
     65     """New plot function with equal aspect ratio axes returned."""
---> 66     ax = plot(*args, **kwargs)
     67     ax.set_aspect("equal")
     68     return ax

File ~\AppData\Local\anaconda3\envs\asdf\lib\site-packages\tidy3d\components\viz.py:52, in add_ax_if_none.<locals>._plot(*args, **kwargs)
     50     ax = make_ax()
     51     kwargs["ax"] = ax
---> 52 return plot(*args, **kwargs)

File ~\AppData\Local\anaconda3\envs\asdf\lib\site-packages\tidy3d\components\base_sim\simulation.py:457, in AbstractSimulation.plot_structures(self, x, y, z, ax, hlim, vlim)
    430 """Plot each of simulation's structures on a plane defined by one nonzero x,y,z coordinate.
    431 
    432 Parameters
   (...)
    450     The supplied or created matplotlib axes.
    451 """
    453 hlim_new, vlim_new = Scene._get_plot_lims(
    454     bounds=self.simulation_bounds, x=x, y=y, z=z, hlim=hlim, vlim=vlim
    455 )
--> 457 return self.scene.plot_structures(x=x, y=y, z=z, ax=ax, hlim=hlim_new, vlim=vlim_new)

File ~\AppData\Local\anaconda3\envs\asdf\lib\site-packages\tidy3d\components\viz.py:66, in equal_aspect.<locals>._plot(*args, **kwargs)
     63 @wraps(plot)
     64 def _plot(*args, **kwargs) -> Ax:
     65     """New plot function with equal aspect ratio axes returned."""
---> 66     ax = plot(*args, **kwargs)
     67     ax.set_aspect("equal")
     68     return ax

File ~\AppData\Local\anaconda3\envs\asdf\lib\site-packages\tidy3d\components\viz.py:52, in add_ax_if_none.<locals>._plot(*args, **kwargs)
     50     ax = make_ax()
     51     kwargs["ax"] = ax
---> 52 return plot(*args, **kwargs)

File ~\AppData\Local\anaconda3\envs\asdf\lib\site-packages\tidy3d\components\scene.py:391, in Scene.plot_structures(self, x, y, z, ax, hlim, vlim)
    389 for medium, shape in medium_shapes:
    390     mat_index = medium_map[medium]
--> 391     ax = self._plot_shape_structure(medium=medium, mat_index=mat_index, shape=shape, ax=ax)
    392 ax = self._set_plot_bounds(bounds=self.bounds, ax=ax, x=x, y=y, z=z, hlim=hlim, vlim=vlim)
    394 # clean up the axis display

File ~\AppData\Local\anaconda3\envs\asdf\lib\site-packages\tidy3d\components\scene.py:404, in Scene._plot_shape_structure(self, medium, mat_index, shape, ax)
    402 """Plot a structure's cross section shape for a given medium."""
    403 plot_params_struct = self._get_structure_plot_params(medium=medium, mat_index=mat_index)
--> 404 ax = self.box.plot_shape(shape=shape, plot_params=plot_params_struct, ax=ax)
    405 return ax

File ~\AppData\Local\anaconda3\envs\asdf\lib\site-packages\tidy3d\components\base.py:51, in cache.<locals>.cached_property_getter(self)
     48 if stored_value is not None:
     49     return stored_value
---> 51 computed_value = prop(self)
     52 self._cached_properties[prop_name] = computed_value
     53 return computed_value

File ~\AppData\Local\anaconda3\envs\asdf\lib\site-packages\tidy3d\components\scene.py:172, in Scene.box(self)
    162 @cached_property
    163 def box(self) -> Box:
    164     """Automatically defined scene's :class:`.Box`.
    165 
    166     Returns
   (...)
    169         Scene's box.
    170     """
--> 172     return Box(center=self.center, size=self.size)

File ~\AppData\Local\anaconda3\envs\asdf\lib\site-packages\tidy3d\components\base.py:135, in Tidy3dBaseModel.__init__(self, **kwargs)
    133 """Init method, includes post-init validators."""
    134 log.begin_capture()
--> 135 super().__init__(**kwargs)
    136 self._post_init_validators()
    137 log.end_capture(self)

File ~\AppData\Local\anaconda3\envs\asdf\lib\site-packages\pydantic\v1\main.py:341, in BaseModel.__init__(__pydantic_self__, **data)
    339 values, fields_set, validation_error = validate_model(__pydantic_self__.__class__, data)
    340 if validation_error:
--> 341     raise validation_error
    342 try:
    343     object_setattr(__pydantic_self__, '__dict__', values)

ValidationError: 1 validation error for Box
size -> 2
  ensure this value is greater than or equal to 0 (type=value_error.number.not_ge; limit_value=0)

Note that without rotation, I expect the plot to look something like the below:

download

@lucas-flexcompute
Copy link
Collaborator

Thanks for the report, @samuelkim16

Geometry transformations can't yet handle inf coordinates. In this case, replacing the box y size with something like 1000, should give the expected result:
image

I'll see if/how we can implement proper handling of inf values. If not possible, we should at least give a better error report.

lucas-flexcompute added a commit that referenced this issue Mar 13, 2024
Signed-off-by: Lucas Heitzmann Gabrielli <lucas@flexcompute.com>
@samuelkim16
Copy link
Author

Thank you for the reply!

For the box, I replaced it with the following lines:

box_geom = td.Box(center=[0, 0, 0], size=[2, 100, 1])
box_geom = box_geom.rotated(axis=[0, 0, 1], angle=np.pi/4)

but I am still getting the AttributeError: 'list' object has no attribute 'tolist' error (listed in my first post) that I saw with the triangle. I am also still seeing the same error with the triangle, although it does not use td.inf.

I am on tidy3d version 2.5.2

@lucas-flexcompute
Copy link
Collaborator

Strange, I can't reproduce this error. Do you have all your dependencies updated? In particular shapely and trimesh are heavily used for transformations.

@samuelkim16
Copy link
Author

Strange. I updated shapely from 2.0.2 to 2.0.3, and trimesh from 4.0.9 to 4.2.0, and I am still getting the same error.

@lucas-flexcompute
Copy link
Collaborator

Ok, I see the problem. Your trimesh is actually too recent! We fix trimesh version at 3.20 (see the requirements file).

It seems to_plannar in trimesh 4+ is incompatible with 3+. If you downgrade trimesh to 3.20 it should work. I've tested here in a virtual environment with 4.2 and got the same error as you.

@samuelkim16
Copy link
Author

Ah perfect, that was it! Thank you for your help!

Are there any plans to upgrade to trimesh 4+? One of the other packages (gdsfactory) lists >=4 as a requirement, although I haven't played around with it yet to see if downgrading to 3.20 breaks anything in gdsfactory.

@lucas-flexcompute
Copy link
Collaborator

That was not on our roadmap, but I've just created the necessary ticket (#1544) to investigate all the places where we use it and how it can be upgraded!

lucas-flexcompute added a commit that referenced this issue Mar 18, 2024
Signed-off-by: Lucas Heitzmann Gabrielli <lucas@flexcompute.com>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants