diff --git a/tests/test_components/test_layerrefinement.py b/tests/test_components/test_layerrefinement.py
index 59efc7551f..8716f17f33 100644
--- a/tests/test_components/test_layerrefinement.py
+++ b/tests/test_components/test_layerrefinement.py
@@ -104,19 +104,19 @@ def test_layerrefinement():
         assert layer.size[axis] == 1
         assert layer.size[(axis + 1) % 3] == td.inf
         assert layer.size[(axis + 2) % 3] == td.inf
-        assert not layer._is_inplane_bounded
+        assert not layer._is_inplane_bounded(layer)
 
     layer = LayerRefinementSpec.from_bounds(axis=axis, rmin=(0, 0, 0), rmax=(1, 2, 3))
     layer = LayerRefinementSpec.from_bounds(rmin=(0, 0, 0), rmax=(1, 2, 3))
     assert layer.axis == 0
     assert np.isclose(layer.length_axis, 1)
     assert np.isclose(layer.center_axis, 0.5)
-    assert layer._is_inplane_bounded
+    assert layer._is_inplane_bounded(layer)
 
     # from structures
     structures = [td.Structure(geometry=td.Box(size=(td.inf, 2, 3)), medium=td.Medium())]
     layer = LayerRefinementSpec.from_structures(structures)
-    assert layer._is_inplane_bounded
+    assert layer._is_inplane_bounded(layer)
     assert layer.axis == 1
 
     with pytest.raises(pydantic.ValidationError):
@@ -138,11 +138,11 @@ def test_layerrefinement():
 def test_layerrefinement_inplane_inside():
     # inplane inside
     layer = LayerRefinementSpec.from_layer_bounds(axis=2, bounds=(0, 1))
-    assert not layer._is_inplane_bounded
-    assert layer._inplane_inside([3e3, 4e4])
+    assert not layer._is_inplane_bounded(layer)
+    assert layer._inplane_inside(layer, [3e3, 4e4])
     layer = LayerRefinementSpec(axis=1, size=(1, 0, 1))
-    assert layer._inplane_inside([0, 0])
-    assert not layer._inplane_inside([2, 0])
+    assert layer._inplane_inside(layer, [0, 0])
+    assert not layer._inplane_inside(layer, [2, 0])
 
 
 def test_layerrefinement_snapping_points():
@@ -263,7 +263,11 @@ def count_grids_within_layer(sim_t):
     assert (
         len(
             sim2.grid_spec.all_override_structures(
-                list(sim2.structures), 1.0, sim2.size, lumped_elements
+                list(sim2.structures),
+                1.0,
+                lumped_elements,
+                sim2._internal_layerrefinement_boundary_types,
+                sim2.bounds,
             )
         )
         == 1
@@ -281,7 +285,11 @@ def count_grids_within_layer(sim_t):
     assert (
         len(
             sim2.grid_spec.all_override_structures(
-                list(sim2.structures), 1.0, sim2.size, lumped_elements
+                list(sim2.structures),
+                1.0,
+                lumped_elements,
+                sim2._internal_layerrefinement_boundary_types,
+                sim2.bounds,
             )
         )
         == 2
@@ -425,7 +433,11 @@ def test_dl_min_from_smallest_feature():
         ),
         medium=td.PECMedium(),
     )
-
+    boundary_types = [[None] * 2] * 3
+    sim_bounds = [
+        [-td.inf] * 3,
+        [td.inf] * 3,
+    ]
     # check expected dl_min
     layer_spec = td.LayerRefinementSpec(
         axis=2,
@@ -434,7 +446,7 @@ def test_dl_min_from_smallest_feature():
             convex_resolution=10,
         ),
     )
-    dl_min = layer_spec._dl_min_from_smallest_feature([structure])
+    dl_min = layer_spec._dl_min_from_smallest_feature([structure], sim_bounds, boundary_types)
     assert np.allclose(0.3 / 10, dl_min)
 
     layer_spec = td.LayerRefinementSpec(
@@ -442,7 +454,7 @@ def test_dl_min_from_smallest_feature():
         size=(td.inf, td.inf, 2),
         corner_finder=td.CornerFinderSpec(mixed_resolution=10),
     )
-    dl_min = layer_spec._dl_min_from_smallest_feature([structure])
+    dl_min = layer_spec._dl_min_from_smallest_feature([structure], sim_bounds, boundary_types)
     assert np.allclose(0.2 / 10, dl_min)
 
     layer_spec = td.LayerRefinementSpec(
@@ -452,7 +464,7 @@ def test_dl_min_from_smallest_feature():
             concave_resolution=10,
         ),
     )
-    dl_min = layer_spec._dl_min_from_smallest_feature([structure])
+    dl_min = layer_spec._dl_min_from_smallest_feature([structure], sim_bounds, boundary_types)
     assert np.allclose(0.1 / 10, dl_min)
 
     # check grid is generated succesfully
@@ -517,7 +529,8 @@ def test_gap_meshing():
             layer_refinement_specs=[
                 td.LayerRefinementSpec(
                     axis=2,
-                    corner_finder=None,
+                    corner_snapping=False,
+                    corner_refinement=None,
                     gap_meshing_iters=num_iters,
                     size=[td.inf, td.inf, 2],
                     center=[0, 0, 1],
@@ -583,7 +596,8 @@ def test_gap_meshing():
                 td.LayerRefinementSpec(
                     axis=1,
                     size=(td.inf, 0.2, td.inf),
-                    corner_finder=None,
+                    corner_snapping=False,
+                    corner_refinement=None,
                     gap_meshing_iters=1,
                     dl_min_from_gap_width=True,
                 )
@@ -626,7 +640,8 @@ def test_gap_meshing():
                 td.LayerRefinementSpec(
                     axis=1,
                     size=(td.inf, 0.2, td.inf),
-                    corner_finder=None,
+                    corner_snapping=False,
+                    corner_refinement=None,
                     gap_meshing_iters=2,
                     dl_min_from_gap_width=True,
                     interior_disjoint_geometries=False,
@@ -658,7 +673,8 @@ def test_gap_meshing():
                 td.LayerRefinementSpec(
                     axis=1,
                     size=(td.inf, 0.2, td.inf),
-                    corner_finder=None,
+                    corner_snapping=False,
+                    corner_refinement=None,
                     gap_meshing_iters=1,
                     dl_min_from_gap_width=True,
                 )
@@ -697,7 +713,8 @@ def test_gap_meshing():
                     axis=1,
                     size=(0.5, 0.2, 0.5),
                     center=(0.5, 0, -0.5),
-                    corner_finder=None,
+                    corner_snapping=False,
+                    corner_refinement=None,
                     gap_meshing_iters=1,
                     dl_min_from_gap_width=True,
                 )
@@ -731,7 +748,8 @@ def test_gap_meshing():
                     axis=1,
                     size=(0.05, 0.2, 0.05),
                     center=(0.05, 0, 0.05),
-                    corner_finder=None,
+                    corner_snapping=False,
+                    corner_refinement=None,
                     gap_meshing_iters=2,
                     dl_min_from_gap_width=True,
                 )
@@ -765,7 +783,8 @@ def test_gap_meshing():
                     axis=1,
                     size=(0.01, 0.2, 0.01),
                     center=(10.0, 0, 10.0),
-                    corner_finder=None,
+                    corner_snapping=False,
+                    corner_refinement=None,
                     gap_meshing_iters=1,
                     dl_min_from_gap_width=True,
                 )
@@ -804,7 +823,8 @@ def test_gap_meshing():
                 td.LayerRefinementSpec(
                     axis=1,
                     size=(td.inf, 0.2, td.inf),
-                    corner_finder=None,
+                    corner_snapping=False,
+                    corner_refinement=None,
                     gap_meshing_iters=1,
                     dl_min_from_gap_width=True,
                 )
@@ -842,7 +862,8 @@ def test_gap_meshing():
                 td.LayerRefinementSpec(
                     axis=0,
                     size=(0.2, td.inf, td.inf),
-                    corner_finder=None,
+                    corner_snapping=False,
+                    corner_refinement=None,
                     gap_meshing_iters=1,
                     dl_min_from_gap_width=True,
                 )
diff --git a/tidy3d/components/geometry/base.py b/tidy3d/components/geometry/base.py
index 441a6bb7ee..41682ae736 100644
--- a/tidy3d/components/geometry/base.py
+++ b/tidy3d/components/geometry/base.py
@@ -31,6 +31,7 @@
 from tidy3d.components.autograd.derivative_utils import DerivativeInfo
 from tidy3d.components.base import Tidy3dBaseModel, cached_property
 from tidy3d.components.geometry.bound_ops import bounds_intersection, bounds_union
+from tidy3d.components.geometry.float_utils import increment_float
 from tidy3d.components.transformation import ReflectionFromPlane, RotationAroundAxis
 from tidy3d.components.types import (
     ArrayFloat2D,
@@ -2173,6 +2174,11 @@ def intersections_with(self, other: Shapely) -> list[Shapely]:
         shapely_box = Geometry.evaluate_inf_shape(shapely_box)
         return [Geometry.evaluate_inf_shape(shape) & shapely_box for shape in shapes_plane]
 
+    def slightly_enlarged_copy(self) -> Box:
+        """Box size slightly enlarged around machine precision."""
+        size = [increment_float(orig_length, 1) for orig_length in self.size]
+        return self.updated_copy(size=size)
+
     def padded_copy(
         self,
         x: Optional[tuple[pydantic.NonNegativeFloat, pydantic.NonNegativeFloat]] = None,
diff --git a/tidy3d/components/geometry/float_utils.py b/tidy3d/components/geometry/float_utils.py
new file mode 100644
index 0000000000..5ab7b438be
--- /dev/null
+++ b/tidy3d/components/geometry/float_utils.py
@@ -0,0 +1,31 @@
+"""Utilities for float manipulation."""
+
+from __future__ import annotations
+
+import numpy as np
+
+from tidy3d.constants import inf
+
+
+def increment_float(val: float, sign: int) -> float:
+    """Applies a small positive or negative shift as though `val` is a 32bit float
+    using numpy.nextafter, but additionally handles some corner cases.
+    """
+    # Infinity is left unchanged
+    if val == inf or val == -inf:
+        return val
+
+    if sign >= 0:
+        sign = 1
+    else:
+        sign = -1
+
+    # Avoid small increments within subnormal values
+    if np.abs(val) <= np.finfo(np.float32).tiny:
+        return val + sign * np.finfo(np.float32).tiny
+
+    # Numpy seems to skip over the increment from -0.0 and +0.0
+    # which is different from c++
+    val_inc = np.nextafter(val, sign * inf, dtype=np.float32)
+
+    return np.float32(val_inc)
diff --git a/tidy3d/components/geometry/utils_2d.py b/tidy3d/components/geometry/utils_2d.py
index 28c12c1380..f7399d775a 100644
--- a/tidy3d/components/geometry/utils_2d.py
+++ b/tidy3d/components/geometry/utils_2d.py
@@ -8,36 +8,13 @@
 import shapely
 
 from tidy3d.components.geometry.base import Box, ClipOperation, Geometry, GeometryGroup
+from tidy3d.components.geometry.float_utils import increment_float
 from tidy3d.components.geometry.polyslab import _MIN_POLYGON_AREA, PolySlab
 from tidy3d.components.grid.grid import Grid
 from tidy3d.components.scene import Scene
 from tidy3d.components.structure import Structure
 from tidy3d.components.types import Axis, Shapely
-from tidy3d.constants import fp_eps, inf
-
-
-def increment_float(val: float, sign: int) -> float:
-    """Applies a small positive or negative shift as though `val` is a 32bit float
-    using numpy.nextafter, but additionally handles some corner cases.
-    """
-    # Infinity is left unchanged
-    if val == inf or val == -inf:
-        return val
-
-    if sign >= 0:
-        sign = 1
-    else:
-        sign = -1
-
-    # Avoid small increments within subnormal values
-    if np.abs(val) <= np.finfo(np.float32).tiny:
-        return val + sign * np.finfo(np.float32).tiny
-
-    # Numpy seems to skip over the increment from -0.0 and +0.0
-    # which is different from c++
-    val_inc = np.nextafter(val, sign * inf, dtype=np.float32)
-
-    return np.float32(val_inc)
+from tidy3d.constants import fp_eps
 
 
 def snap_coordinate_to_grid(grid: Grid, center: float, axis: Axis) -> float:
diff --git a/tidy3d/components/grid/corner_finder.py b/tidy3d/components/grid/corner_finder.py
index cfb613d63e..59386c5595 100644
--- a/tidy3d/components/grid/corner_finder.py
+++ b/tidy3d/components/grid/corner_finder.py
@@ -86,6 +86,7 @@ def _merged_pec_on_plane(
         center: tuple[float, float] = [0, 0, 0],
         size: tuple[float, float, float] = [inf, inf, inf],
         interior_disjoint_geometries: bool = False,
+        keep_metal_only: bool = False,
     ) -> list[tuple[Any, Shapely]]:
         """On a 2D plane specified by axis = `normal_axis` and coordinate `coord`, merge geometries made of PEC.
 
@@ -95,19 +96,20 @@ def _merged_pec_on_plane(
             Axis normal to the 2D plane.
         coord : float
             Position of plane along the normal axis.
-        structure_list : List[Structure]
-            List of structures present in simulation.
-        center : Tuple[float, float] = [0, 0, 0]
+        structure_list : list[Structure]
+            list of structures present in simulation.
+        center : tuple[float, float] = [0, 0, 0]
             Center of the 2D plane (coordinate along ``axis`` is ignored)
-        size : Tuple[float, float, float] = [inf, inf, inf]
+        size : tuple[float, float, float] = [inf, inf, inf]
             Size of the 2D plane (size along ``axis`` is ignored)
         interior_disjoint_geometries: bool = False
             If ``True``, geometries on the plane must not be overlapping.
-
+        keep_metal_only: bool = False
+            If ``True``, drop all other structures that are not made of metal.
         Returns
         -------
-        List[Tuple[Any, Shapely]]
-            List of shapes and their property value on the plane after merging.
+        list[tuple[Any, Shapely]]
+            list of shapes and their property value on the plane after merging.
         """
 
         # Construct plane
@@ -125,6 +127,9 @@ def _merged_pec_on_plane(
         medium_list = [
             PEC if (mat.is_pec or isinstance(mat, LossyMetalMedium)) else mat for mat in medium_list
         ]
+        if keep_metal_only:
+            geometry_list = [geo for geo, mat in zip(geometry_list, medium_list) if mat.is_pec]
+            medium_list = [PEC for _ in geometry_list]
         # merge geometries
         merged_geos = merging_geometries_on_plane(
             geometry_list, plane, medium_list, interior_disjoint_geometries
@@ -134,11 +139,8 @@ def _merged_pec_on_plane(
 
     def _corners_and_convexity(
         self,
-        normal_axis: Axis,
-        coord: float,
-        structure_list: list[Structure],
+        merged_geos: list[tuple[Any, Shapely]],
         ravel: bool,
-        interior_disjoint_geometries: bool = False,
     ) -> tuple[ArrayFloat2D, ArrayFloat1D]:
         """On a 2D plane specified by axis = `normal_axis` and coordinate `coord`, find out corners of merged
         geometries made of PEC.
@@ -146,31 +148,16 @@ def _corners_and_convexity(
 
         Parameters
         ----------
-        normal_axis : Axis
-            Axis normal to the 2D plane.
-        coord : float
-            Position of plane along the normal axis.
-        structure_list : List[Structure]
-            List of structures present in simulation.
+        merged_geos : list[tuple[Any, Shapely]]
+            list of shapes and their property value on the plane after merging.
         ravel : bool
             Whether to put the resulting corners in a single list or per polygon.
-        interior_disjoint_geometries: bool = False
-            If ``True``, geometries made of different materials on the plane must not be overlapping.
 
         Returns
         -------
-        Tuple[ArrayFloat2D, ArrayFloat1D]
+        tuple[ArrayFloat2D, ArrayFloat1D]
             Corner coordinates and their convexity.
         """
-
-        # merge geometries
-        merged_geos = self._merged_pec_on_plane(
-            normal_axis=normal_axis,
-            coord=coord,
-            structure_list=structure_list,
-            interior_disjoint_geometries=interior_disjoint_geometries,
-        )
-
         # corner finder
         corner_list = []
         convexity_list = []
@@ -209,7 +196,10 @@ def corners(
         normal_axis: Axis,
         coord: float,
         structure_list: list[Structure],
+        center: tuple[float, float, float] = [0, 0, 0],
+        size: tuple[float, float, float] = [inf, inf, inf],
         interior_disjoint_geometries: bool = False,
+        keep_metal_only: bool = False,
     ) -> ArrayFloat2D:
         """On a 2D plane specified by axis = `normal_axis` and coordinate `coord`, find out corners of merged
         geometries made of `medium`.
@@ -221,22 +211,33 @@ def corners(
             Axis normal to the 2D plane.
         coord : float
             Position of plane along the normal axis.
-        structure_list : List[Structure]
-            List of structures present in simulation.
+        structure_list : list[Structure]
+            list of structures present in simulation.
         interior_disjoint_geometries: bool = False
             If ``True``, geometries made of different materials on the plane must not be overlapping.
+        center : tuple[float, float, float]=[0, 0, 0]
+            Center of the 2D plane (coordinate along ``axis`` is ignored).
+        size : tuple[float, float, float]=[inf, inf, inf]
+            Size of the 2D plane (size along ``axis`` is ignored).
+        keep_metal_only: bool = False
+            If ``True``, drop all other structures that are not made of metal.
         Returns
         -------
         ArrayFloat2D
             Corner coordinates.
         """
-
-        corner_list, _ = self._corners_and_convexity(
+        merged_geos = self._merged_pec_on_plane(
             normal_axis=normal_axis,
             coord=coord,
             structure_list=structure_list,
-            ravel=True,
+            center=center,
+            size=size,
             interior_disjoint_geometries=interior_disjoint_geometries,
+            keep_metal_only=keep_metal_only,
+        )
+        corner_list, _ = self._corners_and_convexity(
+            merged_geos=merged_geos,
+            ravel=True,
         )
         return corner_list
 
diff --git a/tidy3d/components/grid/grid_spec.py b/tidy3d/components/grid/grid_spec.py
index 6ad1013af9..bcffe6b2f3 100644
--- a/tidy3d/components/grid/grid_spec.py
+++ b/tidy3d/components/grid/grid_spec.py
@@ -10,7 +10,6 @@
 
 from tidy3d.components.base import Tidy3dBaseModel, cached_property, skip_if_fields_missing
 from tidy3d.components.geometry.base import Box, ClipOperation
-from tidy3d.components.geometry.utils_2d import increment_float
 from tidy3d.components.lumped_element import LumpedElementType
 from tidy3d.components.source.utils import SourceType
 from tidy3d.components.structure import MeshOverrideStructure, Structure, StructureType
@@ -22,11 +21,12 @@
     Coordinate,
     CoordinateOptional,
     PriorityMode,
+    Shapely,
     Symmetry,
     Undefined,
     annotate_type,
 )
-from tidy3d.constants import C_0, MICROMETER, dp_eps, fp_eps, inf
+from tidy3d.constants import C_0, MICROMETER, fp_eps, inf
 from tidy3d.exceptions import SetupError
 from tidy3d.log import log
 
@@ -57,6 +57,8 @@ def make_coords(
         wavelength: pd.PositiveFloat,
         num_pml_layers: tuple[pd.NonNegativeInt, pd.NonNegativeInt],
         snapping_points: tuple[CoordinateOptional, ...],
+        parse_structures_interval_coords: np.ndarray = None,
+        parse_structures_max_dl_list: np.ndarray = None,
     ) -> Coords1D:
         """Generate 1D coords to be used as grid boundaries, based on simulation parameters.
         Symmetry, and PML layers will be treated here.
@@ -79,6 +81,10 @@ def make_coords(
             number of layers in the absorber + and - direction along one dimension.
         snapping_points : Tuple[CoordinateOptional, ...]
             A set of points that enforce grid boundaries to pass through them.
+        parse_structures_interval_coords : np.ndarray, optional
+            If not None, pre-computed interval coordinates from parsing structures.
+        parse_structures_max_dl_list : np.ndarray, optional
+            If not None, pre-computed maximum grid spacing list from parsing structures.
 
         Returns
         -------
@@ -98,6 +104,8 @@ def make_coords(
             symmetry=symmetry,
             is_periodic=is_periodic,
             snapping_points=snapping_points,
+            parse_structures_interval_coords=parse_structures_interval_coords,
+            parse_structures_max_dl_list=parse_structures_max_dl_list,
         )
 
         # incorporate symmetries
@@ -589,39 +597,32 @@ def _filtered_dl(self, dl: float, sim_size: tuple[float, 3]) -> float:
         """Grid step size after applying minimal and maximal filtering."""
         return max(min(dl, self._dl_max(sim_size)), self._dl_min)
 
-    def _make_coords_initial(
+    def _compute_symmetry_domain(
         self,
-        axis: Axis,
         structures: list[StructureType],
-        wavelength: float,
         symmetry: Symmetry,
-        is_periodic: bool,
-        snapping_points: tuple[CoordinateOptional, ...],
-    ) -> Coords1D:
-        """Customized 1D coords to be used as grid boundaries.
+    ) -> tuple[Box, list, list, float]:
+        """Compute the symmetry domain from structures and symmetry parameters.
 
         Parameters
         ----------
-        axis : Axis
-            Axis of this direction.
         structures : List[StructureType]
             List of structures present in simulation.
-        wavelength : float
-            Free-space wavelength.
         symmetry : Tuple[Symmetry, Symmetry, Symmetry]
             Reflection symmetry across a plane bisecting the simulation domain
             normal to each of the three axes.
-        is_periodic : bool
-            Apply periodic boundary condition or not.
-        snapping_points : Tuple[CoordinateOptional, ...]
-            A set of points that enforce grid boundaries to pass through them.
 
         Returns
         -------
-        :class:`.Coords1D`:
-            1D coords to be used as grid boundaries.
+        symmetry_domain : Box
+            Simulation domain with symmetry applied.
+        sim_cent : list
+            Center of the simulation domain (with symmetry applied).
+        sim_size : list
+            Size of the simulation domain (with symmetry applied).
+        dl_max : float
+            Upper bound of grid step size based on total sim_size.
         """
-
         sim_cent = list(structures[0].geometry.center)
         sim_size = list(structures[0].geometry.size)
 
@@ -634,6 +635,44 @@ def _make_coords_initial(
                 sim_size[dim] /= 2
         symmetry_domain = Box(center=sim_cent, size=sim_size)
 
+        return symmetry_domain, sim_cent, sim_size, dl_max
+
+    def _parse_structures(
+        self,
+        axis: Axis,
+        structures: list[StructureType],
+        wavelength: float,
+        symmetry: Symmetry,
+        snapping_points: tuple[CoordinateOptional, ...],
+    ) -> tuple[np.ndarray, np.ndarray, Box, list]:
+        """Compute interval coordinates and max dl list from mesher.parse_structures
+
+        Parameters
+        ----------
+        axis : Axis
+            Axis of this direction.
+        structures : List[StructureType]
+            List of structures present in simulation.
+        wavelength : float
+            Free-space wavelength.
+        symmetry : Tuple[Symmetry, Symmetry, Symmetry]
+            Reflection symmetry across a plane bisecting the simulation domain
+            normal to each of the three axes.
+        snapping_points : Tuple[CoordinateOptional, ...]
+            A set of points that enforce grid boundaries to pass through them.
+
+        Returns
+        -------
+        interval_coords : np.ndarray
+            Interval coordinates for meshing.
+        max_dl_list : np.ndarray
+            Maximum grid spacing list for each interval.
+        """
+
+        symmetry_domain, sim_cent, sim_size, dl_max = self._compute_symmetry_domain(
+            structures, symmetry
+        )
+
         # New list of structures with symmetry applied
         struct_list = [Structure(geometry=symmetry_domain, medium=structures[0].medium)]
         rmin_domain, rmax_domain = symmetry_domain.bounds
@@ -657,9 +696,64 @@ def _make_coords_initial(
             self._dl_min,
             dl_max,
         )
+
+        return interval_coords, max_dl_list
+
+    def _make_coords_initial(
+        self,
+        axis: Axis,
+        structures: list[StructureType],
+        wavelength: float,
+        symmetry: Symmetry,
+        is_periodic: bool,
+        snapping_points: tuple[CoordinateOptional, ...],
+        parse_structures_interval_coords: np.ndarray = None,
+        parse_structures_max_dl_list: np.ndarray = None,
+    ) -> Coords1D:
+        """Customized 1D coords to be used as grid boundaries.
+
+        Parameters
+        ----------
+        axis : Axis
+            Axis of this direction.
+        structures : List[StructureType]
+            List of structures present in simulation.
+        wavelength : float
+            Free-space wavelength.
+        symmetry : Tuple[Symmetry, Symmetry, Symmetry]
+            Reflection symmetry across a plane bisecting the simulation domain
+            normal to each of the three axes.
+        is_periodic : bool
+            Apply periodic boundary condition or not.
+        snapping_points : Tuple[CoordinateOptional, ...]
+            A set of points that enforce grid boundaries to pass through them.
+        parse_structures_interval_coords : np.ndarray, optional
+            If not None, pre-computed interval coordinates.
+        parse_structures_max_dl_list : np.ndarray, optional
+            If not None, pre-computed maximum grid spacing list.
+
+        Returns
+        -------
+        :class:`.Coords1D`:
+            1D coords to be used as grid boundaries.
+        """
+        symmetry_domain, sim_cent, sim_size, dl_max = self._compute_symmetry_domain(
+            structures, symmetry
+        )
+
+        # Compute interval_coords and max_dl_list if not provided
+        if parse_structures_interval_coords is None or parse_structures_max_dl_list is None:
+            parse_structures_interval_coords, parse_structures_max_dl_list = self._parse_structures(
+                axis, structures, wavelength, symmetry, snapping_points
+            )
+
         # insert snapping_points
         interval_coords, max_dl_list = self.mesher.insert_snapping_points(
-            self._dl_min, axis, interval_coords, max_dl_list, snapping_points
+            self._dl_min,
+            axis,
+            parse_structures_interval_coords,
+            parse_structures_max_dl_list,
+            snapping_points,
         )
 
         # Put just a single pixel if 2D-like simulation
@@ -931,7 +1025,7 @@ def _grid_size(self, grid_size_in_vacuum: float) -> float:
 
         Parameters
         ----------
-        grid_size_in_vaccum : float
+        grid_size_in_vacuum : float
             Grid step size in vaccum.
 
         Returns
@@ -957,7 +1051,7 @@ def override_structure(
         center : CoordinateOptional
             Center of the override structure. `None` coordinate along an axis means refinement is not
             applied along that axis.
-        grid_size_in_vaccum : float
+        grid_size_in_vacuum : float
             Grid step size in vaccum.
         drop_outside_sim : bool
             Drop override structures outside simulation domain.
@@ -1332,20 +1426,12 @@ def center_axis(self) -> float:
         """Gets the position of the center of the layer along the layer dimension."""
         return self.center[self.axis]
 
-    @cached_property
-    def _is_inplane_bounded(self) -> bool:
-        """Whether the layer is bounded in at least one of the inplane dimensions."""
-        return np.isfinite(self.size[(self.axis + 1) % 3]) or np.isfinite(
-            self.size[(self.axis + 2) % 3]
+    def _is_inplane_bounded(self, geometry: Box) -> bool:
+        """Whether the geometry is bounded in at least one of the inplane dimensions."""
+        return np.isfinite(geometry.size[(self.axis + 1) % 3]) or np.isfinite(
+            geometry.size[(self.axis + 2) % 3]
         )
 
-    @cached_property
-    def _slightly_enlarged_box(self) -> Box:
-        """Slightly enlarged box for robust point containment querying."""
-        # increase size slightly
-        size = [increment_float(orig_length, 1) for orig_length in self.size]
-        return Box(center=self.center, size=size)
-
     def _unpop_axis(self, ax_coord: float, plane_coord: Any) -> CoordinateOptional:
         """Combine coordinate along axis with identical coordinates on the plane tangential to the axis.
 
@@ -1363,13 +1449,32 @@ def _unpop_axis(self, ax_coord: float, plane_coord: Any) -> CoordinateOptional:
         """
         return self.unpop_axis(ax_coord, [plane_coord, plane_coord], self.axis)
 
-    def suggested_dl_min(self, grid_size_in_vacuum: float, structures: list[Structure]) -> float:
+    def suggested_dl_min(
+        self,
+        grid_size_in_vacuum: float,
+        structures: list[Structure],
+        sim_bounds: tuple,
+        boundary_type: tuple,
+        cached_merged_geos=None,
+        cached_corners_and_convexity=None,
+    ) -> float:
         """Suggested lower bound of grid step size for this layer.
 
         Parameters
         ----------
-        grid_size_in_vaccum : float
+        grid_size_in_vacuum : float
             Grid step size in vaccum.
+        structures : list[Structure]
+            List of structures present in simulation.
+        sim_bounds : tuple
+            Bounds of the simulation domain excluding the PML regions, formatted as
+            ``(mins, maxs)`` where each is a 3-tuple of coordinates.
+        boundary_type : tuple
+            Boundary type of the simulation domain.
+        cached_merged_geos : Optional[list[tuple[Any, Shapely]]]
+            Cached merged geometries. If None, will be computed.
+        cached_corners_and_convexity : Optional[tuple[list[ArrayFloat2D], list[ArrayFloat1D]]]
+            Cached corners and convexity data. If None, will be computed.
 
         Returns
         -------
@@ -1397,38 +1502,65 @@ def suggested_dl_min(self, grid_size_in_vacuum: float, structures: list[Structur
 
         # min feature size
         if self.corner_finder is not None and not self.corner_finder._no_min_dl_override:
-            dl_suggested = self._dl_min_from_smallest_feature(structures)
+            dl_suggested = self._dl_min_from_smallest_feature(
+                structures,
+                sim_bounds=sim_bounds,
+                boundary_type=boundary_type,
+                cached_merged_geos=cached_merged_geos,
+                cached_corners_and_convexity=cached_corners_and_convexity,
+            )
             dl_min = min(dl_min, dl_suggested)
 
         return dl_min
 
     def generate_snapping_points(
-        self, structure_list: list[Structure], cached_corners_and_convexity=None
+        self,
+        structure_list: list[Structure],
+        sim_bounds: tuple,
+        boundary_type: tuple,
+        cached_corners_and_convexity=None,
+        cached_merged_geos=None,
     ) -> list[CoordinateOptional]:
         """generate snapping points for mesh refinement."""
         snapping_points = self._snapping_points_along_axis
         if self.corner_snapping:
-            snapping_points += self._corners(structure_list, cached_corners_and_convexity)
+            snapping_points += self._corners(
+                structure_list,
+                sim_bounds=sim_bounds,
+                boundary_type=boundary_type,
+                cached_corners_and_convexity=cached_corners_and_convexity,
+                cached_merged_geos=cached_merged_geos,
+            )
         return snapping_points
 
     def generate_override_structures(
         self,
         grid_size_in_vacuum: float,
         structure_list: list[Structure],
+        sim_bounds: tuple,
+        boundary_type: tuple,
         cached_corners_and_convexity=None,
+        cached_merged_geos=None,
     ) -> list[MeshOverrideStructure]:
         """Generate mesh override structures for mesh refinement."""
         return self._override_structures_along_axis(
             grid_size_in_vacuum
         ) + self._override_structures_inplane(
-            structure_list, grid_size_in_vacuum, cached_corners_and_convexity
+            structure_list,
+            grid_size_in_vacuum,
+            sim_bounds,
+            boundary_type,
+            cached_corners_and_convexity,
+            cached_merged_geos,
         )
 
-    def _inplane_inside(self, point: ArrayFloat2D) -> bool:
-        """On the inplane cross section, whether the point is inside the layer.
+    def _inplane_inside(self, box: Box, point: ArrayFloat2D) -> bool:
+        """On the inplane cross section, whether the point is inside the box.
 
         Parameters
         ----------
+        box : Box
+            Box to check point containment.
         point : ArrayFloat2D
             Point position on inplane plane.
 
@@ -1441,35 +1573,135 @@ def _inplane_inside(self, point: ArrayFloat2D) -> bool:
         point_3d = self.unpop_axis(
             ax_coord=self.center[self.axis], plane_coords=point, axis=self.axis
         )
-        return self._slightly_enlarged_box.inside(point_3d[0], point_3d[1], point_3d[2])
+        return box.inside(point_3d[0], point_3d[1], point_3d[2])
+
+    def _layer_box(self, sim_bounds: tuple, boundary_type: tuple) -> Box:
+        """Layer box with size slightly adjusted to boundary conditions.
+
+        Parameters
+        ----------
+        sim_bounds : tuple
+            Bounds of the simulation domain excluding the PML regions, formatted as
+            ``(mins, maxs)`` where each is a 3-tuple of coordinates.
+        boundary_type : tuple
+            Boundary type of the simulation domain.
+
+        Returns
+        -------
+        Box
+            Layer box with size slightly adjusted to boundary conditions.
+        """
+        layer_geometry = Box(center=self.center, size=self.size)
+
+        # start from the actual bounds of the layer geometry
+        layer_min, layer_max = (list(b) for b in layer_geometry.bounds)
+        sim_min, sim_max = sim_bounds
+
+        for axis2d_ind in range(2):
+            axis_ind = (self.axis + axis2d_ind + 1) % 3
+            bc_minus, bc_plus = boundary_type[axis_ind]
+            # shrink away from lower periodic boundary if the layer touches it
+            if bc_minus == "periodic" and layer_min[axis_ind] <= sim_min[axis_ind] + fp_eps:
+                layer_min[axis_ind] = sim_min[axis_ind] + fp_eps
+
+            # shrink away from upper periodic boundary if the layer touches it
+            if bc_plus == "periodic" and layer_max[axis_ind] >= sim_max[axis_ind] - fp_eps:
+                layer_max[axis_ind] = sim_max[axis_ind] - fp_eps
+
+        new_box = Box.from_bounds(rmin=layer_min, rmax=layer_max)
+        return layer_geometry.updated_copy(center=new_box.center, size=new_box.size)
+
+    def _merged_geos(
+        self, structure_list: list[Structure], sim_bounds: tuple, boundary_type: tuple
+    ) -> list[tuple[Any, Shapely]]:
+        """Merged geometries on the inplane plane.
+
+        Parameters
+        ----------
+        structure_list : list[Structure]
+            List of structures present in simulation.
+        sim_bounds : tuple
+            Bounds of the simulation domain excluding the PML regions.
+        boundary_type : tuple
+            Boundary type of the simulation domain.
+
+        Returns
+        -------
+        list[tuple[Any, Shapely]]
+            Merged geometries on the inplane plane.
+        """
+        if self.corner_finder is None:
+            return []
+
+        layer_geometry = self._layer_box(sim_bounds, boundary_type)
+        # filter structures outside the layer
+        structures_intersect = structure_list
+        if self._is_inplane_bounded(layer_geometry):
+            structures_intersect = [
+                s for s in structure_list if layer_geometry.intersects(s.geometry)
+            ]
+        merged_geos = self.corner_finder._merged_pec_on_plane(
+            normal_axis=self.axis,
+            coord=self.center_axis,
+            structure_list=structures_intersect,
+            center=layer_geometry.center,
+            size=layer_geometry.size,
+            interior_disjoint_geometries=self.interior_disjoint_geometries,
+            keep_metal_only=(
+                self.interior_disjoint_geometries and self.corner_finder.medium == "metal"
+            ),
+        )
+        return merged_geos
 
     def _corners_and_convexity_2d(
-        self, structure_list: list[Structure], ravel: bool
+        self,
+        merged_geos: list[tuple[Any, Shapely]],
+        structure_list: list[Structure],
+        ravel: bool,
+        sim_bounds: tuple,
+        boundary_type: tuple,
     ) -> tuple[list[ArrayFloat2D], list[ArrayFloat1D]]:
-        """Raw inplane corners and their convexity."""
+        """Raw inplane corners and their convexity.
+
+        Parameters
+        ----------
+        merged_geos : list[tuple[Any, Shapely]]
+            Merged geometries on the inplane plane.
+        structure_list : list[Structure]
+            List of structures present in simulation.
+        ravel : bool
+            Whether to put the resulting corners in a single list or per polygon.
+        sim_bounds : tuple
+            Bounds of the simulation domain excluding the PML regions.
+        boundary_type : tuple
+            Boundary type of the simulation domain.
+
+        Returns
+        -------
+        tuple[list[ArrayFloat2D], list[ArrayFloat1D]]
+            Raw inplane corners and their convexity.
+        """
         if self.corner_finder is None:
             return [], []
 
-        # filter structures outside the layer
-        structures_intersect = structure_list
-        if self._is_inplane_bounded:
-            structures_intersect = [s for s in structure_list if self.intersects(s.geometry)]
         inplane_points, convexity = self.corner_finder._corners_and_convexity(
-            self.axis,
-            self.center_axis,
-            structures_intersect,
-            ravel,
-            self.interior_disjoint_geometries,
+            merged_geos=merged_geos,
+            ravel=ravel,
         )
 
         # filter corners outside the inplane bounds
-        if self._is_inplane_bounded and len(inplane_points) > 0:
+        layer_geometry = self._layer_box(sim_bounds, boundary_type)
+        layer_geometry_slightly_enlarged = layer_geometry.slightly_enlarged_copy()
+        if self._is_inplane_bounded(layer_geometry) and len(inplane_points) > 0:
             # flatten temporary list of arrays for faster processing
             if not ravel:
                 split_inds = np.cumsum([len(pts) for pts in inplane_points])[:-1]
                 inplane_points = np.concatenate(inplane_points)
                 convexity = np.concatenate(convexity)
-            inds = [self._inplane_inside(point) for point in inplane_points]
+            inds = [
+                self._inplane_inside(layer_geometry_slightly_enlarged, point)
+                for point in inplane_points
+            ]
             inplane_points = inplane_points[inds]
             convexity = convexity[inds]
             if not ravel:
@@ -1478,12 +1710,30 @@ def _corners_and_convexity_2d(
 
         return inplane_points, convexity
 
-    def _dl_min_from_smallest_feature(self, structure_list: list[Structure]):
+    def _dl_min_from_smallest_feature(
+        self,
+        structure_list: list[Structure],
+        sim_bounds: tuple,
+        boundary_type: tuple,
+        cached_merged_geos=None,
+        cached_corners_and_convexity=None,
+    ):
         """Calculate `dl_min` suggestion based on smallest feature size."""
 
-        inplane_points, convexity = self._corners_and_convexity_2d(
-            structure_list=structure_list, ravel=False
-        )
+        if cached_corners_and_convexity is None:
+            if cached_merged_geos is None:
+                merged_geos = self._merged_geos(structure_list, sim_bounds, boundary_type)
+            else:
+                merged_geos = cached_merged_geos
+            inplane_points, convexity = self._corners_and_convexity_2d(
+                merged_geos=merged_geos,
+                structure_list=structure_list,
+                ravel=False,
+                sim_bounds=sim_bounds,
+                boundary_type=boundary_type,
+            )
+        else:
+            inplane_points, convexity = cached_corners_and_convexity
 
         dl_min = inf
 
@@ -1517,14 +1767,27 @@ def _dl_min_from_smallest_feature(self, structure_list: list[Structure]):
         return dl_min
 
     def _corners(
-        self, structure_list: list[Structure], cached_corners_and_convexity=None
+        self,
+        structure_list: list[Structure],
+        sim_bounds: tuple,
+        boundary_type: tuple,
+        cached_corners_and_convexity=None,
+        cached_merged_geos=None,
     ) -> list[CoordinateOptional]:
         """Inplane corners in 3D coordinate."""
         if self.corner_finder is None:
             return []
         if cached_corners_and_convexity is None:
+            if cached_merged_geos is None:
+                merged_geos = self._merged_geos(structure_list, sim_bounds, boundary_type)
+            else:
+                merged_geos = cached_merged_geos
             inplane_points, _ = self._corners_and_convexity_2d(
-                structure_list=structure_list, ravel=True
+                merged_geos=merged_geos,
+                structure_list=structure_list,
+                ravel=True,
+                sim_bounds=sim_bounds,
+                boundary_type=boundary_type,
             )
         else:
             inplane_points, convexity = cached_corners_and_convexity
@@ -1567,7 +1830,10 @@ def _override_structures_inplane(
         self,
         structure_list: list[Structure],
         grid_size_in_vacuum: float,
+        sim_bounds: tuple,
+        boundary_type: tuple,
         cached_corners_and_convexity=None,
+        cached_merged_geos=None,
     ) -> list[MeshOverrideStructure]:
         """Inplane mesh override structures for refining mesh around corners."""
         if self.corner_refinement is None:
@@ -1577,7 +1843,13 @@ def _override_structures_inplane(
             self.corner_refinement.override_structure(
                 corner, grid_size_in_vacuum, self.refinement_inside_sim_only
             )
-            for corner in self._corners(structure_list, cached_corners_and_convexity)
+            for corner in self._corners(
+                structure_list,
+                sim_bounds=sim_bounds,
+                boundary_type=boundary_type,
+                cached_corners_and_convexity=cached_corners_and_convexity,
+                cached_merged_geos=cached_merged_geos,
+            )
         ]
 
     def _override_structures_along_axis(
@@ -1991,18 +2263,18 @@ def _generate_horizontal_snapping_lines(
         return snapping_lines_y, min_gap_width
 
     def _resolve_gaps(
-        self, structures: list[Structure], grid: Grid, boundary_types: tuple
+        self, grid: Grid, merged_geos: list[tuple[Any, Shapely]], boundary_type: tuple
     ) -> tuple[list[CoordinateOptional], float]:
         """
         Detect underresolved gaps and place snapping lines in them. Also return the detected minimal gap width.
 
         Parameters
         ----------
-        structures : list[Structure]
-            List of structures to consider.
         grid : Grid
             Grid to resolve gaps on.
-        boundary_types : Tuple[Tuple[str, str], Tuple[str, str], Tuple[str, str]] = [[None, None], [None, None], [None, None]]
+        merged_geos : list[tuple[Any, Shapely]]
+            Merged geometries on the inplane plane.
+        boundary_type : tuple[tuple[str, str], tuple[str, str], tuple[str, str]]
             Type of boundary conditions along each dimension: "pec/pmc", "periodic", or
             None for any other. This is relevant only for gap meshing.
 
@@ -2017,13 +2289,15 @@ def _resolve_gaps(
         x = grid.boundaries.to_list[tan_dims[0]]
         y = grid.boundaries.to_list[tan_dims[1]]
 
-        _, boundaries_tan = Box.pop_axis(boundary_types, self.axis)
-
         # restrict to the size of layer spec
         rmin, rmax = self.bounds
         _, rmin = Box.pop_axis(rmin, self.axis)
         _, rmax = Box.pop_axis(rmax, self.axis)
 
+        # extract tangential boundaries
+        _, boundaries_tan = Box.pop_axis(boundary_type, self.axis)
+
+        # new coords are expanded by a grid at both min/max
         new_coords = []
         new_boundaries = []
         for coord, cmin, cmax, bdry in zip([x, y], rmin, rmax, boundaries_tan):
@@ -2032,12 +2306,12 @@ def _resolve_gaps(
             if cmin < coord[0]:
                 ind_min = 0
             else:
-                ind_min = max(0, np.argmax(coord >= cmin) - 1)
+                ind_min = max(0, np.argmax(coord >= cmin) - 2)
 
             if cmax > coord[-1]:
                 ind_max = len(coord) - 1
             else:
-                ind_max = np.argmax(coord >= cmax)
+                ind_max = np.argmax(coord >= cmax) + 1
 
             if ind_min >= ind_max - 1:
                 return [], inf
@@ -2053,36 +2327,6 @@ def _resolve_gaps(
 
         x, y = new_coords
 
-        merging_area_bounds = np.array(
-            [[x[0] - fp_eps, y[0] - fp_eps], [x[-1] + fp_eps, y[-1] + fp_eps]]
-        )
-
-        # restrict size of the plane where pec polygons are found in case of periodic boundary conditions
-        # this is to make sure gaps across periodic boundary conditions are resolved
-        # (if there is a PEC structure going into periodic boundary, now it will generate a grid line
-        # intersection next to that boundary and it will be propagated to the other side)
-        for ind in range(2):
-            if new_boundaries[ind][0] == "periodic":
-                merging_area_bounds[0][ind] += fp_eps + dp_eps
-                merging_area_bounds[1][ind] -= fp_eps + dp_eps
-
-        merging_area_center = 0.5 * (merging_area_bounds[0] + merging_area_bounds[1])
-        merging_area_size = merging_area_bounds[1] - merging_area_bounds[0]
-
-        merging_area_center = Box.unpop_axis(self.center[self.axis], merging_area_center, self.axis)
-        merging_area_size = Box.unpop_axis(self.size[self.axis], merging_area_size, self.axis)
-
-        # get merged pec structures on plane
-        # note that we expect this function to also convert all LossyMetal's into PEC
-        plane_slice = CornerFinderSpec._merged_pec_on_plane(
-            coord=self.center_axis,
-            normal_axis=self.axis,
-            structure_list=structures,
-            center=merging_area_center,
-            size=merging_area_size,
-            interior_disjoint_geometries=self.interior_disjoint_geometries,
-        )
-
         # find intersections of pec polygons with grid lines
         # specifically:
         # 0. cells that contain intersections of vertical grid lines
@@ -2090,7 +2334,7 @@ def _resolve_gaps(
         # 2. cells that contain intersections of horizontal grid lines
         # 3. relative locations of those intersections along x axis
         v_cells_ij, v_cells_dy, h_cells_ij, h_cells_dx = self._process_slice(
-            x, y, plane_slice, new_boundaries
+            x, y, merged_geos, new_boundaries
         )
 
         # generate horizontal snapping lines
@@ -2302,7 +2546,10 @@ def internal_snapping_points(
         self,
         structures: list[Structure],
         lumped_elements: list[LumpedElementType],
+        boundary_types: tuple[tuple[str, str], tuple[str, str], tuple[str, str]],
+        sim_bounds: tuple,
         cached_corners_and_convexity=None,
+        cached_merged_geos=None,
     ) -> list[CoordinateOptional]:
         """Internal snapping points. So far, internal snapping points are generated by
         `layer_refinement_specs` and lumped element.
@@ -2313,8 +2560,14 @@ def internal_snapping_points(
             List of physical structures.
         lumped_elements : List[LumpedElementType]
             List of lumped elements.
+        boundary_types : tuple[tuple[str, str], tuple[str, str], tuple[str, str]]
+            Boundary type of the simulation domain.
         cached_corners_and_convexity : Optional[list[CachedCornersAndConvexity]]
             Cached corners and convexity data.
+        cached_merged_geos : Optional[list[list[tuple[Any, Shapely]]]]
+            Cached merged geometries for each layer. If None, will be computed.
+        sim_bounds : tuple
+            Bounds of the simulation domain excluding the PML regions.
 
         Returns
         -------
@@ -2334,8 +2587,16 @@ def internal_snapping_points(
                     cached_data = cached_corners_and_convexity[ind]
                 else:
                     cached_data = None
+                if cached_merged_geos is not None:
+                    cached_merged = cached_merged_geos[ind]
+                else:
+                    cached_merged = None
                 snapping_points += layer_spec.generate_snapping_points(
-                    list(structures), cached_data
+                    list(structures),
+                    sim_bounds,
+                    boundary_types,
+                    cached_data,
+                    cached_merged,
                 )
         # ) from lumped_elements
         for lumped_element in lumped_elements:
@@ -2346,6 +2607,8 @@ def all_snapping_points(
         self,
         structures: list[Structure],
         lumped_elements: list[LumpedElementType],
+        boundary_types: tuple[tuple[str, str], tuple[str, str], tuple[str, str]],
+        sim_bounds: tuple,
         internal_snapping_points: Optional[list[CoordinateOptional]] = None,
     ) -> list[CoordinateOptional]:
         """Internal and external snapping points. External snapping points take higher priority.
@@ -2357,6 +2620,8 @@ def all_snapping_points(
             List of physical structures.
         lumped_elements : List[LumpedElementType]
             List of lumped elements.
+        boundary_types : tuple[tuple[str, str], tuple[str, str], tuple[str, str]]
+            Boundary type of the simulation domain.
         internal_snapping_points : List[CoordinateOptional]
             If `None`, recomputes internal snapping points.
 
@@ -2367,9 +2632,12 @@ def all_snapping_points(
         """
 
         if internal_snapping_points is None:
-            return self.internal_snapping_points(structures, lumped_elements) + list(
-                self.snapping_points
-            )
+            return self.internal_snapping_points(
+                structures,
+                lumped_elements,
+                boundary_types,
+                sim_bounds,
+            ) + list(self.snapping_points)
         return internal_snapping_points + list(self.snapping_points)
 
     @property
@@ -2381,9 +2649,11 @@ def internal_override_structures(
         self,
         structures: list[Structure],
         wavelength: pd.PositiveFloat,
-        sim_size: tuple[float, 3],
+        sim_bounds: tuple,
         lumped_elements: list[LumpedElementType],
+        boundary_types: tuple[tuple[str, str], tuple[str, str], tuple[str, str]],
         cached_corners_and_convexity=None,
+        cached_merged_geos=None,
     ) -> list[StructureType]:
         """Internal mesh override structures. So far, internal override structures are generated by
         `layer_refinement_specs` and lumped element.
@@ -2394,12 +2664,16 @@ def internal_override_structures(
             List of structures, with the simulation structure being the first item.
         wavelength : pd.PositiveFloat
             Wavelength to use for minimal step size in vaccum.
-        sim_size : Tuple[float, 3]
-            Simulation domain size.
         lumped_elements : List[LumpedElementType]
             List of lumped elements.
+        boundary_types : tuple[tuple[str, str], tuple[str, str], tuple[str, str]]
+            Boundary type of the simulation domain.
         cached_corners_and_convexity : Optional[list[CachedCornersAndConvexity]]
             Cached corners and convexity data.
+        cached_merged_geos : Optional[list[list[tuple[Any, Shapely]]]]
+            Cached merged geometries for each layer. If None, will be computed.
+        sim_bounds : tuple
+            Bounds of the simulation domain excluding the PML regions.
 
         Returns
         -------
@@ -2419,10 +2693,19 @@ def internal_override_structures(
                     cached_data = cached_corners_and_convexity[ind]
                 else:
                     cached_data = None
+                if cached_merged_geos is not None:
+                    cached_merged = cached_merged_geos[ind]
+                else:
+                    cached_merged = None
+                # use simulation size derived from bounds for vacuum dl
+                sim_size_local = tuple(bmax - bmin for bmin, bmax in zip(*sim_bounds))
                 override_structures += layer_spec.generate_override_structures(
-                    self._min_vacuum_dl_in_autogrid(wavelength, sim_size),
+                    self._min_vacuum_dl_in_autogrid(wavelength, sim_size_local),
                     list(structures),
+                    sim_bounds,
+                    boundary_types,
                     cached_data,
+                    cached_merged,
                 )
         # 2) from lumped element
         for lumped_element in lumped_elements:
@@ -2433,8 +2716,9 @@ def all_override_structures(
         self,
         structures: list[Structure],
         wavelength: pd.PositiveFloat,
-        sim_size: tuple[float, 3],
         lumped_elements: list[LumpedElementType],
+        boundary_types: tuple[tuple[str, str], tuple[str, str], tuple[str, str]],
+        sim_bounds: tuple,
         structure_priority_mode: PriorityMode = "equal",
         internal_override_structures: Optional[list[MeshOverrideStructure]] = None,
     ) -> list[StructureType]:
@@ -2447,15 +2731,16 @@ def all_override_structures(
             List of structures, with the simulation structure being the first item.
         wavelength : pd.PositiveFloat
             Wavelength to use for minimal step size in vaccum.
-        sim_size : Tuple[float, 3]
-            Simulation domain size.
         lumped_elements : List[LumpedElementType]
             List of lumped elements.
+        sim_bounds : tuple
+            Bounds of the simulation domain excluding the PML regions.
+        boundary_types : tuple[tuple[str, str], tuple[str, str], tuple[str, str]]
+            Boundary type of the simulation domain.
         structure_priority_mode : PriorityMode
             Structure priority setting.
         internal_override_structures : List[MeshOverrideStructure]
             If `None`, recomputes internal override structures.
-
         Returns
         -------
         List[StructureType]
@@ -2464,11 +2749,59 @@ def all_override_structures(
 
         if internal_override_structures is None:
             internal_override_structures = self.internal_override_structures(
-                structures, wavelength, sim_size, lumped_elements
+                structures,
+                wavelength,
+                sim_bounds,
+                lumped_elements,
+                boundary_types,
             )
         all_structures = internal_override_structures + self.external_override_structures
         return Structure._sort_structures(all_structures, structure_priority_mode)
 
+    def _get_all_structures_affecting_grid(
+        self,
+        structures: list[Structure],
+        wavelength: pd.PositiveFloat,
+        lumped_elements: list[LumpedElementType],
+        boundary_types: tuple[tuple[str, str], tuple[str, str], tuple[str, str]],
+        sim_bounds: tuple,
+        structure_priority_mode: PriorityMode = "equal",
+        internal_override_structures: Optional[list[MeshOverrideStructure]] = None,
+    ) -> list[StructureType]:
+        """Get all structures including original structures and all override structures.
+
+        Parameters
+        ----------
+        structures : List[Structure]
+            List of structures, with the simulation structure being the first item.
+        wavelength : pd.PositiveFloat
+            Wavelength to use for minimal step size in vacuum.
+        lumped_elements : List[LumpedElementType]
+            List of lumped elements.
+        boundary_types : tuple[tuple[str, str], tuple[str, str], tuple[str, str]]
+            Boundary type of the simulation domain.
+        sim_bounds : tuple
+            Bounds of the simulation domain excluding the PML regions.
+        structure_priority_mode : PriorityMode
+            Structure priority setting.
+        internal_override_structures : List[MeshOverrideStructure]
+            If `None`, recomputes internal override structures.
+
+        Returns
+        -------
+        List[StructureType]
+            List of all structures including original structures and override structures.
+        """
+        return list(structures) + self.all_override_structures(
+            list(structures),
+            wavelength,
+            lumped_elements,
+            boundary_types,
+            sim_bounds,
+            structure_priority_mode,
+            internal_override_structures,
+        )
+
     def _min_vacuum_dl_in_autogrid(self, wavelength: float, sim_size: tuple[float, 3]) -> float:
         """Compute grid step size in vacuum for Autogrd. If AutoGrid is applied along more than 1 dimension,
         return the minimal.
@@ -2483,8 +2816,10 @@ def _dl_min(
         self,
         wavelength: float,
         structure_list: list[StructureType],
-        sim_size: tuple[float, 3],
+        sim_bounds: tuple,
         lumped_elements: list[LumpedElementType],
+        boundary_types: tuple[tuple[str, str], tuple[str, str], tuple[str, str]],
+        cached_merged_geos: Optional[list[list[tuple[Any, Shapely]]]] = None,
     ) -> float:
         """Lower bound of grid size to be applied to dimensions where AutoGrid with unset
         `dl_min` (0 or None) is applied.
@@ -2504,15 +2839,28 @@ def _dl_min(
             for dl in structure.dl:
                 if dl is not None and dl < min_dl:
                     min_dl = dl
+        # local simulation size derived from bounds
+        sim_size_local = tuple(bmax - bmin for bmin, bmax in zip(*sim_bounds))
+
         # from mesh specification
         for grid in [self.grid_x, self.grid_y, self.grid_z]:
-            min_dl = min(min_dl, grid.estimated_min_dl(wavelength, structures, sim_size))
+            min_dl = min(min_dl, grid.estimated_min_dl(wavelength, structures, sim_size_local))
 
         # from layer refinement specifications
         if self.layer_refinement_used:
-            min_vacuum_dl = self._min_vacuum_dl_in_autogrid(wavelength, sim_size)
-            for layer in self.layer_refinement_specs:
-                min_dl = min(min_dl, layer.suggested_dl_min(min_vacuum_dl, structures))
+            min_vacuum_dl = self._min_vacuum_dl_in_autogrid(wavelength, sim_size_local)
+            for ind, layer in enumerate(self.layer_refinement_specs):
+                cached_merged = cached_merged_geos[ind] if cached_merged_geos is not None else None
+                min_dl = min(
+                    min_dl,
+                    layer.suggested_dl_min(
+                        min_vacuum_dl,
+                        structures,
+                        sim_bounds,
+                        boundary_types,
+                        cached_merged,
+                    ),
+                )
         # from lumped elements
         for lumped_element in lumped_elements:
             for override_structure in lumped_element.to_mesh_overrides():
@@ -2609,9 +2957,10 @@ def _make_grid_and_snapping_lines(
             [None, None],
         ],
         structure_priority_mode: PriorityMode = "equal",
+        cached_merged_geos: Optional[list[list[tuple[Any, Shapely]]]] = None,
     ) -> tuple[Grid, list[CoordinateOptional]]:
         """Make the entire simulation grid based on some simulation parameters.
-        Also return snappiung point resulted from iterative gap meshing.
+        Also return snapping point resulted from iterative gap meshing.
 
         Parameters
         ----------
@@ -2639,6 +2988,8 @@ def _make_grid_and_snapping_lines(
             None for any other. This is relevant only for gap meshing.
         structure_priority_mode : PriorityMode
             Structure priority setting.
+        cached_merged_geos : Optional[list[list[tuple[Any, Shapely]]]]
+            Cached merged geometries for each layer. If None, will be computed.
 
         Returns
         -------
@@ -2646,6 +2997,45 @@ def _make_grid_and_snapping_lines(
             Entire simulation grid and snapping points generated during iterative gap meshing.
         """
 
+        # Pre-compute results from parse_structures
+        wavelength = self.get_wavelength(sources)
+        sim_bounds = structures[0].geometry.bounds
+        all_structures = self._get_all_structures_affecting_grid(
+            structures,
+            wavelength,
+            lumped_elements,
+            boundary_types,
+            sim_bounds,
+            structure_priority_mode,
+            internal_override_structures,
+        )
+
+        parse_structures_interval_coords = []
+        parse_structures_max_dl_list = []
+        grids_1d = [self.grid_x, self.grid_y, self.grid_z]
+
+        for idim, grid_1d in enumerate(grids_1d):
+            if isinstance(grid_1d, AbstractAutoGrid):
+                interval_coords, max_dl_list = grid_1d._parse_structures(
+                    axis=idim,
+                    structures=all_structures,
+                    wavelength=wavelength,
+                    symmetry=symmetry,
+                    snapping_points=self.all_snapping_points(
+                        structures,
+                        lumped_elements,
+                        boundary_types,
+                        sim_bounds,
+                        internal_snapping_points,
+                    ),
+                )
+                parse_structures_interval_coords.append(interval_coords)
+                parse_structures_max_dl_list.append(max_dl_list)
+            else:
+                # For non-AutoGrid, append None since they don't use _parse_structures
+                parse_structures_interval_coords.append(None)
+                parse_structures_max_dl_list.append(None)
+
         old_grid = self._make_grid_one_iteration(
             structures=structures,
             symmetry=symmetry,
@@ -2656,8 +3046,13 @@ def _make_grid_and_snapping_lines(
             internal_override_structures=internal_override_structures,
             internal_snapping_points=internal_snapping_points,
             structure_priority_mode=structure_priority_mode,
+            boundary_types=boundary_types,
+            parse_structures_interval_coords=parse_structures_interval_coords,
+            parse_structures_max_dl_list=parse_structures_max_dl_list,
         )
 
+        # gap refinement only place snapping points, and decrease dl_min that only affects
+        # snapping points insertion.
         snapping_lines = []
         if len(self.layer_refinement_specs) > 0:
             num_iters = max(
@@ -2667,11 +3062,21 @@ def _make_grid_and_snapping_lines(
             min_gap_width = inf
             for ind in range(num_iters):
                 new_snapping_lines = []
-                for layer_spec in self.layer_refinement_specs:
+                sim_bounds = structures[0].geometry.bounds
+                for ind_layer, layer_spec in enumerate(self.layer_refinement_specs):
                     if layer_spec.gap_meshing_iters > ind:
+                        # use cached merged geometries if available, otherwise compute
+                        if cached_merged_geos is not None:
+                            merged_geos = cached_merged_geos[ind_layer]
+                        else:
+                            merged_geos = layer_spec._merged_geos(
+                                structures,
+                                sim_bounds,
+                                boundary_types,
+                            )
                         one_layer_snapping_lines, gap_width = layer_spec._resolve_gaps(
-                            structures,
                             old_grid,
+                            merged_geos,
                             boundary_types,
                         )
                         new_snapping_lines = new_snapping_lines + one_layer_snapping_lines
@@ -2698,6 +3103,9 @@ def _make_grid_and_snapping_lines(
                     internal_snapping_points=snapping_lines + internal_snapping_points,
                     dl_min_from_gaps=0.45 * min_gap_width,
                     structure_priority_mode=structure_priority_mode,
+                    boundary_types=boundary_types,
+                    parse_structures_interval_coords=parse_structures_interval_coords,
+                    parse_structures_max_dl_list=parse_structures_max_dl_list,
                 )
 
                 same = old_grid == new_grid
@@ -2720,11 +3128,14 @@ def _make_grid_one_iteration(
         periodic: tuple[bool, bool, bool],
         sources: list[SourceType],
         num_pml_layers: list[tuple[pd.NonNegativeInt, pd.NonNegativeInt]],
+        boundary_types: tuple[tuple[str, str], tuple[str, str], tuple[str, str]],
         lumped_elements: list[LumpedElementType] = (),
         internal_override_structures: Optional[list[MeshOverrideStructure]] = None,
         internal_snapping_points: Optional[list[CoordinateOptional]] = None,
         dl_min_from_gaps: pd.PositiveFloat = inf,
         structure_priority_mode: PriorityMode = "equal",
+        parse_structures_interval_coords: Optional[list[np.ndarray]] = None,
+        parse_structures_max_dl_list: Optional[list[np.ndarray]] = None,
     ) -> Grid:
         """Make the entire simulation grid based on some simulation parameters.
 
@@ -2753,6 +3164,12 @@ def _make_grid_one_iteration(
             Minimal grid size computed based on autodetected gaps.
         structure_priority_mode : PriorityMode
             Structure priority setting.
+        parse_structures_interval_coords : Optional[List[np.ndarray]]
+            If not None, pre-computed interval coordinates from parsing structures for each dimension.
+            List of length 3, one for each axis (x, y, z).
+        parse_structures_max_dl_list : Optional[List[np.ndarray]]
+            If not None, pre-computed maximum grid spacing list from parsing structures for each dimension.
+            List of length 3, one for each axis (x, y, z).
 
         Returns
         -------
@@ -2809,12 +3226,13 @@ def _make_grid_one_iteration(
                 "an instance without any autograd tracers."
             )
 
-        sim_size = list(structures[0].geometry.size)
-        all_structures = list(structures) + self.all_override_structures(
-            list(structures),
+        sim_bounds = structures[0].geometry.bounds
+        all_structures = self._get_all_structures_affecting_grid(
+            structures,
             wavelength,
-            sim_size,
             lumped_elements,
+            boundary_types,
+            sim_bounds,
             structure_priority_mode,
             internal_override_structures,
         )
@@ -2829,8 +3247,9 @@ def _make_grid_one_iteration(
             new_dl_min = self._dl_min(
                 wavelength,
                 list(structures) + self.external_override_structures,
-                sim_size,
+                sim_bounds,
                 lumped_elements,
+                boundary_types,
             )
             new_dl_min = min(new_dl_min, dl_min_from_gaps)
             for ind, grid in enumerate(grids_1d):
@@ -2847,7 +3266,21 @@ def _make_grid_one_iteration(
                 wavelength=wavelength,
                 num_pml_layers=num_pml_layers[idim],
                 snapping_points=self.all_snapping_points(
-                    structures, lumped_elements, internal_snapping_points
+                    structures,
+                    lumped_elements,
+                    boundary_types,
+                    sim_bounds,
+                    internal_snapping_points,
+                ),
+                parse_structures_interval_coords=(
+                    parse_structures_interval_coords[idim]
+                    if parse_structures_interval_coords is not None
+                    else None
+                ),
+                parse_structures_max_dl_list=(
+                    parse_structures_max_dl_list[idim]
+                    if parse_structures_max_dl_list is not None
+                    else None
                 ),
             )
 
diff --git a/tidy3d/components/simulation.py b/tidy3d/components/simulation.py
index b82b57d36d..86c57d71a1 100644
--- a/tidy3d/components/simulation.py
+++ b/tidy3d/components/simulation.py
@@ -126,6 +126,7 @@
     FreqBound,
     InterpMethod,
     PermittivityComponent,
+    Shapely,
     Symmetry,
     annotate_type,
 )
@@ -978,16 +979,48 @@ def pml_thicknesses(self) -> list[tuple[float, float]]:
 
         return pml_thicknesses
 
+    @cached_property
+    def _internal_layerrefinement_boundary_types(self):
+        """Boundary types for layer refinement."""
+        boundary_types = [[None, None], [None, None], [None, None]]
+        for dim, boundary in enumerate(self.boundary_spec.to_list):
+            for side, edge in enumerate(boundary):
+                if isinstance(edge, (PECBoundary, PMCBoundary)):
+                    boundary_types[dim][side] = "pec/pmc"
+                elif isinstance(edge, (Periodic, BlochBoundary)):
+                    boundary_types[dim][side] = "periodic"
+        return boundary_types
+
+    @cached_property
+    def _internal_layerrefinement_merged_geos(self) -> list[tuple[Any, Shapely]]:
+        """Merged geometries on the plane for each layer refinement specification."""
+        cached_data = []
+        for layer in self.grid_spec.layer_refinement_specs:
+            cached_data.append(
+                layer._merged_geos(
+                    structure_list=self.scene.all_structures,
+                    sim_bounds=self.bounds,
+                    boundary_type=self._internal_layerrefinement_boundary_types,
+                )
+            )
+        return cached_data
+
     @cached_property
     def _internal_layerfinement_corners_and_convexity_2d(
         self,
     ) -> list[tuple[list[ArrayFloat2D], list[ArrayFloat1D]]]:
         """Internal inplane corners and their convexity for each layer_refinement_specs."""
         cached_data = []
-        for layer in self.grid_spec.layer_refinement_specs:
+        for merged_geos, layer in zip(
+            self._internal_layerrefinement_merged_geos, self.grid_spec.layer_refinement_specs
+        ):
             cached_data.append(
                 layer._corners_and_convexity_2d(
-                    structure_list=self.scene.all_structures, ravel=False
+                    merged_geos=merged_geos,
+                    structure_list=self.scene.all_structures,
+                    ravel=False,
+                    sim_bounds=self.bounds,
+                    boundary_type=self._internal_layerrefinement_boundary_types,
                 )
             )
         return cached_data
@@ -1005,9 +1038,11 @@ def internal_override_structures(self) -> list[MeshOverrideStructure]:
         return self.grid_spec.internal_override_structures(
             self.scene.all_structures,
             wavelength,
-            self.geometry.size,
+            self.bounds,
             self.lumped_elements,
+            self._internal_layerrefinement_boundary_types,
             self._internal_layerfinement_corners_and_convexity_2d,
+            self._internal_layerrefinement_merged_geos,
         )
 
     @cached_property
@@ -1022,7 +1057,10 @@ def internal_snapping_points(self) -> list[CoordinateOptional]:
         return self.grid_spec.internal_snapping_points(
             self.scene.all_structures,
             self.lumped_elements,
+            self._internal_layerrefinement_boundary_types,
+            self.bounds,
             self._internal_layerfinement_corners_and_convexity_2d,
+            self._internal_layerrefinement_merged_geos,
         )
 
     @equal_aspect
@@ -1389,6 +1427,7 @@ def _grid_and_snapping_lines(self) -> tuple[Grid, list[CoordinateOptional]]:
             internal_override_structures=self.internal_override_structures,
             boundary_types=boundary_types,
             structure_priority_mode=self.scene.structure_priority_mode,
+            cached_merged_geos=self._internal_layerrefinement_merged_geos,
         )
 
         # This would AutoGrid the in-plane directions of the 2D materials
@@ -5422,6 +5461,7 @@ def _grid_corrections_2dmaterials(self, grid: Grid) -> Grid:
             lumped_elements=self.lumped_elements,
             internal_snapping_points=self.internal_snapping_points,
             internal_override_structures=self.internal_override_structures,
+            cached_merged_geos=self._internal_layerrefinement_merged_geos,
         )
 
         # Handle 2D materials if ``AutoGrid`` is used for in-plane directions