Learning FDTD through Examples / How to define...

How to define arbitrary 3D geometries using STL file import in Tidy3D FDTD¶

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.

[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"
)

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()

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()

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)

[15:00:20] Created task 'stl_box' with task_id                                  webapi.py:139
           'fdve-08a21e81-6a03-4f08-ad26-ba713e43fdbcv1'.

Output()

Output()

[15:00:25] status = queued                                                      webapi.py:269

Output()

[15:00:37] status = preprocess                                                  webapi.py:263

[15:00:41] Maximum FlexCredit cost: 0.101. Use 'web.real_cost(task_id)' to get  webapi.py:286
           the billed FlexCredit cost after a simulation run.

           starting up solver                                                   webapi.py:290

           running solver                                                       webapi.py:300

Output()

[15:00:52] early shutoff detected, exiting.                                     webapi.py:313

           status = postprocess                                                 webapi.py:330

Output()

[15:00:56] status = success                                                     webapi.py:337

Output()

[15:01:00] loading SimulationData from data/stl_box.hdf5                        webapi.py:509

[15:01:00] Created task 'stl_box_ref' with task_id                              webapi.py:139
           'fdve-1af9a1ee-a01d-4dd0-a333-60953c584be7v1'.

Output()

[15:01:02] status = queued                                                      webapi.py:269

Output()

[15:01:06] status = preprocess                                                  webapi.py:263

[15:01:11] Maximum FlexCredit cost: 0.101. Use 'web.real_cost(task_id)' to get  webapi.py:286
           the billed FlexCredit cost after a simulation run.

           starting up solver                                                   webapi.py:290

           running solver                                                       webapi.py:300

Output()

[15:01:20] early shutoff detected, exiting.                                     webapi.py:313

           status = postprocess                                                 webapi.py:330

Output()

[15:01:26] status = success                                                     webapi.py:337

Output()

[15:01:32] loading SimulationData from data/stl_box_ref.hdf5                    webapi.py:509

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()

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()

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:

[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()

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: