diff --git a/tests/test_meshgenerate.py b/tests/test_meshgenerate.py
index 20df9224fe..cf3db601b4 100644
--- a/tests/test_meshgenerate.py
+++ b/tests/test_meshgenerate.py
@@ -329,11 +329,7 @@ def test_grid_refinement():
 def test_parse_structures():
     """Test some aspects of the structure parsing."""
 
-    source = td.PointDipole(
-        source_time=td.GaussianPulse(freq0=1e14, fwidth=1e13),
-        size=(0, 0, 0),
-        polarization="Ex",
-    )
+    wavelength = 2.9
 
     box1 = td.Structure(
         geometry=td.Box(center=(0, 0, 0), size=(2, 2, 2)), medium=td.Medium(permittivity=9)
@@ -354,10 +350,9 @@ def test_parse_structures():
     # Test that the box2 permittivity is used along z in the region where it fully covers box1
     sim = td.Simulation(
         size=(3, 3, 3),
-        grid_spec=td.GridSpec.auto(),
+        grid_spec=td.GridSpec.auto(wavelength=wavelength),
         run_time=1e-13,
         structures=[box1, box2],
-        sources=[source],
     )
     sizes = sim.grid.sizes.to_list[2]
     assert sizes[sizes.size // 2] > 0.1
@@ -365,10 +360,9 @@ def test_parse_structures():
     # Test that the box3 permittivity is not used along z as it doesn't fully cover box1
     sim = td.Simulation(
         size=(3, 3, 3),
-        grid_spec=td.GridSpec.auto(),
+        grid_spec=td.GridSpec.auto(wavelength=wavelength),
         run_time=1e-13,
         structures=[box1, box3],
-        sources=[source],
     )
     sizes = sim.grid.sizes.to_list[2]
     assert sizes[sizes.size // 2] < 0.1
@@ -378,10 +372,41 @@ def test_parse_structures():
     assert 1.0 in boundaries
     sim = td.Simulation(
         size=(3, 3, 3),
-        grid_spec=td.GridSpec.auto(),
+        grid_spec=td.GridSpec.auto(wavelength=wavelength),
         run_time=1e-13,
         structures=[box1, box4],
-        sources=[source],
     )
     boundaries = sim.grid.boundaries.to_list[1]
     assert 1.0 not in boundaries
+
+    # Test high-index background medium
+    sim = td.Simulation(
+        size=(3, 3, 6),
+        grid_spec=td.GridSpec.auto(wavelength=wavelength),
+        run_time=1e-13,
+        structures=[box1, box2],
+        medium=td.Medium(permittivity=5**2),
+    )
+    sizes = sim.grid.sizes.to_list[2]
+    assert sizes[0] < wavelength / 50
+
+    # Test high-index background with override box
+    sim = td.Simulation(
+        size=(3, 3, 6),
+        grid_spec=td.GridSpec.auto(
+            wavelength=wavelength,
+            override_structures=[
+                td.Structure(
+                    geometry=td.Box(size=(td.inf, td.inf, td.inf)),
+                    medium=td.Medium(permittivity=1),
+                ),
+                box1,
+                box2,
+            ],
+        ),
+        run_time=1e-13,
+        structures=[box1, box2],
+        medium=td.Medium(permittivity=5**2),
+    )
+    sizes = sim.grid.sizes.to_list[2]
+    assert np.isclose(sizes[0], wavelength / 10)
diff --git a/tidy3d/components/grid/mesher.py b/tidy3d/components/grid/mesher.py
index c7b4feab58..5e44668e05 100644
--- a/tidy3d/components/grid/mesher.py
+++ b/tidy3d/components/grid/mesher.py
@@ -99,9 +99,9 @@ def parse_structures(
         # List of indexes of structures which are present in each interval
         interval_structs = [[0]]  # will have len equal to len(interval_coords) - 1
         # List of indexes of structures that every structre contains in 2D
-        struct_contains = []  # will have len equal to len(structures)
+        struct_contains = [[]]  # will have len equal to len(structures)
 
-        for struct_ind in range(len(structures) - 1, 0, -1):
+        for struct_ind in range(1, len(structures)):
             bbox = struct_bbox[struct_ind]
 
             # indexes of structures that the current structure contains in 2D
@@ -135,11 +135,6 @@ def parse_structures(
             for interval_ind in range(indmin, indmax):
                 interval_structs[interval_ind].append(struct_ind)
 
-        # The simulation domain is the lowest structure and doesn't contain anything
-        struct_contains.append([])
-        # Reverse to match the order in the structures list
-        struct_contains = struct_contains[::-1]
-
         # Truncate intervals to domain bounds
         b_array = np.array(interval_coords)
         in_domain = np.argwhere((b_array >= domain_bounds[0]) * (b_array <= domain_bounds[1]))
@@ -164,7 +159,7 @@ def parse_structures(
             struct_list = interval_structs[coord_ind]
             struct_list_filter = []
             # Handle containment
-            for ind, struct_ind in enumerate(struct_list):
+            for ind, struct_ind in enumerate(struct_list[::-1]):
                 if ind >= len(struct_list):
                     # This can happen because we modify struct_list in the loop
                     break