Importing STL files#

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.

[1]:
# 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.

[2]:
# 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.

[3]:
# 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.

[4]:
# 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()

../_images/notebooks_STLImport_7_0.png

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.

[5]:
# 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()

../_images/notebooks_STLImport_9_0.png

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.

[6]:
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.

[7]:
# 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)

[16:30:34] Created task 'stl_box' with task_id                     webapi.py:188
           'fdve-3248cfc0-6e6d-4fe5-a439-d2dc00099eb8v1'.                       
[16:30:36] status = queued                                         webapi.py:361
[16:30:45] status = preprocess                                     webapi.py:355
[16:30:49] Maximum FlexCredit cost: 0.101. Use                     webapi.py:341
           'web.real_cost(task_id)' to get the billed FlexCredit                
           cost after a simulation run.                                         
           starting up solver                                      webapi.py:377
           running solver                                          webapi.py:386
           To cancel the simulation, use 'web.abort(task_id)' or   webapi.py:387
           'web.delete(task_id)' or abort/delete the task in the                
           web UI. Terminating the Python script will not stop the              
           job running on the cloud.                                            
[16:30:58] early shutoff detected, exiting.                        webapi.py:404
[16:30:59] status = postprocess                                    webapi.py:419
[16:31:03] status = success                                        webapi.py:426
[16:31:05] loading SimulationData from data/stl_box.hdf5           webapi.py:590
           Created task 'stl_box_ref' with task_id                 webapi.py:188
           'fdve-3da0e018-2c2a-4fee-8b96-b55906af3c35v1'.                       
[16:31:06] status = queued                                         webapi.py:361
[16:31:15] status = preprocess                                     webapi.py:355
[16:31:19] Maximum FlexCredit cost: 0.101. Use                     webapi.py:341
           'web.real_cost(task_id)' to get the billed FlexCredit                
           cost after a simulation run.                                         
           starting up solver                                      webapi.py:377
           running solver                                          webapi.py:386
           To cancel the simulation, use 'web.abort(task_id)' or   webapi.py:387
           'web.delete(task_id)' or abort/delete the task in the                
           web UI. Terminating the Python script will not stop the              
           job running on the cloud.                                            
[16:31:27] status = postprocess                                    webapi.py:419
[16:31:32] status = success                                        webapi.py:426
[16:31:34] loading SimulationData from data/stl_box_ref.hdf5       webapi.py:590

Visualize and compare#

First, for both simulation objects, let’s plot the permittivity data recorded by the PermittivityMonitor.

[8]:
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()

../_images/notebooks_STLImport_15_0.png

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.

[9]:
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()

../_images/notebooks_STLImport_17_0.png
../_images/notebooks_STLImport_17_1.png

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 Tidy3D Geometry; - if the STL contains multiple objects stored as different STL solids, then each solid can be imported individually by index using the solid_index argument to td.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:

3D plot of the first object 3D plot of the second object 3D plot of the third object

[10]:
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()

../_images/notebooks_STLImport_20_0.png
../_images/notebooks_STLImport_20_1.png

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:

Unionized object Nonunionized object Wire frame rendering

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.

[11]:
# 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()

../_images/notebooks_STLImport_23_0.png

Run Simulations#

We’ll run both simulations to make sure the results match.

[12]:
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,
)

[16:31:40] Created task 'stl_icecream_union' with task_id          webapi.py:188
           'fdve-255a8559-9566-415e-b959-127bb7ec40aev1'.                       
[16:31:41] status = queued                                         webapi.py:361
[16:31:51] status = preprocess                                     webapi.py:355
[16:31:55] Maximum FlexCredit cost: 0.118. Use                     webapi.py:341
           'web.real_cost(task_id)' to get the billed FlexCredit                
           cost after a simulation run.                                         
           starting up solver                                      webapi.py:377
           running solver                                          webapi.py:386
           To cancel the simulation, use 'web.abort(task_id)' or   webapi.py:387
           'web.delete(task_id)' or abort/delete the task in the                
           web UI. Terminating the Python script will not stop the              
           job running on the cloud.                                            
[16:32:07] early shutoff detected, exiting.                        webapi.py:404
           status = postprocess                                    webapi.py:419
[16:32:13] status = success                                        webapi.py:426
[16:32:15] loading SimulationData from                             webapi.py:590
           data/stl_icecream_union.hdf5                                         
           Created task 'stl_icecream_nounion' with task_id        webapi.py:188
           'fdve-8eb8c00b-d561-4689-b2a1-a06c0cf19e28v1'.                       
[16:32:16] status = queued                                         webapi.py:361
[16:32:25] status = preprocess                                     webapi.py:355
[16:32:29] Maximum FlexCredit cost: 0.118. Use                     webapi.py:341
           'web.real_cost(task_id)' to get the billed FlexCredit                
           cost after a simulation run.                                         
           starting up solver                                      webapi.py:377
           running solver                                          webapi.py:386
           To cancel the simulation, use 'web.abort(task_id)' or   webapi.py:387
           'web.delete(task_id)' or abort/delete the task in the                
           web UI. Terminating the Python script will not stop the              
           job running on the cloud.                                            
[16:32:43] early shutoff detected, exiting.                        webapi.py:404
           status = postprocess                                    webapi.py:419
[16:32:47] status = success                                        webapi.py:426
[16:32:49] loading SimulationData from                             webapi.py:590
           data/stl_icecream_nounion.hdf5                                       

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.

[13]:
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()

../_images/notebooks_STLImport_27_0.png

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.

[14]:
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()

../_images/notebooks_STLImport_29_0.png
../_images/notebooks_STLImport_29_1.png

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.

[ ]: