--- a/src/Grids/manifolds.jl	Fri Feb 07 08:44:36 2025 +0100
+++ b/src/Grids/manifolds.jl	Fri Feb 07 09:12:41 2025 +0100
@@ -68,23 +68,20 @@
 function connections(a::CartesianAtlas)
     c = Tuple{MultiBlockBoundary, MultiBlockBoundary}[]
 
-    N,M = size(a.charts)
-    for j ∈ 1:M, i ∈ 1:N-1
-        push!(c,
-            (
-                MultiBlockBoundary{(i,j), CartesianBoundary{1,UpperBoundary}}(),
-                MultiBlockBoundary{(i+1,j), CartesianBoundary{1,LowerBoundary}}(),
-            ),
-        )
-    end
-
-    for i ∈ 1:N, j ∈ 1:M-1
-        push!(c,
-            (
-                MultiBlockBoundary{(i,j), CartesianBoundary{2,UpperBoundary}}(),
-                MultiBlockBoundary{(i,j+1), CartesianBoundary{2,LowerBoundary}}(),
-            ),
-        )
+    for d ∈ 1:ndims(charts(a))
+        Is = eachslice(CartesianIndices(charts(a)); dims=d)
+        for i ∈ 1:length(Is)-1 # For each interface between slices
+            for jk ∈ eachindex(Is[i]) # For each block in slice
+                Iᵢⱼₖ = Tuple(Is[i][jk])
+                Iᵢ₊₁ⱼₖ = Tuple(Is[i+1][jk])
+                push!(c,
+                    (
+                        MultiBlockBoundary{Iᵢⱼₖ,   CartesianBoundary{d,UpperBoundary}}(),
+                        MultiBlockBoundary{Iᵢ₊₁ⱼₖ, CartesianBoundary{d,LowerBoundary}}(),
+                    )
+                )
+            end
+        end
     end
 
     return c

--- a/test/Grids/manifolds_test.jl	Fri Feb 07 08:44:36 2025 +0100
+++ b/test/Grids/manifolds_test.jl	Fri Feb 07 09:12:41 2025 +0100
@@ -36,6 +36,7 @@
 
     @testset "connections" begin
         a = CartesianAtlas(fill(c, 2,3))
+        # 2D
         west = CartesianBoundary{1,LowerBoundary}
         east = CartesianBoundary{1,UpperBoundary}
         south = CartesianBoundary{2,LowerBoundary}
@@ -50,6 +51,37 @@
             (MultiBlockBoundary{(2,2), north}(), MultiBlockBoundary{(2,3), south}()),
             (MultiBlockBoundary{(1,3), east}(), MultiBlockBoundary{(2,3), west}()),
         ])
+
+        # 3D
+        a = CartesianAtlas(fill(c, 2,2,3))
+        bottom = CartesianBoundary{3, LowerBoundary}
+        top = CartesianBoundary{3, UpperBoundary}
+        @test Set(connections(a)) == Set([
+            (MultiBlockBoundary{(1,1,1), east}(),  MultiBlockBoundary{(2,1,1), west}()),
+            (MultiBlockBoundary{(1,1,1), north}(), MultiBlockBoundary{(1,2,1), south}()),
+            (MultiBlockBoundary{(2,1,1), north}(), MultiBlockBoundary{(2,2,1), south}()),
+            (MultiBlockBoundary{(1,2,1), east}(),  MultiBlockBoundary{(2,2,1), west}()),
+
+            (MultiBlockBoundary{(1,1,2), east}(),  MultiBlockBoundary{(2,1,2), west}()),
+            (MultiBlockBoundary{(1,1,2), north}(), MultiBlockBoundary{(1,2,2), south}()),
+            (MultiBlockBoundary{(2,1,2), north}(), MultiBlockBoundary{(2,2,2), south}()),
+            (MultiBlockBoundary{(1,2,2), east}(),  MultiBlockBoundary{(2,2,2), west}()),
+
+            (MultiBlockBoundary{(1,1,3), east}(),  MultiBlockBoundary{(2,1,3), west}()),
+            (MultiBlockBoundary{(1,1,3), north}(), MultiBlockBoundary{(1,2,3), south}()),
+            (MultiBlockBoundary{(2,1,3), north}(), MultiBlockBoundary{(2,2,3), south}()),
+            (MultiBlockBoundary{(1,2,3), east}(),  MultiBlockBoundary{(2,2,3), west}()),
+
+            (MultiBlockBoundary{(1,1,1), top}(), MultiBlockBoundary{(1,1,2), bottom}()),
+            (MultiBlockBoundary{(2,1,1), top}(), MultiBlockBoundary{(2,1,2), bottom}()),
+            (MultiBlockBoundary{(1,2,1), top}(), MultiBlockBoundary{(1,2,2), bottom}()),
+            (MultiBlockBoundary{(2,2,1), top}(), MultiBlockBoundary{(2,2,2), bottom}()),
+
+            (MultiBlockBoundary{(1,1,2), top}(), MultiBlockBoundary{(1,1,3), bottom}()),
+            (MultiBlockBoundary{(2,1,2), top}(), MultiBlockBoundary{(2,1,3), bottom}()),
+            (MultiBlockBoundary{(1,2,2), top}(), MultiBlockBoundary{(1,2,3), bottom}()),
+            (MultiBlockBoundary{(2,2,2), top}(), MultiBlockBoundary{(2,2,3), bottom}()),
+        ])
     end
 end