Enable fixpoints in grid generation

tests/test_components/test_grid_spec.py

@@ -30,6 +30,7 @@ def test_make_coords():
         periodic=(True, False, False),
         wavelength=1.0,
         num_pml_layers=(10, 4),
+        fixpoints=(),
     )
 
 
@@ -45,6 +46,7 @@ def test_make_coords_2d():
         periodic=(True, True, False),
         wavelength=1.0,
         num_pml_layers=(10, 4),
+        fixpoints=(),
     )
 
 

tidy3d/__init__.py

@@ -8,7 +8,7 @@
 
 # geometry
 from .components.geometry.base import Box, Transformed, ClipOperation, GeometryGroup
-from .components.geometry.primitives import Sphere, Cylinder
+from .components.geometry.primitives import Sphere, Cylinder, Point
 from .components.geometry.mesh import TriangleMesh
 from .components.geometry.polyslab import PolySlab
 
@@ -179,6 +179,7 @@ def set_logging_level(level: str) -> None:
     "AutoGrid",
     "Box",
     "Sphere",
+    "Point",
     "Cylinder",
     "PolySlab",
     "GeometryGroup",

tidy3d/components/geometry/primitives.py

@@ -10,7 +10,7 @@
 import shapely
 
 from ..base import cached_property, skip_if_fields_missing
-from ..types import Axis, Bound, Coordinate, MatrixReal4x4, Shapely, Tuple
+from ..types import Axis, Bound, Coordinate, MatrixReal4x4, Shapely, Tuple, Literal
 from ...exceptions import SetupError, ValidationError
 from ...constants import MICROMETER, LARGE_NUMBER
 from ...packaging import verify_packages_import
@@ -166,6 +166,19 @@ def _surface_area(self, bounds: Bound) -> float:
         return area
 
 
+class Point(Sphere):
+    """A point: a sphere of radius = 0.
+
+    Example
+    -------
+    >>> b = Point(center=(1,2,3))
+    """
+
+    radius: Literal[0] = pydantic.Field(
+        0, title="Radius", description="Radius of geometry.", units=MICROMETER
+    )
+
+
 class Cylinder(base.Centered, base.Circular, base.Planar):
     """Cylindrical geometry with optional sidewall angle along axis
     direction. When ``sidewall_angle`` is nonzero, the shape is a

tidy3d/components/grid/grid_spec.py

@@ -14,6 +14,7 @@
 from ..source import SourceType
 from ..structure import Structure, StructureType
 from ..geometry.base import Box
+from ..geometry.primitives import Point
 from ...log import log
 from ...exceptions import SetupError
 from ...constants import MICROMETER, C_0, fp_eps
@@ -86,6 +87,7 @@ def make_coords(
         periodic: bool,
         wavelength: pd.PositiveFloat,
         num_pml_layers: Tuple[pd.NonNegativeInt, pd.NonNegativeInt],
+        fixpoints: List[Point],
     ) -> Coords1D:
         """Generate 1D coords to be used as grid boundaries, based on simulation parameters.
         Symmetry, and PML layers will be treated here.
@@ -103,6 +105,8 @@ def make_coords(
             Free-space wavelength.
         num_pml_layers : Tuple[int, int]
             number of layers in the absorber + and - direction along one dimension.
+        fixpoints : List[Point]
+            A set of points that enforce grid boundaries to pass through them.
 
         Returns
         -------
@@ -121,6 +125,7 @@ def make_coords(
             wavelength=wavelength,
             symmetry=symmetry,
             is_periodic=is_periodic,
+            fixpoints=fixpoints,
         )
 
         # incorporate symmetries
@@ -423,6 +428,7 @@ def _make_coords_initial(
         wavelength: float,
         symmetry: Symmetry,
         is_periodic: bool,
+        fixpoints: List[Point],
     ) -> Coords1D:
         """Customized 1D coords to be used as grid boundaries.
 
@@ -439,6 +445,8 @@ def _make_coords_initial(
             normal to each of the three axes.
         is_periodic : bool
             Apply periodic boundary condition or not.
+        fixpoints : List[Point]
+            A set of points that enforce grid boundaries to pass through them.
 
         Returns
         -------
@@ -468,6 +476,11 @@ def _make_coords_initial(
             self.min_steps_per_wvl,
             self.dl_min,
         )
+        # insert fixpoints
+        interval_coords, max_dl_list = self.mesher.insert_fixpoints(
+            axis, interval_coords, max_dl_list, fixpoints
+        )
+
         # Put just a single pixel if 2D-like simulation
         if interval_coords.size == 1:
             dl = wavelength / self.min_steps_per_wvl
@@ -572,6 +585,15 @@ class GridSpec(Tidy3dBaseModel):
         "uses :class:`.AutoGrid`.",
     )
 
+    fixpoints: Tuple[Point, ...] = pd.Field(
+        (),
+        title="Grid specification fixpoints",
+        description="A set of points that enforce grid boundaries to pass through them. "
+        "However, some fixpoints might be skipped if they are too close. "
+        "Note: it only takes effect when at least one of the three dimensions "
+        "uses :class:`.AutoGrid`.",
+    )
+
     @property
     def auto_grid_used(self) -> bool:
         """True if any of the three dimensions uses :class:`.AutoGrid`."""
@@ -688,6 +710,7 @@ def make_grid(
                 periodic=periodic[idim],
                 wavelength=wavelength,
                 num_pml_layers=num_pml_layers[idim],
+                fixpoints=list(self.fixpoints),
             )
 
         coords = Coords(**coords_dict)
@@ -700,6 +723,7 @@ def auto(
         min_steps_per_wvl: pd.PositiveFloat = 10.0,
         max_scale: pd.PositiveFloat = 1.4,
         override_structures: List[StructureType] = (),
+        fixpoints: List[Point] = (),
         dl_min: pd.NonNegativeFloat = 0.0,
         mesher: MesherType = GradedMesher(),
     ) -> GridSpec:
@@ -719,6 +743,8 @@ def auto(
             A list of structures that is added on top of the simulation structures in
             the process of generating the grid. This can be used to refine the grid or make it
             coarser depending than the expected need for higher/lower resolution regions.
+        fixpoints : List[Point]
+            A set of points that enforce grid boundaries to pass through them.
         dl_min: pd.NonNegativeFloat
             Lower bound of grid size.
         mesher : MesherType = GradedMesher()
@@ -742,6 +768,7 @@ def auto(
             grid_y=grid_1d,
             grid_z=grid_1d,
             override_structures=override_structures,
+            fixpoints=fixpoints,
         )
 
     @classmethod

tidy3d/components/grid/mesher.py

@@ -17,6 +17,7 @@
 from ..types import Axis, ArrayFloat1D
 from ..structure import Structure, MeshOverrideStructure, StructureType
 from ..medium import AnisotropicMedium, Medium2D, PECMedium
+from ..geometry.primitives import Point
 from ...exceptions import SetupError, ValidationError
 from ...constants import C_0, fp_eps
 from ...log import log
@@ -43,6 +44,16 @@ def parse_structures(
     ) -> Tuple[ArrayFloat1D, ArrayFloat1D]:
         """Calculate the positions of all bounding box interfaces along a given axis."""
 
+    @abstractmethod
+    def insert_fixpoints(
+        self,
+        axis: Axis,
+        interval_coords: ArrayFloat1D,
+        max_dl_list: ArrayFloat1D,
+        fixpoints: List[Point],
+    ) -> Tuple[ArrayFloat1D, ArrayFloat1D]:
+        """Insert fixpoints to the intervals."""
+
     @abstractmethod
     def make_grid_multiple_intervals(
         self,
@@ -58,6 +69,63 @@ class GradedMesher(Mesher):
     """Implements automatic nonuniform meshing with a set minimum steps per wavelength and
     a graded mesh expanding from higher- to lower-resolution regions."""
 
+    def insert_fixpoints(
+        self,
+        axis: Axis,
+        interval_coords: ArrayFloat1D,
+        max_dl_list: ArrayFloat1D,
+        fixpoints: List[Point],
+    ) -> Tuple[ArrayFloat1D, ArrayFloat1D]:
+        """Insert fixpoints to the intervals.
+
+        Parameters
+        ----------
+        axis : Axis
+            Axis index along which to operate.
+        interval_coords : ArrayFloat1D
+            Coordinate of interval boundaries.
+        max_dl_list : ArrayFloat1D
+            Maximal allowed step size of each interval generated from `parse_structures`.
+        fixpoints : List[Point]
+            A set of points that enforce grid boundaries to pass through them.
+
+        Returns
+        -------
+        interval_coords : Array
+            An array of coordinates, where the first element is the simulation min boundary, the
+            last element is the simulation max boundary, and the intermediate coordinates are all
+            locations with a fixpoint or a structure has a bounding box edge along the specified axis.
+            The boundaries are filtered such that no interval is smaller than the smallest
+            of the ``max_steps``.
+        max_steps : array_like
+            An array of size ``interval_coords.size - 1`` giving the maximum grid step required in
+            each ``interval_coords[i]:interval_coords[i+1]`` interval, depending on the materials
+            in that interval, the supplied wavelength, and the minimum required step per wavelength.
+
+        """
+        # Don't do anything for a 2D-like simulation, or no fix points
+        if interval_coords.size == 1 or len(fixpoints) < 1:
+            return interval_coords, max_dl_list
+
+        min_step = np.amin(max_dl_list) * 0.5
+        for point in fixpoints:
+            new_coord = point.center[axis]
+            # Skip if the point is outside the domain
+            if new_coord >= interval_coords[-1] or new_coord <= interval_coords[0]:
+                continue
+            # search insertion location
+            ind = np.searchsorted(interval_coords, new_coord, side="left")
+            # Skip fixpoints if the distance to the existing interval boundarires are
+            # smaller than half of the minimal maximal step size
+            if abs(new_coord - interval_coords[ind]) < min_step:
+                continue
+            if abs(new_coord - interval_coords[ind - 1]) < min_step:
+                continue
+            # insertion
+            interval_coords = np.insert(interval_coords, ind, new_coord)
+            max_dl_list = np.insert(max_dl_list, ind - 1, max_dl_list[ind - 1])
+        return interval_coords, max_dl_list
+
     def parse_structures(
         self,
         axis: Axis,

tidy3d/components/simulation.py

@@ -712,6 +712,12 @@ def plot_grid(
             )
             ax.add_patch(rect)
 
+        # Plot fixpoints
+        for point in self.grid_spec.fixpoints:
+            _, (x_point, y_point) = point.pop_axis(point.center, axis=axis)
+            x_point, y_point = (self._evaluate_inf(v) for v in (x_point, y_point))
+            ax.scatter(x_point, y_point, color=plot_params.edgecolor)
+
         ax = Scene._set_plot_bounds(
             bounds=self.simulation_bounds, ax=ax, x=x, y=y, z=z, hlim=hlim, vlim=vlim
         )