This notebook will give a tutorial on importing and working with .stl surface mesh files in Tidy3D.
To use this functionality, remember to install Tidy3D as pip install "tidy3d[trimesh]", which will install optional dependencies needed for processing surface meshes.
We also provide a conprehensive list of other tutorials such as how to define boundary conditions, how to compute the S-matrix of a device, how to interact with tidy3d's web API, and how to define self-intersecting polygons.
If you are new to the finite-difference time-domain (FDTD) method, we highly recommend going through our FDTD101 tutorials.
# standard python imports
import numpy as np
import matplotlib.pyplot as plt
# tidy3d imports
import tidy3d as td
import tidy3d.web as web
Getting started: a simple box mesh¶
We'll start with importing an STL file representing a simple slab. We need to make sure we understand the units associated with the data in the STL file. Here, we'll assume the STL data is stored in microns. The STL files used in this tutorial can be downloaded from our documentation repo.
As a reference case, we'll also make a box with the same dimensions using the standard Tidy3D approach of making shapes: using td.Box.
The STL box has size 0.8 um x 1.3 um x 0.3 um.
# make the geometry object representing the STL solid from the STL file stored on disk
box = td.TriangleMesh.from_stl(
filename="./misc/box.stl",
scale=1, # the units are already microns as desired, but this parameter can be used to change units [default: 1]
origin=(
0,
0,
0,
), # this can be used to set a custom origin for the stl solid [default: (0, 0, 0)]
solid_index=None, # sometimes, there may be more than one solid in the file; use this to select a specific one by index
)
# define material properties of the box
medium = td.Medium(permittivity=2)
# create a structure composed of the geometry and the medium
structure = td.Structure(geometry=box, medium=medium)
# to make sure the simulation runs correctly, let's also make a reference box the usual way
box_ref = td.Box(center=(0, 0, 0), size=(0.8, 1.3, 0.3))
# make the reference structure
structure_ref = td.Structure(geometry=box_ref, medium=medium)
wavelength = 0.3
f0 = td.C_0 / wavelength / np.sqrt(medium.permittivity)
# set the domain size in x, y, and z
domain_size = 2.5
# construct simulation size array
sim_size = (domain_size, domain_size, domain_size)
# Bandwidth in Hz
fwidth = f0 / 40.0
# Gaussian source offset; the source peak is at time t = offset/fwidth
offset = 4.0
# time dependence of sources
source_time = td.GaussianPulse(freq0=f0, fwidth=fwidth, offset=offset)
# Simulation run time past the source decay (around t=2*offset/fwidth)
run_time = 40 / fwidth
Create sources and monitors¶
To study the effect of the various boundary conditions, we'll define a plane wave source and a set of frequency-domain monitors in the volume of the simulation domain.
# create a plane wave source
source = td.PlaneWave(
center=(0, 0, -1),
source_time=source_time,
size=(td.inf, td.inf, 0),
direction="+",
)
# these monitors will be used to plot fields on planes through the middle of the domain in the frequency domain
monitor_xz = td.FieldMonitor(
center=(0, 0, 0), size=(domain_size, 0, domain_size), freqs=[f0], name="xz"
)
monitor_yz = td.FieldMonitor(
center=(0, 0, 0), size=(0, domain_size, domain_size), freqs=[f0], name="yz"
)
monitor_xy = td.FieldMonitor(
center=(0, 0, 0), size=(domain_size, domain_size, 0), freqs=[f0], name="xy"
)
[16:30:33] WARNING: Default value for the field monitor monitor.py:261 'colocate' setting has changed to 'True' in Tidy3D 2.4.0. All field components will be colocated to the grid boundaries. Set to 'False' to get the raw fields on the Yee grid instead.
WARNING: Default value for the field monitor monitor.py:261 'colocate' setting has changed to 'True' in Tidy3D 2.4.0. All field components will be colocated to the grid boundaries. Set to 'False' to get the raw fields on the Yee grid instead.
WARNING: Default value for the field monitor monitor.py:261 'colocate' setting has changed to 'True' in Tidy3D 2.4.0. All field components will be colocated to the grid boundaries. Set to 'False' to get the raw fields on the Yee grid instead.
Create the simulation objects¶
We'll make two simulation objects: one for the STL box, and the other for the Tidy3D box, in order to compare the fields later on.
# STL simulation
sim = td.Simulation(
size=sim_size,
grid_spec=td.GridSpec.auto(min_steps_per_wvl=20),
sources=[source],
structures=[structure],
monitors=[monitor_xz, monitor_yz, monitor_xy],
run_time=run_time,
boundary_spec=td.BoundarySpec.all_sides(td.PML()),
)
# reference simulation
sim_ref = td.Simulation(
size=sim_size,
grid_spec=td.GridSpec.auto(min_steps_per_wvl=20),
sources=[source],
structures=[structure_ref],
monitors=[monitor_xz, monitor_yz, monitor_xy],
run_time=run_time,
boundary_spec=td.BoundarySpec.all_sides(td.PML()),
)
# plot both simulations to make sure everything is set up correctly
_, (ax1, ax2) = plt.subplots(1, 2, figsize=(9, 4))
sim.plot(y=0, ax=ax1)
sim_ref.plot(y=0, ax=ax2)
plt.show()
We immediately notice a problem: the boxes are not centered in the same way! The reason is that we have not taken into account the definition of the origin in the STL file. In our Tidy3D box, the box's center coincides with (0, 0, 0). However, in the STL file, the (min, min, min) corner of the box happens to have coordinates (-1, -0.3, 0). We need to use the origin argument of from_stl to take this into account, so that STL box is centered on the simulation's origin.
Note that information regarding the STL file units and origin should be known by the user beforehand - it should be readily available in most CAD tools used for generating and manipulating STL files.
# in the STL's coordinate system, the box's center along each dim is (local_origin[dim] + size[dim] / 2)
# so in Tidy3D's coordinate system, we need and offset that is the negative of the above
box = td.TriangleMesh.from_stl(
filename="./misc/box.stl",
origin=(-(-1 + 0.4), -(-0.3 + 0.65), -(0 + 0.15)),
)
# create the structure with the updated box
structure = td.Structure(geometry=box, medium=medium)
# update the simulation object with the new structure
sim = sim.copy(update={"structures": [structure]})
# plot both simulations again
_, (ax1, ax2) = plt.subplots(1, 2, figsize=(9, 4))
sim.plot(y=0, ax=ax1)
sim_ref.plot(y=0, ax=ax2)
plt.show()
This looks much better!
To make sure that the STL geometry is correctly parsed and processed by the solver, we can add a PermittivityMonitor to the simulation to plot the permittivity profile as seen by the solver. One could also use sim.plot_eps(), but the PermittivityMonitor will take into account the use of subpixel averaging at the edges of the solid, where applicable.
monitor_eps_xz = td.PermittivityMonitor(
center=(0, 0, 0), size=(domain_size, 0, domain_size), freqs=[f0], name="xz_eps"
)
# update the simulation objects to add in the new monitor
sim = sim.copy(update={"monitors": list(sim.monitors) + [monitor_eps_xz]})
sim_ref = sim_ref.copy(update={"monitors": list(sim_ref.monitors) + [monitor_eps_xz]})
Run Simulations¶
We can now run both simulations and make sure the results match.
# STL simulation
sim_data = web.run(sim, task_name="stl_box", path="data/stl_box.hdf5", verbose=True)
# reference simulation
sim_data_ref = web.run(sim_ref, task_name="stl_box_ref", path="data/stl_box_ref.hdf5", verbose=True)
View task using web UI at webapi.py:190 'https://tidy3d.simulation.cloud/workbench?taskId=fdve- 3248cfc0-6e6d-4fe5-a439-d2dc00099eb8v1'.
View task using web UI at webapi.py:190 'https://tidy3d.simulation.cloud/workbench?taskId=fdve- 3da0e018-2c2a-4fee-8b96-b55906af3c35v1'.
Visualize and compare¶
First, for both simulation objects, let's plot the permittivity data recorded by the PermittivityMonitor.
fig, (ax1, ax2) = plt.subplots(1, 2, tight_layout=True, figsize=(8, 3))
sim_data["xz_eps"].eps_xx.real.plot(x="x", y="z", ax=ax1, cmap="binary")
sim_data_ref["xz_eps"].eps_xx.real.plot(x="x", y="z", ax=ax2, cmap="binary")
plt.show()
First, the permittivity profiles look identical in the two cases, which reassures us that the STL geometry is equivalent to the one natively built with Tidy3D. Second, we see the effects of subpixel averaging at the edges of the box in the vertical direction. This means that the top and bottom edges of the box lie partway through an FDTD cell, so an average of the box and background permittivities is used for those cells.
Next, we can plot the frequency-domain fields for both simulations and make sure they match.
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, tight_layout=True, figsize=(12, 3))
sim_data.plot_field(field_monitor_name="xz", field_name="Ex", y=0, val="real", ax=ax1)
sim_data.plot_field(field_monitor_name="xz", field_name="Ey", y=0, val="real", ax=ax2)
sim_data.plot_field(field_monitor_name="xz", field_name="Ez", y=0, val="real", ax=ax3)
plt.show()
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, tight_layout=True, figsize=(12, 3))
sim_data_ref.plot_field(
field_monitor_name="xz", field_name="Ex", y=0, val="real", ax=ax1
)
sim_data_ref.plot_field(
field_monitor_name="xz", field_name="Ey", y=0, val="real", ax=ax2
)
sim_data_ref.plot_field(
field_monitor_name="xz", field_name="Ez", y=0, val="real", ax=ax3
)
plt.show()
As expected, the fields look identical.
Multiple STL solids and multiple files¶
Next, we'll demonstrate how multiple STL models can be imported into the same simulation. Also, we demonstrate the case when the same STL file contains more than one disjoint object.
We'll consider slightly more complicated geometries here, such as a box with a hole, and a box with a concave surface in the form of an indent. Our STL import and preprocessing functionality can handle all of these cases.
Note that:
- if the same STL file contains multiple disjoint objects stored as a single STL solid, it is assumed that those objects will all have the same material properties, because they are treated as a single
Tidy3DGeometry; - if the STL contains multiple objects stored as different STL solids, then each solid can be imported individually by index using the
solid_indexargument totd.TriangleMesh.from_stl, in which case different material properties can be assigned to different STL solids.
Here are 3D plots of the different STL objects we'll consider:
objs = [None] * 3
# first object: a box with a hole
objs[0] = td.TriangleMesh.from_stl(
filename="./misc/box_with_hole.stl",
origin=(0.1, -0.2, 0),
)
# second object: a box with an indent
objs[1] = td.TriangleMesh.from_stl(
filename="./misc/box_with_indent.stl",
origin=(0, 0, -0.7),
)
# third object: two disjoint boxes in the same STL file
objs[2] = td.TriangleMesh.from_stl(
filename="./misc/two_boxes.stl",
origin=(0.9, -0.5, 1),
)
# update the simulation with these new structures; for simplicity we assume they're all made of the same material
structures = [td.Structure(geometry=i, medium=medium) for i in objs]
sim = sim.copy(update={"structures": structures})
# we've placed the objects at three different elevations along z; let's plot and make sure they are set up correctly
_, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(12, 3))
sim.plot(z=0, ax=ax1) # STL with two boxes
sim.plot(z=-0.7, ax=ax2) # STL with a box with a hole
sim.plot(z=1, ax=ax3) # STL with a box with an indent
# let's also make sure that the permittivity profile makes sense
_, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(14, 3))
sim.plot_eps(z=0, ax=ax1) # STL with two boxes
sim.plot_eps(z=-0.7, ax=ax2) # STL with a box with a hole
sim.plot_eps(z=1, ax=ax3) # STL with a box with an indent
plt.show()
The plots and permittivity profiles all match what we expect based on the STL geometries.
More complicated shapes and mesh considerations¶
Finally, let's consider the case of a sphere intersected by a cone. When different shapes touch or overlap, we have to be more careful. For example, consider the two images below:
In the first image, the sphere and cone were unionized in the CAD software before exporting to STL. Therefore, the two shapes are stitched together so that there is no overlap, and the meshes of the sphere and the cone line up perfectly; i.e., it is one unionized object.
In the second case, the sphere and cone were not unionized, so the meshes for the sphere and the cone are essentially superimposed on one another and intersect each other; this is more clearly seen in the third figure, which is a wireframe of the second one.
For best results, we strongly recommend that all solids in a given STL are unionized prior to export, as in the first image above. All objects must be water-tight for the STL handling to work correctly, and non-unionized intersecting meshes may break the water-tightness of the surface mesh.
When water-tightness is maintained, our solver can handle non-unionized geometries such as the second and third images, but some features may not work as expected: in particular, plotting of the geometry and the client-side mode solver may not work correctly.
In this example, we'll demonstrate the handling of unionized and non-unionized geometries.
# import the unionized sphere-cone STL
obj_union = td.TriangleMesh.from_stl(
filename="./misc/icecream_unionized.stl",
origin=(0, 0, 0.4),
)
# import the non-unionized sphere-cone STL
obj_nounion = td.TriangleMesh.from_stl(
filename="./misc/icecream_nonunionized.stl",
origin=(0, 0, 0.4),
)
# make two simulation objects, one with the unionized shape and one with the non-union one
sim_union = sim.copy(
update={"structures": [td.Structure(geometry=obj_union, medium=medium)]}
)
sim_nounion = sim.copy(
update={"structures": [td.Structure(geometry=obj_nounion, medium=medium)]}
)
# plot both simulations
_, (ax1, ax2) = plt.subplots(1, 2, figsize=(9, 3))
sim_union.plot_eps(y=0, ax=ax1)
sim_nounion.plot_eps(y=0, ax=ax2)
plt.show()
Run Simulations¶
We'll run both simulations to make sure the results match.
sim_data_union = web.run(
sim_union, task_name="stl_icecream_union", path="data/stl_icecream_union.hdf5", verbose=True,
)
sim_data_nounion = web.run(
sim_nounion, task_name="stl_icecream_nounion", path="data/stl_icecream_nounion.hdf5", verbose=True,
)
View task using web UI at webapi.py:190 'https://tidy3d.simulation.cloud/workbench?taskId=fdve- 255a8559-9566-415e-b959-127bb7ec40aev1'.
View task using web UI at webapi.py:190 'https://tidy3d.simulation.cloud/workbench?taskId=fdve- 8eb8c00b-d561-4689-b2a1-a06c0cf19e28v1'.
Visualize and compare¶
We can take a look at the permittivity monitors again to make sure the permittivity profiles are properly interpreted and averaged by the solver.
fig, (ax1, ax2) = plt.subplots(1, 2, tight_layout=True, figsize=(8, 3))
sim_data_union["xz_eps"].eps_xx.real.plot(x="x", y="z", ax=ax1, cmap="binary")
sim_data_nounion["xz_eps"].eps_xx.real.plot(x="x", y="z", ax=ax2, cmap="binary")
plt.show()
The permittivity profiles look correct, and as expected, we can see the effect of subpixel averaging at the edges of the shape.
Let's plot the frequency-domain fields for both simulations and make sure they match.
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, tight_layout=True, figsize=(12, 3))
sim_data_union.plot_field(field_monitor_name="xz", field_name="Ex", val="real", ax=ax1)
sim_data_union.plot_field(field_monitor_name="yz", field_name="Ex", val="real", ax=ax2)
sim_data_union.plot_field(field_monitor_name="xy", field_name="Ex", val="real", ax=ax3)
plt.show()
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, tight_layout=True, figsize=(12, 3))
sim_data_nounion.plot_field(
field_monitor_name="xz", field_name="Ex", val="real", ax=ax1
)
sim_data_nounion.plot_field(
field_monitor_name="yz", field_name="Ex", val="real", ax=ax2
)
sim_data_nounion.plot_field(
field_monitor_name="xy", field_name="Ex", val="real", ax=ax3
)
plt.show()
The fields match extremely well for both meshes. Although the solver works correctly for both the unionized and non-unionized case, the safest approach is to ensure all touching objects are unionized prior to exporting the STL.
