diff --git a/tests/test_components/test_layerrefinement.py b/tests/test_components/test_layerrefinement.py
index dadc92602..a9a79e0bd 100644
--- a/tests/test_components/test_layerrefinement.py
+++ b/tests/test_components/test_layerrefinement.py
@@ -402,3 +402,89 @@ def test_dl_min_from_smallest_feature():
     )
 
     _ = sim.grid
+
+
+def test_gap_meshing():
+    w = 1
+    l = 10
+
+    l_shape_1 = td.Structure(
+        medium=td.PECMedium(),
+        geometry=td.PolySlab(
+            axis=2,
+            slab_bounds=[0, 2],
+            vertices=[
+                [0, 0],
+                [l, 0],
+                [l, w],
+                [w, w],
+                [w, l],
+                [0, l],
+            ],
+        ),
+    )
+
+    gap = 0.1
+    l_shape_2 = td.Structure(
+        medium=td.PECMedium(),
+        geometry=td.PolySlab(
+            axis=2,
+            slab_bounds=[0, 2],
+            vertices=[
+                [w + gap, w + gap],
+                [w + gap + l, w + gap],
+                [w + gap + l, w + gap + w],
+                [w + gap + w, w + gap + w],
+                [w + gap + w, w + gap + l],
+                [w + gap, w + gap + l],
+                [0, w + gap + l],
+                [0, gap + l],
+                [w + gap, gap + l],
+            ],
+        ),
+    )
+
+    # ax = l_shape_1.plot(z=1)
+    # l_shape_2.plot(z=1, ax=ax)
+
+    for num_iters in range(2):
+        grid_spec = td.GridSpec(
+            grid_x=td.AutoGrid(min_steps_per_wvl=10),
+            grid_y=td.AutoGrid(min_steps_per_wvl=10),
+            grid_z=td.AutoGrid(min_steps_per_wvl=10),
+            layer_refinement_specs=[
+                td.LayerRefinementSpec(
+                    axis=2,
+                    corner_finder=None,
+                    gap_meshing_iters=num_iters,
+                    size=[td.inf, td.inf, 2],
+                    center=[0, 0, 1],
+                )
+            ],
+            wavelength=7,
+        )
+
+        sim = td.Simulation(
+            structures=[l_shape_1, l_shape_2],
+            grid_spec=grid_spec,
+            size=(1.2 * l, 1.2 * l, 2),
+            center=(0.5 * l, 0.5 * l, 0),
+            run_time=1e-15,
+        )
+
+        # ax = sim.plot(z=1)
+        # sim.plot_grid(z=1, ax=ax)
+        # ax.set_xlim([0, 2])
+        # ax.set_ylim([0, 2])
+
+        resolved_x = np.any(
+            np.logical_and(sim.grid.boundaries.x > w, sim.grid.boundaries.x < w + gap)
+        )
+        resolved_y = np.any(
+            np.logical_and(sim.grid.boundaries.y > w, sim.grid.boundaries.y < w + gap)
+        )
+
+        if num_iters == 0:
+            assert (not resolved_x) and (not resolved_y)
+        else:
+            assert resolved_x and resolved_y
diff --git a/tidy3d/components/grid/corner_finder.py b/tidy3d/components/grid/corner_finder.py
index 2e61d4fd0..7b3d7448a 100644
--- a/tidy3d/components/grid/corner_finder.py
+++ b/tidy3d/components/grid/corner_finder.py
@@ -75,16 +75,14 @@ def _no_min_dl_override(self):
             )
         )
 
-    def _corners_and_convexity(
-        self,
+    @classmethod
+    def _merged_pec_on_plane(
+        cls,
         normal_axis: Axis,
         coord: float,
         structure_list: List[Structure],
-        ravel: bool,
     ) -> Tuple[ArrayFloat2D, ArrayFloat1D]:
-        """On a 2D plane specified by axis = `normal_axis` and coordinate `coord`, find out corners of merged
-        geometries made of `medium`.
-
+        """On a 2D plane specified by axis = `normal_axis` and coordinate `coord`, merge geometries made of PEC.
 
         Parameters
         ----------
@@ -94,8 +92,6 @@ def _corners_and_convexity(
             Position of plane along the normal axis.
         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.
 
         Returns
         -------
@@ -121,6 +117,41 @@ def _corners_and_convexity(
         # merge geometries
         merged_geos = merging_geometries_on_plane(geometry_list, plane, medium_list)
 
+        return merged_geos
+
+    def _corners_and_convexity(
+        self,
+        normal_axis: Axis,
+        coord: float,
+        structure_list: List[Structure],
+        ravel: bool,
+    ) -> Tuple[ArrayFloat2D, ArrayFloat1D]:
+        """On a 2D plane specified by axis = `normal_axis` and coordinate `coord`, find out corners of merged
+        geometries made of PEC.
+
+
+        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.
+        ravel : bool
+            Whether to put the resulting corners in a single list or per polygon.
+
+        Returns
+        -------
+        ArrayFloat2D
+            Corner coordinates.
+        """
+
+        # merge geometries
+        merged_geos = self._merged_pec_on_plane(
+            normal_axis=normal_axis, coord=coord, structure_list=structure_list
+        )
+
         # corner finder
         corner_list = []
         convexity_list = []
diff --git a/tidy3d/components/grid/grid_spec.py b/tidy3d/components/grid/grid_spec.py
index 8942a8c34..2c99d8a6a 100644
--- a/tidy3d/components/grid/grid_spec.py
+++ b/tidy3d/components/grid/grid_spec.py
@@ -12,7 +12,7 @@
 from ...exceptions import SetupError
 from ...log import log
 from ..base import Tidy3dBaseModel, cached_property, skip_if_fields_missing
-from ..geometry.base import Box
+from ..geometry.base import Box, ClipOperation
 from ..lumped_element import LumpedElementType
 from ..source.utils import SourceType
 from ..structure import MeshOverrideStructure, Structure, StructureType
@@ -1024,6 +1024,19 @@ class LayerRefinementSpec(Box):
         "and the projection of the simulation domain overlaps.",
     )
 
+    gap_meshing_iters: pd.NonNegativeInt = pd.Field(
+        1,
+        title="Gap Meshing Iterations",
+        description="Number of recursive iterations for resolving thin gaps. "
+        "The underlying algorithm detects gaps contained in a single cell and places a snapping plane at the gaps's centers.",
+    )
+
+    dl_min_from_gap_width: bool = pd.Field(
+        True,
+        title="Set ``dl_min`` from Gap Width",
+        description="Take into account autodetected minimal PEC gap width when determining ``dl_min``.",
+    )
+
     @pd.validator("axis", always=True)
     @skip_if_fields_missing(["size"])
     def _finite_size_along_axis(cls, val, values):
@@ -1044,6 +1057,8 @@ def from_layer_bounds(
         corner_snapping: bool = True,
         corner_refinement: GridRefinement = GridRefinement(),
         refinement_inside_sim_only: bool = True,
+        gap_meshing_iters: pd.NonNegativeInt = 1,
+        dl_min_from_gap_width: bool = True,
     ):
         """Constructs a :class:`LayerRefiementSpec` that is unbounded in inplane dimensions from bounds along
         layer thickness dimension.
@@ -1069,6 +1084,10 @@ def from_layer_bounds(
             Inplane mesh refinement factor around corners.
         refinement_inside_sim_only : bool = True
             Apply refinement only to features inside simulation domain.
+        gap_meshing_iters : bool = True
+            Number of recursive iterations for resolving thin gaps.
+        dl_min_from_gap_width : bool = True
+            Take into account autodetected minimal PEC gap width when determining ``dl_min``.
 
 
         Example
@@ -1090,6 +1109,8 @@ def from_layer_bounds(
             corner_snapping=corner_snapping,
             corner_refinement=corner_refinement,
             refinement_inside_sim_only=refinement_inside_sim_only,
+            gap_meshing_iters=gap_meshing_iters,
+            dl_min_from_gap_width=dl_min_from_gap_width,
         )
 
     @classmethod
@@ -1105,6 +1126,8 @@ def from_bounds(
         corner_snapping: bool = True,
         corner_refinement: GridRefinement = GridRefinement(),
         refinement_inside_sim_only: bool = True,
+        gap_meshing_iters: pd.NonNegativeInt = 1,
+        dl_min_from_gap_width: bool = True,
     ):
         """Constructs a :class:`LayerRefiementSpec` from minimum and maximum coordinate bounds.
 
@@ -1132,6 +1155,10 @@ def from_bounds(
             Inplane mesh refinement factor around corners.
         refinement_inside_sim_only : bool = True
             Apply refinement only to features inside simulation domain.
+        gap_meshing_iters : bool = True
+            Number of recursive iterations for resolving thin gaps.
+        dl_min_from_gap_width : bool = True
+            Take into account autodetected minimal PEC gap width when determining ``dl_min``.
 
 
         Example
@@ -1153,6 +1180,8 @@ def from_bounds(
             corner_snapping=corner_snapping,
             corner_refinement=corner_refinement,
             refinement_inside_sim_only=refinement_inside_sim_only,
+            gap_meshing_iters=gap_meshing_iters,
+            dl_min_from_gap_width=dl_min_from_gap_width,
         )
 
     @classmethod
@@ -1167,6 +1196,8 @@ def from_structures(
         corner_snapping: bool = True,
         corner_refinement: GridRefinement = GridRefinement(),
         refinement_inside_sim_only: bool = True,
+        gap_meshing_iters: pd.NonNegativeInt = 1,
+        dl_min_from_gap_width: bool = True,
     ):
         """Constructs a :class:`LayerRefiementSpec` from the bounding box of a list of structures.
 
@@ -1192,6 +1223,10 @@ def from_structures(
             Inplane mesh refinement factor around corners.
         refinement_inside_sim_only : bool = True
             Apply refinement only to features inside simulation domain.
+        gap_meshing_iters : bool = True
+            Number of recursive iterations for resolving thin gaps.
+        dl_min_from_gap_width : bool = True
+            Take into account autodetected minimal PEC gap width when determining ``dl_min``.
 
         """
 
@@ -1213,6 +1248,8 @@ def from_structures(
             corner_snapping=corner_snapping,
             corner_refinement=corner_refinement,
             refinement_inside_sim_only=refinement_inside_sim_only,
+            gap_meshing_iters=gap_meshing_iters,
+            dl_min_from_gap_width=dl_min_from_gap_width,
         )
 
     @cached_property
@@ -1497,6 +1534,253 @@ def _override_structures_along_axis(
                 override_structures += refinement_structures
         return override_structures
 
+    def _find_vertical_intersections(self, x, y, vertices):
+        """Detect intersection points of single polygon and vertical grid lines."""
+
+        cells_ij = []
+        cells_dy = []
+
+        # for each polygon vertex find the index of the first grid line on the right
+        grid_lines_on_right = np.argmax(x[:, None] > vertices[None, :, 0], axis=0)
+        grid_lines_on_right[vertices[:, 0] > x[-1]] = len(x)
+        # once we know these indices then we can find grid lines intersected by the i-th
+        # segment of the polygon as
+        # [grid_lines_on_right[i], grid_lines_on_right[i+1]) for grid_lines_on_right[i] > grid_lines_on_right[i+1]
+        # or
+        # [grid_lines_on_right[i+1], grid_lines_on_right[i]) for grid_lines_on_right[i] < grid_lines_on_right[i+1]
+
+        # loop over segments of the polygon and determine in which cells and where exactly they cross grid lines
+        # v_beg and v_end are the starting and ending points of the segment
+        # ind_beg and ind_end are starting and ending indices of vertical grid lines that the segment intersects
+        # as described above
+        for ind_beg, ind_end, v_beg, v_end in zip(
+            grid_lines_on_right,
+            np.roll(grid_lines_on_right, -1),
+            vertices,
+            np.roll(vertices, axis=0, shift=-1),
+        ):
+            # no intersections
+            if ind_end == ind_beg:
+                continue
+
+            # sort vertices in ascending order to make treatmeant unifrom
+            if ind_beg > ind_end:
+                ind_beg, ind_end, v_beg, v_end = ind_end, ind_beg, v_end, v_beg
+
+            # x coordinates are simply x coordinates of intersected vertical grid lines
+            intersections_x = x[ind_beg:ind_end]
+
+            # y coordinates can be found from line equation
+            intersections_y = v_beg[1] + (v_end[1] - v_beg[1]) / (v_end[0] - v_beg[0]) * (
+                intersections_x - v_beg[0]
+            )
+
+            # however, some of the vertical lines might be crossed
+            # outside of computational domain
+            # so we need to see which ones are actually inside along y axis
+            inds_inside_grid = np.logical_and(intersections_y >= y[0], intersections_y <= y[-1])
+
+            intersections_y = intersections_y[inds_inside_grid]
+
+            # find i and j indices of cells which contain these intersections
+
+            # i indices are simply indices of crossed vertical grid lines
+            cell_is = np.arange(ind_beg, ind_end)[inds_inside_grid]
+
+            # j indices can be computed by finding insertion indices
+            # of y coordinates of intersection points into array of y coordinates
+            # of the grid lines that preserve sorting
+            cell_js = np.searchsorted(y, intersections_y) - 1
+
+            # find local dy, that is, the distance between the intersection point
+            # and the bottom edge of the cell
+            dy = intersections_y - y[cell_js]
+
+            # record info
+            cells_ij.append(np.transpose([cell_is, cell_js]))
+            cells_dy.append(dy)
+
+        if len(cells_ij) > 0:
+            cells_ij = np.concatenate(cells_ij)
+
+        if len(cells_dy) > 0:
+            cells_dy = np.concatenate(cells_dy)
+
+        return cells_ij, cells_dy
+
+    def _process_poly(self, x, y, vertices):
+        """Detect intersection points of single polygon and grid lines."""
+
+        # find cells that contain intersections of vertical grid lines
+        # and locations of those intersections (along y axis)
+        v_cells_ij, v_cells_dy = self._find_vertical_intersections(x, y, vertices)
+
+        # find cells that contain intersections of horizontal grid lines
+        # and locations of those intersections (along x axis)
+        # reuse the same command but flip dimensions
+        h_cells_ij, h_cells_dx = self._find_vertical_intersections(y, x, np.flip(vertices, axis=1))
+        if len(h_cells_ij) > 0:
+            h_cells_ij = np.roll(h_cells_ij, axis=1, shift=1)
+
+        return v_cells_ij, v_cells_dy, h_cells_ij, h_cells_dx
+
+    def _process_slice(self, x, y, merged_geos):
+        """Detect intersection points of geometries boundaries and grid lines."""
+
+        # cells that contain intersections of vertical grid lines
+        v_cells_ij = []
+        # locations of those intersections (along y axis)
+        v_cells_dy = []
+
+        # cells that contain intersections of horizontal grid lines
+        h_cells_ij = []
+        # locations of those intersections (along x axis)
+        h_cells_dx = []
+
+        # loop over all shapes
+        for mat, shapes in merged_geos:
+            if not mat.is_pec:
+                continue
+            polygon_list = ClipOperation.to_polygon_list(shapes)
+            for poly in polygon_list:
+                poly = poly.normalize().buffer(0)
+
+                # find intersections of a polygon with grid lines
+                # specifically:
+                # 0. cells that contain intersections of vertical grid lines
+                # 1. locations of those intersections along y axis
+                # 2. cells that contain intersections of horizontal grid lines
+                # 3. locations of those intersections along x axis
+                data = self._process_poly(x, y, np.array(poly.exterior.coords)[:-1])
+
+                if len(data[0]) > 0:
+                    v_cells_ij.append(data[0])
+                    v_cells_dy.append(data[1])
+
+                if len(data[2]) > 0:
+                    h_cells_ij.append(data[2])
+                    h_cells_dx.append(data[3])
+
+                # in case the polygon has holes
+                for poly_inner in poly.interiors:
+                    data = self._process_poly(x, y, np.array(poly_inner.coords)[:-1])
+                    if len(data[0]) > 0:
+                        v_cells_ij.append(data[0])
+                        v_cells_dy.append(data[1])
+
+                    if len(data[2]) > 0:
+                        h_cells_ij.append(data[2])
+                        h_cells_dx.append(data[3])
+
+        if len(v_cells_ij) > 0:
+            v_cells_ij = np.concatenate(v_cells_ij)
+
+        if len(v_cells_dy) > 0:
+            v_cells_dy = np.concatenate(v_cells_dy)
+
+        if len(h_cells_ij) > 0:
+            h_cells_ij = np.concatenate(h_cells_ij)
+
+        if len(h_cells_dx) > 0:
+            h_cells_dx = np.concatenate(h_cells_dx)
+
+        return v_cells_ij, v_cells_dy, h_cells_ij, h_cells_dx
+
+    def _resolve_gaps(self, structures: List, grid: Grid):
+        """Detect underresolved gaps and place snapping lines in them. Also return the detected minimal gap width."""
+
+        # get merged pec structures on plane
+        plane_slice = CornerFinderSpec._merged_pec_on_plane(
+            coord=self.center_axis, normal_axis=self.axis, structure_list=structures
+        )
+
+        # get x and y coordinates of grid lines
+        _, tan_dims = Box.pop_axis([0, 1, 2], self.axis)
+        x = grid.boundaries.to_list[tan_dims[0]]
+        y = grid.boundaries.to_list[tan_dims[1]]
+
+        # find intersections of pec polygons with grid lines
+        # specifically:
+        # 0. cells that contain intersections of vertical grid lines
+        # 1. locations of those intersections along y axis
+        # 2. cells that contain intersections of horizontal grid lines
+        # 3. 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)
+
+        detected_gap_width = inf
+
+        # generate horizontal snapping lines
+        snapping_lines_y = []
+        if len(v_cells_ij) > 0:
+            # count intersections of vertical grid lines in each cell
+            ind_unique, counts = np.unique(
+                v_cells_ij[:, 0] * len(y) + v_cells_ij[:, 1], return_counts=True
+            )
+            # when we count intersections we use linearized 2d index because we really
+            # need to count intersections in each cell separately
+
+            # but when we need to decide about refinement, due to cartesian nature of grid
+            # we will need to consider all cells with a given j index at a time
+
+            # so, let's compute j index for each cell in the inique list
+            ind_unique_j = ind_unique % len(y)
+
+            # loop through all j rows that contain intersections
+            for ind_j in np.unique(ind_unique_j):
+                # we need to refine between two grid lines corresponding to index j
+                # if at least one cell with given j contains > 1 intersections
+                max_count = np.max(counts[ind_unique_j == ind_j])
+                if max_count > 1:
+                    # This could use some improvement in future:
+                    # currently we just calculate the position of the lowest intersection
+                    # and the highest intersection in row j
+                    # and just place a snapping line in between.
+                    min_dy = np.min(v_cells_dy[v_cells_ij[:, 1] == ind_j])
+                    max_dy = np.max(v_cells_dy[v_cells_ij[:, 1] == ind_j])
+                    snapping_lines_y.append(y[ind_j] + 0.5 * (min_dy + max_dy))
+                    detected_gap_width = min(detected_gap_width, max_dy - min_dy)
+                    # This is expected to work fine when we have a simple gap/strip passing
+                    # through a number of cells. But not optimal if we have slanted/curved
+                    # boundaries and/or more than 2 intersections.
+
+        # generate vertical snapping lines
+        snapping_lines_x = []
+        if len(h_cells_ij) > 0:
+            # count intersections of horizontal grid lines in each cell
+            ind_unique, counts = np.unique(
+                h_cells_ij[:, 0] * len(y) + h_cells_ij[:, 1], return_counts=True
+            )
+            # when we count intersections we use linearized 2d index because we really
+            # need to count intersections in each cell separately
+
+            # but when we need to decide about refinement, due to cartesian nature of grid
+            # we will need to consider all cells with a given i index at a time
+
+            # so, let's compute i index for each cell in the inique list
+            ind_unique_i = ind_unique // len(y)
+
+            # loop through all i columns that contain intersections
+            for ind_i in np.unique(ind_unique_i):
+                # we need to refine between two grid lines corresponding to index i
+                # if at least one cell with given i contains > 1 intersections
+                max_count = np.max(counts[ind_unique_i == ind_i])
+                if max_count > 1:
+                    # This could use some improvement in future:
+                    # currently we just calculate the position of the leftmost intersection
+                    # and the rightmost intersection in column i
+                    # and just place a snapping line in between.
+                    min_dx = np.min(h_cells_dx[h_cells_ij[:, 0] == ind_i])
+                    max_dx = np.max(h_cells_dx[h_cells_ij[:, 0] == ind_i])
+                    snapping_lines_x.append(x[ind_i] + 0.5 * (min_dx + max_dx))
+                    detected_gap_width = min(detected_gap_width, max_dx - min_dx)
+                    # This is expected to work fine when we have a simple gap/strip passing
+                    # through a number of cells. But not optimal if we have slanted/curved
+                    # boundaries and/or more than 2 intersections.
+
+        return [Box.unpop_axis(X, (None, None), axis=tan_dims[0]) for X in snapping_lines_x] + [
+            Box.unpop_axis(Y, (None, None), axis=tan_dims[1]) for Y in snapping_lines_y
+        ], detected_gap_width
+
 
 class GridSpec(Tidy3dBaseModel):
     """Collective grid specification for all three dimensions.
@@ -1921,6 +2205,108 @@ def make_grid(
             Entire simulation grid.
         """
 
+        old_grid = self._make_grid_one_iteration(
+            structures=structures,
+            symmetry=symmetry,
+            periodic=periodic,
+            sources=sources,
+            num_pml_layers=num_pml_layers,
+            lumped_elements=lumped_elements,
+            internal_override_structures=internal_override_structures,
+            internal_snapping_points=internal_snapping_points,
+        )
+
+        if len(self.layer_refinement_specs) > 0:
+            num_iters = max(
+                layer_spec.gap_meshing_iters for layer_spec in self.layer_refinement_specs
+            )
+
+            snapping_lines = []
+            min_gap_width = inf
+            for ind in range(num_iters):
+                new_snapping_lines = []
+                for layer_spec in self.layer_refinement_specs:
+                    if layer_spec.gap_meshing_iters > ind:
+                        one_layer_snapping_lines, gap_width = layer_spec._resolve_gaps(
+                            structures, old_grid
+                        )
+                        new_snapping_lines = new_snapping_lines + one_layer_snapping_lines
+                        if layer_spec.dl_min_from_gap_width:
+                            min_gap_width = min(min_gap_width, gap_width)
+
+                if len(new_snapping_lines) == 0:
+                    log.info(
+                        "Grid is no longer changing. "
+                        f"Stopping iterative gap meshing after {ind+1}/{num_iters} iterations."
+                    )
+                    break
+
+                snapping_lines = snapping_lines + new_snapping_lines
+
+                new_grid = self._make_grid_one_iteration(
+                    structures=structures,
+                    symmetry=symmetry,
+                    periodic=periodic,
+                    sources=sources,
+                    num_pml_layers=num_pml_layers,
+                    lumped_elements=lumped_elements,
+                    internal_override_structures=internal_override_structures,
+                    internal_snapping_points=internal_snapping_points + snapping_lines,
+                    min_dl_from_gaps=0.4 * min_gap_width,
+                )
+
+                same = old_grid == new_grid
+
+                if same:
+                    log.info(
+                        "Grid is no longer changing. "
+                        f"Stopping iterative gap meshing after {ind+1}/{num_iters} iterations."
+                    )
+                    break
+
+                old_grid = new_grid
+
+        return old_grid
+
+    def _make_grid_one_iteration(
+        self,
+        structures: List[Structure],
+        symmetry: Tuple[Symmetry, Symmetry, Symmetry],
+        periodic: Tuple[bool, bool, bool],
+        sources: List[SourceType],
+        num_pml_layers: List[Tuple[pd.NonNegativeInt, pd.NonNegativeInt]],
+        lumped_elements: List[LumpedElementType] = (),
+        internal_override_structures: List[MeshOverrideStructure] = None,
+        internal_snapping_points: List[CoordinateOptional] = None,
+        min_dl_from_gaps: pd.PositiveFloat = inf,
+    ) -> Grid:
+        """Make the entire simulation grid based on some simulation parameters.
+
+        Parameters
+        ----------
+        structures : List[Structure]
+            List of structures present in the simulation. The first structure must be the
+            simulation geometry with the simulation background medium.
+        symmetry : Tuple[Symmetry, Symmetry, Symmetry]
+            Reflection symmetry across a plane bisecting the simulation domain
+            normal to each of the three axes.
+        sources : List[SourceType]
+            List of sources.
+        num_pml_layers : List[Tuple[float, float]]
+            List containing the number of absorber layers in - and + boundaries.
+        lumped_elements : List[LumpedElementType]
+            List of lumped elements.
+        internal_override_structures : List[MeshOverrideStructure]
+            If `None`, recomputes internal override structures.
+        internal_snapping_points : List[CoordinateOptional]
+            If `None`, recomputes internal snapping points.
+
+        Returns
+        -------
+        Grid:
+            Entire simulation grid.
+        """
+
         # Set up wavelength for automatic mesh generation if needed.
         wavelength = self.get_wavelength(sources)
 
@@ -1992,6 +2378,7 @@ def make_grid(
                 sim_size,
                 lumped_elements,
             )
+            new_dl_min = min(new_dl_min, min_dl_from_gaps)
             for ind, grid in enumerate(grids_1d):
                 if isinstance(grid, AutoGrid) and grid._undefined_dl_min:
                     grids_1d[ind] = grid.updated_copy(dl_min=new_dl_min)