diff --git a/Project.toml b/Project.toml
index e3c5eb3fbd..f03c4aa344 100644
--- a/Project.toml
+++ b/Project.toml
@@ -54,6 +54,7 @@ SHA = "ea8e919c-243c-51af-8825-aaa63cd721ce"
 TaskLocalValues = "ed4db957-447d-4319-bfb6-7fa9ae7ecf34"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
+AllocCheck = "9b6a8646-10ed-4001-bbdc-1d2f46dfbb1a"
 
 [targets]
-test = ["BlockArrays", "Downloads", "FerriteGmsh", "ForwardDiff", "Gmsh", "IterativeSolvers", "Metis", "NBInclude", "OhMyThreads", "ProgressMeter", "Random", "SHA", "TaskLocalValues", "Test", "TimerOutputs", "Logging"]
+test = ["BlockArrays", "Downloads", "FerriteGmsh", "ForwardDiff", "Gmsh", "IterativeSolvers", "Metis", "NBInclude", "OhMyThreads", "ProgressMeter", "Random", "SHA", "TaskLocalValues", "Test", "TimerOutputs", "Logging", "AllocCheck"]
diff --git a/src/Dofs/DofHandler.jl b/src/Dofs/DofHandler.jl
index 84f59e62b7..3a8e7572e3 100644
--- a/src/Dofs/DofHandler.jl
+++ b/src/Dofs/DofHandler.jl
@@ -12,7 +12,7 @@ Access some grid representation for the dof handler.
 """
 get_grid(dh::AbstractDofHandler)
 
-mutable struct SubDofHandler{DH} <: AbstractDofHandler
+mutable struct SubDofHandler{DH, CT} <: AbstractDofHandler
     # From constructor
     const dh::DH
     const cellset::OrderedSet{Int}
@@ -68,7 +68,7 @@ function SubDofHandler(dh::DH, cellset::AbstractVecOrSet{Int}) where {DH <: Abst
         end
     end
     # Construct and insert into the parent dh
-    sdh = SubDofHandler{typeof(dh)}(dh, convert_to_orderedset(cellset), Symbol[], Interpolation[], Int[], -1)
+    sdh = SubDofHandler{typeof(dh), CT}(dh, convert_to_orderedset(cellset), Symbol[], Interpolation[], Int[], -1)
     push!(dh.subdofhandlers, sdh)
     return sdh
 end
diff --git a/src/Dofs/DofRenumbering.jl b/src/Dofs/DofRenumbering.jl
index 0d1a2f70b0..83e44615fb 100644
--- a/src/Dofs/DofRenumbering.jl
+++ b/src/Dofs/DofRenumbering.jl
@@ -196,7 +196,7 @@ function compute_renumber_permutation(dh::DofHandler, _, order::DofOrder.Compone
     for sdh in dh.subdofhandlers
         dof_ranges = [dof_range(sdh, f) for f in eachindex(sdh.field_names)]
         global_idxs = [findfirst(x -> x === f, dh.field_names) for f in sdh.field_names]
-        for cell in CellIterator(dh, sdh.cellset, flags)
+        for cell in CellIterator(sdh, flags)
             cdofs = celldofs(cell)
             for (local_idx, global_idx) in pairs(global_idxs)
                 rng = dof_ranges[local_idx]
diff --git a/src/Dofs/sparsity_pattern.jl b/src/Dofs/sparsity_pattern.jl
index 95f80841cc..bc1d14aef1 100644
--- a/src/Dofs/sparsity_pattern.jl
+++ b/src/Dofs/sparsity_pattern.jl
@@ -446,34 +446,39 @@ function _coupling_to_local_dof_coupling(dh::DofHandler, coupling::AbstractMatri
     sz == size(coupling, 2) || error("coupling not square")
 
     # Return one matrix per (potential) sub-domain
-    outs = Matrix{Bool}[]
+    nsdh = length(dh.subdofhandlers)
+    outs = Matrix{Matrix{Bool}}(undef, nsdh, nsdh)
     field_dims = map(fieldname -> n_components(dh, fieldname), dh.field_names)
 
-    for sdh in dh.subdofhandlers
-        out = zeros(Bool, ndofs_per_cell(sdh), ndofs_per_cell(sdh))
-        push!(outs, out)
+    for (sdh1_idx, sdh1) in enumerate(dh.subdofhandlers)
+        for (sdh2_idx, sdh2) in enumerate(dh.subdofhandlers)
+            out = zeros(Bool, ndofs_per_cell(sdh1), ndofs_per_cell(sdh2))
+            outs[sdh1_idx, sdh2_idx] = out
 
-        dof_ranges = [dof_range(sdh, f) for f in sdh.field_names]
-        global_idxs = [findfirst(x -> x === f, dh.field_names) for f in sdh.field_names]
+            dof_ranges_here = [dof_range(sdh1, f) for f in sdh1.field_names]
+            dof_ranges_there = [dof_range(sdh2, f) for f in sdh2.field_names]
+            global_idxs_here = [findfirst(x -> x === f, dh.field_names) for f in sdh1.field_names]
+            global_idxs_there = [findfirst(x -> x === f, dh.field_names) for f in sdh2.field_names]
 
-        if sz == length(dh.field_names) # Coupling given by fields
-            for (j, jrange) in pairs(dof_ranges), (i, irange) in pairs(dof_ranges)
-                out[irange, jrange] .= coupling[global_idxs[i], global_idxs[j]]
-            end
-        elseif sz == sum(field_dims) # Coupling given by components
-            component_offsets = pushfirst!(cumsum(field_dims), 0)
-            for (jf, jrange) in pairs(dof_ranges), (j, J) in pairs(jrange)
-                jc = mod1(j, field_dims[global_idxs[jf]]) + component_offsets[global_idxs[jf]]
-                for (i_f, irange) in pairs(dof_ranges), (i, I) in pairs(irange)
-                    ic = mod1(i, field_dims[global_idxs[i_f]]) + component_offsets[global_idxs[i_f]]
-                    out[I, J] = coupling[ic, jc]
+            if sz == length(dh.field_names) # Coupling given by fields
+                for (j, jrange) in pairs(dof_ranges_there), (i, irange) in pairs(dof_ranges_here)
+                    out[irange, jrange] .= coupling[global_idxs_here[i], global_idxs_there[j]]
+                end
+            elseif sz == sum(field_dims) # Coupling given by components
+                component_offsets = pushfirst!(cumsum(field_dims), 0)
+                for (jf, jrange) in pairs(dof_ranges_there), (j, J) in pairs(jrange)
+                    jc = mod1(j, field_dims[global_idxs_there[jf]]) + component_offsets[global_idxs_there[jf]]
+                    for (i_f, irange) in pairs(dof_ranges_here), (i, I) in pairs(irange)
+                        ic = mod1(i, field_dims[global_idxs_here[i_f]]) + component_offsets[global_idxs_here[i_f]]
+                        out[I, J] = coupling[ic, jc]
+                    end
                 end
+            elseif sz == ndofs_per_cell(sdh1) # Coupling given by template local matrix
+                # TODO: coupling[fieldhandler_idx] if different template per subddomain
+                out .= coupling
+            else
+                error("could not create coupling")
             end
-        elseif sz == ndofs_per_cell(sdh) # Coupling given by template local matrix
-            # TODO: coupling[fieldhandler_idx] if different template per subddomain
-            out .= coupling
-        else
-            error("could not create coupling")
         end
     end
     return outs
@@ -481,16 +486,13 @@ end
 
 function _add_cell_entries!(
         sp::AbstractSparsityPattern, dh::DofHandler, ch::Union{ConstraintHandler, Nothing},
-        keep_constrained::Bool, coupling::Union{Vector{<:AbstractMatrix{Bool}}, Nothing},
+        keep_constrained::Bool, coupling::Union{<:AbstractMatrix{<:AbstractMatrix{Bool}}, Nothing},
     )
     # Add all connections between dofs for every cell while filtering based
     # on a) constraints, and b) field/dof coupling.
-    cc = CellCache(dh)
     for (sdhi, sdh) in pairs(dh.subdofhandlers)
-        set = BitSet(sdh.cellset)
         coupling === nothing || (coupling_sdh = coupling[sdhi])
-        for cell_id in set
-            reinit!(cc, cell_id)
+        for cc in CellIterator(sdh)
             for (i, row) in pairs(cc.dofs)
                 # a) check constraint for row
                 !keep_constrained && haskey(ch.dofmapping, row) && continue
@@ -606,41 +608,39 @@ function _add_interface_entries!(
         interface_coupling::AbstractMatrix{Bool},
     )
     couplings = _coupling_to_local_dof_coupling(dh, interface_coupling)
-    for ic in InterfaceIterator(dh, topology)
-        # TODO: This looks like it can be optimized for the common case where
-        #       the cells are in the same subdofhandler
-        sdhs_idx = dh.cell_to_subdofhandler[cellid.([ic.a, ic.b])]
-        sdhs = dh.subdofhandlers[sdhs_idx]
-        to_check = Dict{Int, Vector{Int}}()
-        for (i, sdh) in pairs(sdhs)
-            sdh_idx = sdhs_idx[i]
-            coupling_sdh = couplings[sdh_idx]
-            for cell_field in sdh.field_names
-                dofrange1 = dof_range(sdh, cell_field)
-                cell_dofs = celldofs(sdh_idx == 1 ? ic.a : ic.b)
-                cell_field_dofs = @view cell_dofs[dofrange1]
-                for neighbor_field in sdh.field_names
-                    sdh2 = sdhs[i == 1 ? 2 : 1]
-                    neighbor_field ∈ sdh2.field_names || continue
-                    dofrange2 = dof_range(sdh2, neighbor_field)
-                    neighbor_dofs = celldofs(sdh_idx == 2 ? ic.a : ic.b)
-                    neighbor_field_dofs = @view neighbor_dofs[dofrange2]
-
-                    empty!(to_check)
-                    for (j, dof_j) in enumerate(dofrange2)
-                        for (i, dof_i) in enumerate(dofrange1)
-                            coupling_sdh[dof_i, dof_j] || continue
-                            push!(get!(Vector{Int}, to_check, j), i)
+    for (_i, sdh1) in enumerate(dh.subdofhandlers)
+        for (_j, sdh2) in enumerate(dh.subdofhandlers)
+            for ic in InterfaceIterator(sdh1, sdh2, topology, false)
+                # TODO: This looks like it can be optimized for the common case where
+                #       the cells are in the same subdofhandler
+                to_check = Dict{Int, Vector{Int}}()
+                coupling_sdh = couplings[_i, _j]
+                for cell_field in sdh1.field_names
+                    dofrange1 = dof_range(sdh1, cell_field)
+                    cell_dofs = celldofs(ic.a)
+                    cell_field_dofs = @view cell_dofs[dofrange1]
+                    for neighbor_field in sdh1.field_names
+                        neighbor_field ∈ sdh2.field_names || continue
+                        dofrange2 = dof_range(sdh2, neighbor_field)
+                        neighbor_dofs = celldofs(ic.b)
+                        neighbor_field_dofs = @view neighbor_dofs[dofrange2]
+
+                        empty!(to_check)
+                        for (j, dof_j) in enumerate(dofrange2)
+                            for (i, dof_i) in enumerate(dofrange1)
+                                coupling_sdh[dof_i, dof_j] || continue
+                                push!(get!(Vector{Int}, to_check, j), i)
+                            end
                         end
-                    end
 
-                    for (j, is) in to_check
-                        # Avoid coupling the shared dof in continuous interpolations as cross-element. They're coupled in the local coupling matrix.
-                        neighbor_field_dofs[j] ∈ cell_dofs && continue
-                        for i in is
-                            cell_field_dofs[i] ∈ neighbor_dofs && continue
-                            _add_interface_entry(sp, cell_field_dofs, neighbor_field_dofs, i, j, keep_constrained, ch)
-                            _add_interface_entry(sp, neighbor_field_dofs, cell_field_dofs, j, i, keep_constrained, ch)
+                        for (j, is) in to_check
+                            # Avoid coupling the shared dof in continuous interpolations as cross-element. They're coupled in the local coupling matrix.
+                            neighbor_field_dofs[j] ∈ cell_dofs && continue
+                            for i in is
+                                cell_field_dofs[i] ∈ neighbor_dofs && continue
+                                _add_interface_entry(sp, cell_field_dofs, neighbor_field_dofs, i, j, keep_constrained, ch)
+                                _add_interface_entry(sp, neighbor_field_dofs, cell_field_dofs, j, i, keep_constrained, ch)
+                            end
                         end
                     end
                 end
diff --git a/src/Ferrite.jl b/src/Ferrite.jl
index 7b358463be..3b220ac5c1 100644
--- a/src/Ferrite.jl
+++ b/src/Ferrite.jl
@@ -107,6 +107,13 @@ struct FacetIndex <: BoundaryIndex
     idx::Tuple{Int, Int} # cell and side
 end
 
+"""
+An `InterfaceIndex` wraps an (Int, Int, Int, Int) and defines an interface by pointing to a (cell_here, facet_here, cell_there, facet_there).
+"""
+struct InterfaceIndex <: BoundaryIndex
+    idx::Tuple{Int, Int, Int, Int} # cell - side - cell - side
+end
+
 const AbstractVecOrSet{T} = Union{AbstractSet{T}, AbstractVector{T}}
 const IntegerCollection = AbstractVecOrSet{<:Integer}
 
diff --git a/src/Grid/grid.jl b/src/Grid/grid.jl
index 3d27431bd1..de38b68992 100644
--- a/src/Grid/grid.jl
+++ b/src/Grid/grid.jl
@@ -641,20 +641,20 @@ boundaryfunction(::Type{EdgeIndex}) = edges
 boundaryfunction(::Type{VertexIndex}) = vertices
 boundaryfunction(::Type{FacetIndex}) = facets
 
-for INDEX in (:VertexIndex, :EdgeIndex, :FaceIndex, :FacetIndex)
+for INDEX in (:VertexIndex, :EdgeIndex, :FaceIndex, :FacetIndex, :InterfaceIndex)
     @eval begin
         #Constructor
-        ($INDEX)(a::Int, b::Int) = ($INDEX)((a, b))
+        ($INDEX)(args...) = ($INDEX)((args...,))
 
         Base.getindex(I::($INDEX), i::Int) = I.idx[i]
 
         #To be able to do a,b = faceidx
-        Base.iterate(I::($INDEX), state::Int = 1) = (state == 3) ? nothing : (I[state], state + 1)
+        Base.iterate(I::($INDEX), state::Int = 1) = (state == length(I.idx) + 1) ? nothing : (I[state], state + 1)
 
         # Necessary to check if, e.g. `(cellid, faceidx) in faceset`
         Base.isequal(x::$INDEX, y::$INDEX) = x.idx == y.idx
-        Base.isequal(x::Tuple{Int, Int}, y::$INDEX) = x[1] == y.idx[1] && x[2] == y.idx[2]
-        Base.isequal(y::$INDEX, x::Tuple{Int, Int}) = x[1] == y.idx[1] && x[2] == y.idx[2]
+        Base.isequal(x::NTuple{N, Int}, y::$INDEX) where {N} = all(i -> x[i] == y.idx[i], 1:N)
+        Base.isequal(y::$INDEX, x::NTuple{N, Int}) where {N} = all(i -> x[i] == y.idx[i], 1:N)
         Base.hash(x::$INDEX, h::UInt) = hash(x.idx, h)
     end
 end
diff --git a/src/Grid/utils.jl b/src/Grid/utils.jl
index ce5937010e..04400ae765 100644
--- a/src/Grid/utils.jl
+++ b/src/Grid/utils.jl
@@ -165,6 +165,7 @@ function _create_boundaryset(f::Function, grid::AbstractGrid, top::ExclusiveTopo
             cell_idx = ff_nh_idx[1]
             facet_nr = ff_nh_idx[2]
             cell = getcells(grid, cell_idx)
+            length(facets(cell)) < facet_nr && continue
             facet_nodes = facets(cell)[facet_nr]
             for (subentity_idx, subentity_nodes) in pairs(boundaryfunction(BI)(cell))
                 if Base.all(n -> n in facet_nodes, subentity_nodes)
diff --git a/src/exports.jl b/src/exports.jl
index c710861c15..04afb094aa 100644
--- a/src/exports.jl
+++ b/src/exports.jl
@@ -77,6 +77,7 @@ export
     EdgeIndex,
     VertexIndex,
     FacetIndex,
+    InterfaceIndex,
     geometric_interpolation,
     ExclusiveTopology,
     getneighborhood,
diff --git a/src/iterators.jl b/src/iterators.jl
index 0abf2972b1..169ccd48d8 100644
--- a/src/iterators.jl
+++ b/src/iterators.jl
@@ -61,11 +61,16 @@ function CellCache(dh::DofHandler{dim}, flags::UpdateFlags = UpdateFlags()) wher
     return CellCache(flags, get_grid(dh), -1, nodes, coords, dh, celldofs)
 end
 
-function CellCache(sdh::SubDofHandler, flags::UpdateFlags = UpdateFlags())
-    Tv = get_coordinate_type(sdh.dh.grid)
-    return CellCache(flags, sdh.dh.grid, -1, Int[], Tv[], sdh, Int[])
+function CellCache(sdh::SubDofHandler{<:DofHandler{dim}}, flags::UpdateFlags = UpdateFlags()) where {dim}
+    n = ndofs_per_cell(sdh) # dofs and coords will be resized in `reinit!`
+    N = nnodes_per_cell(get_grid(sdh.dh), sdh.cellset[1])
+    nodes = zeros(Int, N)
+    coords = zeros(Vec{dim, get_coordinate_eltype(get_grid(sdh.dh))}, N)
+    celldofs = zeros(Int, n)
+    return CellCache(flags, sdh.dh.grid, -1, nodes, coords, sdh, celldofs)
 end
 
+# TODO: Find a better way to make AllocCheck.jl happy
 function reinit!(cc::CellCache, i::Int)
     cc.cellid = i
     if cc.flags.nodes
@@ -83,6 +88,23 @@ function reinit!(cc::CellCache, i::Int)
     return cc
 end
 
+function reinit!(cc::CellCache{<:Any, <:AbstractGrid, <:SubDofHandler{<:Any, CT}}, i::Int) where {CT}
+    # If we have a DofHandler the cells must be of the same type -> no need to resize
+    cc.cellid = i
+    if cc.flags.nodes
+        cell = getcells(cc.grid, i)::CT
+        _cellnodes!(cc.nodes, cell)
+    end
+    if cc.flags.coords
+        cell = getcells(cc.grid, i)::CT
+        getcoordinates!(cc.coords, cc.grid, cell)
+    end
+    if cc.dh !== nothing && cc.flags.dofs
+        celldofs!(cc.dofs, cc.dh, i)
+    end
+    return cc
+end
+
 # reinit! FEValues with CellCache
 reinit!(cv::CellValues, cc::CellCache) = reinit!(cv, cc.coords)
 reinit!(fv::FacetValues, cc::CellCache, f::Int) = reinit!(fv, cc.coords, f)
@@ -164,22 +186,27 @@ interface. The cache is updated for a new cell by calling `reinit!(cache, facet_
 
 See also [`InterfaceIterator`](@ref).
 """
-struct InterfaceCache{FC <: FacetCache}
-    a::FC
-    b::FC
+struct InterfaceCache{FC_HERE <: FacetCache, FC_THERE <: FacetCache}
+    a::FC_HERE
+    b::FC_THERE
     dofs::Vector{Int}
 end
 
 function InterfaceCache(gridordh::Union{AbstractGrid, AbstractDofHandler})
     fc_a = FacetCache(gridordh)
     fc_b = FacetCache(gridordh)
-    return InterfaceCache(fc_a, fc_b, Int[])
+    return InterfaceCache(fc_a, fc_b, zeros(Int, length(celldofs(fc_a)) + length(celldofs(fc_b))))
+end
+
+function InterfaceCache(sdh_here::SubDofHandler, sdh_there::SubDofHandler)
+    fc_a = FacetCache(sdh_here)
+    fc_b = FacetCache(sdh_there)
+    return InterfaceCache(fc_a, fc_b, zeros(Int, length(celldofs(fc_a)) + length(celldofs(fc_b))))
 end
 
-function reinit!(cache::InterfaceCache, facet_a::BoundaryIndex, facet_b::BoundaryIndex)
-    reinit!(cache.a, facet_a)
-    reinit!(cache.b, facet_b)
-    resize!(cache.dofs, length(celldofs(cache.a)) + length(celldofs(cache.b)))
+function reinit!(cache::InterfaceCache, interface::InterfaceIndex)
+    reinit!(cache.a, FacetIndex(interface.idx[1], interface.idx[2]))
+    reinit!(cache.b, FacetIndex(interface.idx[3], interface.idx[4]))
     for (i, d) in pairs(cache.a.dofs)
         cache.dofs[i] = d
     end
@@ -261,9 +288,12 @@ function CellIterator(gridordh::Union{Grid, DofHandler}, flags::UpdateFlags)
     return CellIterator(gridordh, nothing, flags)
 end
 function CellIterator(sdh::SubDofHandler, flags::UpdateFlags = UpdateFlags())
-    return CellIterator(CellCache(sdh, flags), sdh.cellset)
+    return CellIterator(sdh, sdh.cellset, flags)
 end
 
+function CellIterator(sdh::SubDofHandler, set::IntegerCollection, flags::UpdateFlags = UpdateFlags())
+    return CellIterator(CellCache(sdh, flags), set)
+end
 @inline _getset(ci::CellIterator) = ci.set
 @inline _getcache(ci::CellIterator) = ci.cc
 
@@ -293,9 +323,9 @@ for faceindex in facetset
     # ...
 end
 """
-struct FacetIterator{FC <: FacetCache}
+struct FacetIterator{FC <: FacetCache, SetType <: AbstractVecOrSet{FacetIndex}}
     fc::FC
-    set::OrderedSet{FacetIndex}
+    set::SetType
 end
 
 function FacetIterator(
@@ -309,6 +339,25 @@ function FacetIterator(
     return FacetIterator(FacetCache(gridordh, flags), set)
 end
 
+function FacetIterator(
+        gridordh::Union{<:AbstractGrid, <:DofHandler}, topology::ExclusiveTopology,
+        flags::UpdateFlags = UpdateFlags()
+    )
+    grid = gridordh isa DofHandler ? get_grid(gridordh) : gridordh
+    set = Set(create_boundaryfacetset(grid, topology, _ -> true))
+    return FacetIterator(gridordh, set, flags)
+end
+
+function FacetIterator(
+        sdh::SubDofHandler, topology::ExclusiveTopology,
+        flags::UpdateFlags = UpdateFlags()
+    )
+    grid = get_grid(sdh.dh)
+    set_unfiltered = create_boundaryfacetset(grid, topology, _ -> true)
+    set = Set(filter!(facet -> facet[1] ∈ sdh.cellset, set_unfiltered))
+    return FacetIterator(sdh, set, flags)
+end
+
 @inline _getcache(fi::FacetIterator) = fi.fc
 @inline _getset(fi::FacetIterator) = fi.set
 
@@ -341,40 +390,101 @@ end
     `InterfaceIterator` is stateful and should not be used for things other than `for`-looping
     (e.g. broadcasting over, or collecting the iterator may yield unexpected results).
 """
-struct InterfaceIterator{IC <: InterfaceCache, G <: Grid}
+struct InterfaceIterator{IC <: InterfaceCache, SetType <: AbstractSet{InterfaceIndex}}
     cache::IC
-    grid::G
-    topology::ExclusiveTopology
+    set::SetType
 end
 
+@inline _getcache(ii::InterfaceIterator) = ii.cache
+@inline _getset(ii::InterfaceIterator) = ii.set
+
 function InterfaceIterator(
-        gridordh::Union{Grid, AbstractDofHandler},
+        gridordh::Union{Grid, DofHandler},
+        set::AbstractVecOrSet{InterfaceIndex}
+    )
+    grid = gridordh isa Grid ? gridordh : get_grid(gridordh)
+    if gridordh isa DofHandler
+        # Keep here to maintain same settings as for CellIterator
+        _check_same_celltype(grid, set)
+    end
+    return InterfaceIterator(InterfaceCache(gridordh), set)
+end
+function InterfaceIterator(
+        sdh_here::SubDofHandler,
+        sdh_there::SubDofHandler,
+        set::AbstractVecOrSet{InterfaceIndex}
+    )
+    grid = get_grid(sdh_here.dh)
+    _check_same_celltype(grid, set, getcelltype(grid, first(sdh_here.cellset)), getcelltype(grid, first(sdh_there.cellset)))
+    return InterfaceIterator(InterfaceCache(sdh_here, sdh_there), set)
+end
+function InterfaceIterator(
+        gridordh::Union{Grid, DofHandler},
         topology::ExclusiveTopology = ExclusiveTopology(gridordh isa Grid ? gridordh : get_grid(gridordh))
     )
     grid = gridordh isa Grid ? gridordh : get_grid(gridordh)
-    return InterfaceIterator(InterfaceCache(gridordh), grid, topology)
+    if gridordh isa DofHandler
+        # Keep here to maintain same settings as for CellIterator
+        _check_same_celltype(grid, 1:getncells(grid))
+        @assert length(gridordh.subdofhandlers) == 1 "Use InterfaceIterator(::SubDofHandler, ::SubDofHandler) for subdomain support"
+    end
+    neighborhood = get_facet_facet_neighborhood(topology, grid)
+    fs = facetskeleton(topology, grid)
+    ninterfaces = count(facet -> !isempty(neighborhood[facet[1], facet[2]]), fs)
+    set = Set{InterfaceIndex}()
+    sizehint!(set, ninterfaces)
+    for facet in fs
+        isempty(neighborhood[facet[1], facet[2]]) && continue
+        push!(set, InterfaceIndex(facet[1], facet[2], neighborhood[facet[1], facet[2]][][1], neighborhood[facet[1], facet[2]][][2]))
+    end
+    return InterfaceIterator(gridordh, set)
 end
 
-# Iterator interface
-function Base.iterate(ii::InterfaceIterator{<:Any, <:Grid{sdim}}, state...) where {sdim}
-    neighborhood = get_facet_facet_neighborhood(ii.topology, ii.grid) # TODO: This could be moved to InterfaceIterator constructor (potentially type-instable for non-union or mixed grids)
-    while true
-        it = iterate(facetskeleton(ii.topology, ii.grid), state...)
-        it === nothing && return nothing
-        facet_a, state = it
-        if isempty(neighborhood[facet_a[1], facet_a[2]])
+function InterfaceIterator(
+        sdh_here::SubDofHandler,
+        sdh_there::SubDofHandler,
+        topology::ExclusiveTopology = ExclusiveTopology(get_grid(sdh_here.dh)),
+        # TODO: better name?
+        allow_all::Bool = true # allow interfaces that belong to the other cell according to the cell numbering to be iterated over as if they were "here" not "there"
+    )
+    grid = get_grid(sdh_here.dh)
+    neighborhood = get_facet_facet_neighborhood(topology, grid)
+    fs = facetskeleton(topology, grid)
+    ninterfaces = 0
+    for facet in fs
+        if facet[1] ∈ sdh_here.cellset
+            neighbors = neighborhood[facet[1], facet[2]]
+            isempty(neighbors) && continue
+            neighbors[][1] ∈ sdh_there.cellset || continue
+        elseif allow_all && facet[1] ∈ sdh_there.cellset
+            neighbors = neighborhood[facet[1], facet[2]]
+            isempty(neighbors) && continue
+            neighbors[][1] ∈ sdh_here.cellset || continue
+        else
             continue
         end
-        neighbors = neighborhood[facet_a[1], facet_a[2]]
-        length(neighbors) > 1 && error("multiple neighboring faces not supported yet")
-        facet_b = neighbors[1]
-        reinit!(ii.cache, facet_a, facet_b)
-        return (ii.cache, state)
+        ninterfaces += 1
     end
-    return
+    set = Set{InterfaceIndex}()
+    sizehint!(set, ninterfaces)
+    for facet in fs
+        if facet[1] ∈ sdh_here.cellset
+            neighbors = neighborhood[facet[1], facet[2]]
+            isempty(neighbors) && continue
+            neighbors[][1] ∈ sdh_there.cellset || continue
+            push!(set, InterfaceIndex(facet[1], facet[2], neighborhood[facet[1], facet[2]][][1], neighborhood[facet[1], facet[2]][][2]))
+        elseif allow_all && facet[1] ∈ sdh_there.cellset
+            neighbors = neighborhood[facet[1], facet[2]]
+            isempty(neighbors) && continue
+            neighbors[][1] ∈ sdh_here.cellset || continue
+            push!(set, InterfaceIndex(neighborhood[facet[1], facet[2]][][1], neighborhood[facet[1], facet[2]][][2], facet[1], facet[2]))
+        else
+            continue
+        end
+    end
+    return InterfaceIterator(sdh_here, sdh_there, set)
 end
 
-
 # Iterator interface for CellIterator/FacetIterator
 const GridIterators{C} = Union{CellIterator{C}, FacetIterator{C}, InterfaceIterator{C}}
 
@@ -408,3 +518,15 @@ function _check_same_celltype(grid::AbstractGrid, facetset::AbstractVecOrSet{<:B
     end
     return
 end
+
+function _check_same_celltype(
+        grid::AbstractGrid, interfaceset::AbstractVecOrSet{InterfaceIndex},
+        celltype_here::Type{<:AbstractCell} = getcelltype(grid, first(interfaceset)[1]),
+        celltype_there::Type{<:AbstractCell} = getcelltype(grid, first(interfaceset)[3])
+    )
+    isconcretetype(getcelltype(grid)) && return nothing # Short circuit check
+    if !all(getcelltype(grid, interface[1]) == celltype_here && getcelltype(grid, interface[3]) == celltype_there for interface in interfaceset)
+        error("The cells in the set (set of InterfaceIndex) are not all of the same celltype on each side.")
+    end
+    return
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 57cd82d8a7..0f339c752c 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -11,6 +11,7 @@ using StaticArrays
 using OrderedCollections
 using WriteVTK
 import Metis
+using AllocCheck
 
 include("test_utils.jl")
 
@@ -20,6 +21,7 @@ include("test_interpolations.jl")
 include("test_cellvalues.jl")
 include("test_facevalues.jl")
 include("test_interfacevalues.jl")
+include("test_iterators.jl")
 include("test_quadrules.jl")
 include("test_assemble.jl")
 include("test_dofs.jl")
diff --git a/test/test_grid_addboundaryset.jl b/test/test_grid_addboundaryset.jl
index ed3ca549ad..ed85602f50 100644
--- a/test/test_grid_addboundaryset.jl
+++ b/test/test_grid_addboundaryset.jl
@@ -200,4 +200,21 @@
         addboundaryfacetset!(grid, topology, "test_boundary_facetset", filter_function)
         @test getfacetset(grid, "test_boundary_facetset") == Ferrite.create_boundaryfacetset(grid, topology, filter_function)
     end
+
+    @testset "addboundaryset Mixed grid" begin
+        nodes = [Node((-1.0, 0.0)), Node((0.0, 0.0)), Node((1.0, 0.0)), Node((-1.0, 1.0)), Node((0.0, 1.0))]
+        cells = [
+            Quadrilateral((1, 2, 5, 4)),
+            Triangle((3, 5, 2)),
+        ]
+        grid = Grid(cells, nodes)
+        topology = ExclusiveTopology(grid)
+        @test extractboundary(grid, topology) == extractboundarycheck(grid)
+
+        filter_function(x) = x[1] > 0
+        addboundaryvertexset!(grid, topology, "test_boundary_vertexset", filter_function)
+        @test getvertexset(grid, "test_boundary_vertexset") == Ferrite.create_boundaryvertexset(grid, topology, filter_function)
+        addboundaryfacetset!(grid, topology, "test_boundary_facetset", filter_function)
+        @test getfacetset(grid, "test_boundary_facetset") == Ferrite.create_boundaryfacetset(grid, topology, filter_function)
+    end
 end
diff --git a/test/test_grid_dofhandler_vtk.jl b/test/test_grid_dofhandler_vtk.jl
index 9f1ad073cd..839114cf50 100644
--- a/test/test_grid_dofhandler_vtk.jl
+++ b/test/test_grid_dofhandler_vtk.jl
@@ -232,12 +232,12 @@ end
     # InterfaceCache
     grid = generate_grid(Quadrilateral, (2, 1))
     ic = InterfaceCache(grid)
-    reinit!(ic, FaceIndex(1, 2), FaceIndex(2, 4))
+    reinit!(ic, InterfaceIndex(1, 2, 2, 4))
     @test interfacedofs(ic) == Int[] # Empty because no DofHandler given
     ip = DiscontinuousLagrange{RefQuadrilateral, 1}()
     dh = DofHandler(grid); add!(dh, :u, ip); close!(dh)
     ic = InterfaceCache(dh)
-    reinit!(ic, FaceIndex(1, 2), FaceIndex(2, 4))
+    reinit!(ic, InterfaceIndex(1, 2, 2, 4))
     @test interfacedofs(ic) == collect(1:8)
     # Mixed Elements
     dim = 2
@@ -253,8 +253,8 @@ end
     sdh1 = SubDofHandler(dh, Set([1])); add!(sdh1, :u, ip1)
     sdh2 = SubDofHandler(dh, Set([2])); add!(sdh2, :u, ip2)
     close!(dh)
-    ic = InterfaceCache(dh)
-    reinit!(ic, FaceIndex(1, 2), FaceIndex(2, 3))
+    ic = InterfaceCache(sdh1, sdh2)
+    reinit!(ic, InterfaceIndex(1, 2, 2, 3))
     @test interfacedofs(ic) == collect(1:7)
     # Unit test of some utilities
     mixed_grid = Grid(
diff --git a/test/test_interfacevalues.jl b/test/test_interfacevalues.jl
index 8fb31dc741..6e0aced047 100644
--- a/test/test_interfacevalues.jl
+++ b/test/test_interfacevalues.jl
@@ -1,13 +1,12 @@
 @testset "InterfaceValues" begin
-    function test_interfacevalues(grid::Ferrite.AbstractGrid, iv::InterfaceValues; tol = 0)
+    function test_interfacevalues(iv::InterfaceValues, args...; tol = 0)
         ip_here = Ferrite.function_interpolation(iv.here)
         ip_there = Ferrite.function_interpolation(iv.there)
         rdim = Ferrite.getrefdim(ip_here)
         n_basefuncs = getnbasefunctions(ip_here) + getnbasefunctions(ip_there)
 
         @test getnbasefunctions(iv) == n_basefuncs
-
-        for ic in InterfaceIterator(grid)
+        for ic in InterfaceIterator(args...)
             reinit!(iv, ic)
             coords_here, coords_there = getcoordinates(ic)
             nqp = getnquadpoints(iv)
@@ -152,7 +151,7 @@
         for func_interpol in (scalar_interpol, VectorizedInterpolation(scalar_interpol))
             iv = cell_shape ∈ (QuadraticLine, QuadraticQuadrilateral, QuadraticTriangle, QuadraticTetrahedron) ?
                 InterfaceValues(quad_rule, func_interpol, ip) : InterfaceValues(quad_rule, func_interpol)
-            test_interfacevalues(grid, iv)
+            test_interfacevalues(iv, grid)
         end
     end
     # Custom quadrature
@@ -182,23 +181,45 @@
         end
         for func_interpol in (scalar_interpol, VectorizedInterpolation(scalar_interpol))
             iv = InterfaceValues(quad_rule, func_interpol)
-            test_interfacevalues(grid, iv; tol = 5 * eps(Float64))
+            test_interfacevalues(iv, grid; tol = 5 * eps(Float64))
         end
     end
-    # @testset "Mixed elements 2D grids" begin # TODO: this shouldn't work because it should change the FacetValues object
-    #     dim = 2
-    #     nodes = [Node((-1.0, 0.0)), Node((0.0, 0.0)), Node((1.0, 0.0)), Node((-1.0, 1.0)), Node((0.0, 1.0))]
-    #     cells = [
-    #                 Quadrilateral((1,2,5,4)),
-    #                 Triangle((3,5,2)),
-    #             ]
+    @testset "Mixed elements 2D grids" begin
+        dim = 2
+        nodes = [Node((-1.0, 0.0)), Node((0.0, 0.0)), Node((1.0, 0.0)), Node((-1.0, 1.0)), Node((0.0, 1.0))]
+        cells = [
+            Quadrilateral((1, 2, 5, 4)),
+            Triangle((3, 5, 2)),
+        ]
 
-    #     grid = Grid(cells, nodes)
-    #     topology = ExclusiveTopology(grid)
-    #     test_interfacevalues(grid,
-    #     DiscontinuousLagrange{RefQuadrilateral, 1}(), FacetQuadratureRule{RefQuadrilateral}(2),
-    #     DiscontinuousLagrange{RefTriangle, 1}(), FacetQuadratureRule{RefTriangle}(2))
-    # end
+        grid = Grid(cells, nodes)
+        topology = ExclusiveTopology(grid)
+        qr_facet_quad = FacetQuadratureRule{RefQuadrilateral}(2)
+        qr_facet_tri = FacetQuadratureRule{RefTriangle}(2)
+        ip_quad = DiscontinuousLagrange{RefQuadrilateral, 1}()
+        ip_tri = DiscontinuousLagrange{RefTriangle, 1}()
+        dh = DofHandler(grid)
+        sdh_quad = SubDofHandler(dh, Set([1]))
+        sdh_tri = SubDofHandler(dh, Set([2]))
+        add!(sdh_quad, :u, ip_quad)
+        add!(sdh_tri, :u, ip_tri)
+        qr_collection = (qr_facet_quad, qr_facet_tri)
+        ip_collection = (ip_quad, ip_tri)
+        close!(dh)
+        for (sdh_here_idx, sdh_here) in enumerate(dh.subdofhandlers)
+            for (sdh_there_idx, sdh_there) in enumerate(dh.subdofhandlers)
+                sdh_there_idx < sdh_here_idx && continue
+                iv = InterfaceValues(
+                    qr_collection[sdh_here_idx], ip_collection[sdh_here_idx],
+                    qr_collection[sdh_there_idx], ip_collection[sdh_there_idx]
+                )
+                test_interfacevalues(
+                    iv, sdh_here, sdh_there
+                )
+            end
+        end
+
+    end
     @testset "Unordered nodes 3D" begin
         dim = 2
         nodes = [
@@ -214,8 +235,8 @@
 
         grid = Grid(cells, nodes)
         test_interfacevalues(
-            grid,
-            InterfaceValues(FacetQuadratureRule{RefHexahedron}(2), DiscontinuousLagrange{RefHexahedron, 1}())
+            InterfaceValues(FacetQuadratureRule{RefHexahedron}(2), DiscontinuousLagrange{RefHexahedron, 1}()),
+            grid
         )
     end
     @testset "Interface dof_range" begin
diff --git a/test/test_iterators.jl b/test/test_iterators.jl
new file mode 100644
index 0000000000..cd011799b7
--- /dev/null
+++ b/test/test_iterators.jl
@@ -0,0 +1,302 @@
+@testset "Iterators" begin
+    ###############################################################################
+    # Function to test for possible allocations/dynamic dispatch when iterating   #
+    ###############################################################################
+    @check_allocs function testallocs_iterate(iterator)
+        for _ in iterator end
+        return nothing
+    end
+    @testset "CellIterator" begin
+        @testset "Single domain" begin
+            grid = generate_grid(Hexahedron, (2, 2, 2))
+            ip = Lagrange{RefHexahedron, 1}()
+            dh = DofHandler(grid)
+            add!(dh, :u, ip)
+            close!(dh)
+
+            ci_dh = CellIterator(dh)
+            ci_grid = CellIterator(grid)
+            @test ci_dh.set == ci_grid.set
+            # Test for allocations
+            # TODO: find a way to avoid resize for single sdh
+            # @test (testallocs_iterate(ci_dh)) === nothing
+            # Iterators over grid can allocate due to potential resize!
+            @test_throws AllocCheckFailure testallocs_iterate(ci_grid)
+        end
+
+        @testset "subdomains" begin
+            @testset "same cell type" begin
+                grid = generate_grid(Hexahedron, (2, 2, 2))
+                ip = Lagrange{RefHexahedron, 1}()
+                dh = DofHandler(grid)
+                sdh1 = SubDofHandler(dh, Set(collect(1:4)))
+                sdh2 = SubDofHandler(dh, Set(collect(5:8)))
+                add!(sdh1, :u, ip^3)
+                add!(sdh2, :p, ip)
+                close!(dh)
+
+                ci_sdh1 = CellIterator(sdh1, Set(collect(1:4)))
+                ci_sdh2 = CellIterator(sdh2, Set(collect(5:8)))
+                ci_dh = CellIterator(dh)
+
+                # Test for allocations
+                @test (testallocs_iterate(ci_sdh1)) === nothing
+                @test (testallocs_iterate(ci_sdh2)) === nothing
+                @test_throws AllocCheckFailure testallocs_iterate(ci_dh)
+            end
+
+            @testset "mixed cell types" begin
+                nodes = [Node((-1.0, 0.0)), Node((0.0, 0.0)), Node((1.0, 0.0)), Node((-1.0, 1.0)), Node((0.0, 1.0))]
+                cells = [
+                    Quadrilateral((1, 2, 5, 4)),
+                    Triangle((3, 5, 2)),
+                ]
+                grid = Grid(cells, nodes)
+                ip_quad = Lagrange{RefQuadrilateral, 1}()
+                ip_tri = Lagrange{RefTriangle, 1}()
+                dh = DofHandler(grid)
+                sdh_quad = SubDofHandler(dh, Set([1]))
+                sdh_tri = SubDofHandler(dh, Set([2]))
+                add!(sdh_quad, :u, ip_quad)
+                add!(sdh_tri, :u, ip_tri)
+                close!(dh)
+
+                ci_quad = CellIterator(sdh_quad, Set([1]))
+                ci_tri = CellIterator(sdh_tri, Set([2]))
+                ci_dh = CellIterator(dh, Set([1]))
+
+                # Test for allocations
+                try
+                    testallocs_iterate(ci_quad)
+
+                catch err
+                    @info err.errors[3]
+                end
+                @test (testallocs_iterate(ci_quad)) === nothing
+                @test (testallocs_iterate(ci_tri)) === nothing
+                @test_throws AllocCheckFailure testallocs_iterate(ci_dh)
+            end
+        end
+
+        @testset "error paths" begin
+            nodes = [Node((-1.0, 0.0)), Node((0.0, 0.0)), Node((1.0, 0.0)), Node((-1.0, 1.0)), Node((0.0, 1.0))]
+            cells = [
+                Quadrilateral((1, 2, 5, 4)),
+                Triangle((3, 5, 2)),
+            ]
+            grid = Grid(cells, nodes)
+            ip_quad = Lagrange{RefQuadrilateral, 1}()
+            ip_tri = Lagrange{RefTriangle, 1}()
+            dh = DofHandler(grid)
+            sdh_quad = SubDofHandler(dh, Set([1]))
+            sdh_tri = SubDofHandler(dh, Set([2]))
+            add!(sdh_quad, :u, ip_quad)
+            add!(sdh_tri, :u, ip_tri)
+            close!(dh)
+            @test_throws ErrorException("The cells in the cellset are not all of the same celltype.") CellIterator(dh)
+        end
+    end
+
+    @testset "FacetIterator" begin
+        @testset "Single domain" begin
+            grid = generate_grid(Hexahedron, (2, 2, 2))
+            ip = Lagrange{RefHexahedron, 1}()
+            dh = DofHandler(grid)
+            add!(dh, :u, ip)
+            close!(dh)
+            topology = ExclusiveTopology(grid)
+            fi_dh = FacetIterator(dh, topology)
+            fi_grid = FacetIterator(grid, topology)
+            @test fi_dh.set == fi_grid.set
+
+            # Test for allocations
+            # TODO: find a way to avoid resize for single sdh
+            # @test (testallocs_iterate(fi_dh)) === nothing
+            # Iterators over grid can allocate due to potential resize!
+            @test_throws AllocCheckFailure testallocs_iterate(fi_grid)
+        end
+
+        @testset "subdomains" begin
+            @testset "same cell type" begin
+                grid = generate_grid(Hexahedron, (2, 2, 2))
+                ip = Lagrange{RefHexahedron, 1}()
+                dh = DofHandler(grid)
+                sdh1 = SubDofHandler(dh, Set(collect(1:4)))
+                sdh2 = SubDofHandler(dh, Set(collect(5:8)))
+                add!(sdh1, :u, ip^3)
+                add!(sdh2, :p, ip)
+                close!(dh)
+                topology = ExclusiveTopology(grid)
+
+                fi_sdh1 = FacetIterator(sdh1, topology)
+                fi_sdh2 = FacetIterator(sdh2, topology)
+                @test isempty(fi_sdh1.set ∩ fi_sdh2.set)
+                @test length(fi_sdh1.set) == 12
+                # Test for allocations
+                @test (testallocs_iterate(fi_sdh1)) === nothing
+                @test (testallocs_iterate(fi_sdh2)) === nothing
+            end
+
+            @testset "mixed cell types" begin
+                nodes = [Node((-1.0, 0.0)), Node((0.0, 0.0)), Node((1.0, 0.0)), Node((-1.0, 1.0)), Node((0.0, 1.0))]
+                cells = [
+                    Quadrilateral((1, 2, 5, 4)),
+                    Triangle((3, 5, 2)),
+                ]
+                grid = Grid(cells, nodes)
+                topology = ExclusiveTopology(grid)
+                ip_quad = Lagrange{RefQuadrilateral, 1}()
+                ip_tri = Lagrange{RefTriangle, 1}()
+                dh = DofHandler(grid)
+                sdh_quad = SubDofHandler(dh, Set([1]))
+                sdh_tri = SubDofHandler(dh, Set([2]))
+                add!(sdh_quad, :u, ip_quad)
+                add!(sdh_tri, :u, ip_tri)
+                close!(dh)
+
+                fi_quad = FacetIterator(sdh_quad, topology)
+                fi_tri = FacetIterator(sdh_tri, topology)
+                @test isempty(fi_quad.set ∩ fi_tri.set)
+                # Test for allocations
+                @test (testallocs_iterate(fi_quad)) === nothing
+                @test (testallocs_iterate(fi_tri)) === nothing
+            end
+        end
+
+        @testset "FacetIndex" begin
+            # checkmate, codecov bot!
+            @test isequal(FacetIndex(1, 2), (1, 2))
+            @test isequal((1, 2), FacetIndex(1, 2))
+        end
+
+        @testset "error paths" begin
+            nodes = [Node((-1.0, 0.0)), Node((0.0, 0.0)), Node((1.0, 0.0)), Node((-1.0, 1.0)), Node((0.0, 1.0))]
+            cells = [
+                Quadrilateral((1, 2, 5, 4)),
+                Triangle((3, 5, 2)),
+            ]
+            grid = Grid(cells, nodes)
+            topology = ExclusiveTopology(grid)
+            ip_quad = Lagrange{RefQuadrilateral, 1}()
+            ip_tri = Lagrange{RefTriangle, 1}()
+            dh = DofHandler(grid)
+            sdh_quad = SubDofHandler(dh, Set([1]))
+            sdh_tri = SubDofHandler(dh, Set([2]))
+            add!(sdh_quad, :u, ip_quad)
+            add!(sdh_tri, :u, ip_tri)
+            close!(dh)
+            @test_throws ErrorException("The cells in the set (set of FacetIndex) are not all of the same celltype.") FacetIterator(dh, topology)
+        end
+    end
+
+    @testset "InterfaceIterator" begin
+        @testset "Single domain" begin
+            grid = generate_grid(Hexahedron, (1000, 2, 1))
+            ip = Lagrange{RefHexahedron, 1}()
+            dh = DofHandler(grid)
+            add!(dh, :u, ip)
+            close!(dh)
+            topology = ExclusiveTopology(grid)
+
+            ii_dh_top = InterfaceIterator(dh, topology)
+            ii_dh = InterfaceIterator(dh)
+            ii_grid_top = InterfaceIterator(grid, topology)
+            ii_grid = InterfaceIterator(grid)
+            @test ii_dh.set == ii_dh_top.set == ii_grid.set == ii_grid_top.set
+            # Test that topology has no effect on iterator type
+            @test typeof(ii_dh) == typeof(ii_dh_top)
+            @test typeof(ii_grid) == typeof(ii_grid_top)
+            # Test for allocations
+            # TODO: find a way to avoid resize for single sdh
+            # @test (testallocs_iterate(ii_dh_top)) === nothing
+            # @test (testallocs_iterate(ii_dh)) === nothing
+            # Iterators over grid can allocate due to potential resize!
+            @test_throws AllocCheckFailure testallocs_iterate(ii_grid)
+            @test_throws AllocCheckFailure testallocs_iterate(ii_grid_top)
+        end
+
+        @testset "subdomains" begin
+            @testset "same cell type" begin
+                grid = generate_grid(Hexahedron, (1000, 2, 1))
+                ip = Lagrange{RefHexahedron, 1}()
+                dh = DofHandler(grid)
+                sdh1 = SubDofHandler(dh, Set([collect(1:500)..., collect(1501:2000)...]))
+                sdh2 = SubDofHandler(dh, Set([collect(501:1500)...]))
+                add!(sdh1, :u, ip)
+                add!(sdh2, :u, ip)
+                close!(dh)
+                topology = ExclusiveTopology(grid)
+
+                ii_sdh_top = InterfaceIterator(sdh1, sdh2, topology)
+                ii_sdh = InterfaceIterator(sdh1, sdh2)
+                ii_same_sdh = InterfaceIterator(sdh1, sdh1)
+                @test ii_sdh.set == ii_sdh_top.set
+                @test length(ii_sdh.set) == 1002
+                @test count(interface_index -> interface_index[1] == interface_index[3] - 1000, ii_sdh.set) == 500
+                @test count(interface_index -> interface_index[1] == interface_index[3] + 1000, ii_sdh.set) == 500
+                ii_sdh_flipped = InterfaceIterator(sdh2, sdh1)
+                @test all(interface -> InterfaceIndex(interface[3], interface[4], interface[1], interface[2]) ∈ ii_sdh_flipped.set, ii_sdh.set)
+                # Test that topology has no effect on iterator type
+                @test typeof(ii_sdh) == typeof(ii_sdh_top)
+                # Test for allocations
+                @test (testallocs_iterate(ii_sdh)) === nothing
+                @test (testallocs_iterate(ii_same_sdh)) === nothing
+            end
+
+            @testset "mixed cell types" begin
+                nodes = [Node((-1.0, 0.0)), Node((0.0, 0.0)), Node((1.0, 0.0)), Node((-1.0, 1.0)), Node((0.0, 1.0))]
+                cells = [
+                    Quadrilateral((1, 2, 5, 4)),
+                    Triangle((3, 5, 2)),
+                ]
+                grid = Grid(cells, nodes)
+                topology = ExclusiveTopology(grid)
+                ip_quad = Lagrange{RefQuadrilateral, 1}()
+                ip_tri = Lagrange{RefTriangle, 1}()
+                dh = DofHandler(grid)
+                sdh_quad = SubDofHandler(dh, Set([1]))
+                sdh_tri = SubDofHandler(dh, Set([2]))
+                add!(sdh_quad, :u, ip_quad)
+                add!(sdh_tri, :u, ip_tri)
+                close!(dh)
+
+                ii_sdh_top = InterfaceIterator(sdh_quad, sdh_tri, topology)
+                ii_sdh = InterfaceIterator(sdh_quad, sdh_tri)
+                @test ii_sdh.set == ii_sdh_top.set
+                @test ii_sdh.set == Set([InterfaceIndex(1, 2, 2, 2)])
+                ii_sdh_flipped = InterfaceIterator(sdh_tri, sdh_quad)
+                @test ii_sdh_flipped.set == Set([InterfaceIndex(2, 2, 1, 2)])
+                # Test that topology has no effect on iterator type
+                @test typeof(ii_sdh) == typeof(ii_sdh_top)
+                # Test for allocations
+                @test (testallocs_iterate(ii_sdh_top)) === nothing
+                @test (testallocs_iterate(ii_sdh)) === nothing
+            end
+        end
+
+        @testset "InterfaceIndex" begin
+            # checkmate, codecov bot!
+            @test isequal(InterfaceIndex(1, 2, 3, 4), (1, 2, 3, 4))
+            @test isequal((1, 2, 3, 4), InterfaceIndex(1, 2, 3, 4))
+        end
+
+        @testset "error paths" begin
+            nodes = [Node((-1.0, 0.0)), Node((0.0, 0.0)), Node((1.0, 0.0)), Node((-1.0, 1.0)), Node((0.0, 1.0))]
+            cells = [
+                Quadrilateral((1, 2, 5, 4)),
+                Triangle((3, 5, 2)),
+            ]
+            grid = Grid(cells, nodes)
+            topology = ExclusiveTopology(grid)
+            ip_quad = Lagrange{RefQuadrilateral, 1}()
+            ip_tri = Lagrange{RefTriangle, 1}()
+            dh = DofHandler(grid)
+            sdh_quad = SubDofHandler(dh, Set([1]))
+            sdh_tri = SubDofHandler(dh, Set([2]))
+            add!(sdh_quad, :u, ip_quad)
+            add!(sdh_tri, :u, ip_tri)
+            close!(dh)
+            @test_throws ErrorException("The cells in the set (set of InterfaceIndex) are not all of the same celltype on each side.") InterfaceIterator(sdh_quad, sdh_tri, Set([InterfaceIndex(2, 2, 1, 2)]))
+        end
+    end
+end
diff --git a/test/test_mixeddofhandler.jl b/test/test_mixeddofhandler.jl
index ec4a7206b8..18eb2127aa 100644
--- a/test/test_mixeddofhandler.jl
+++ b/test/test_mixeddofhandler.jl
@@ -714,13 +714,19 @@ function test_celliterator_on_true_subdomain_smoketest()
 
     for cell in CellIterator(dh, [1, 2, 3])
     end
-
-    for cell in CellIterator(dh)
-        if cellid(cell) <= 3
-            @test length(celldofs(cell)) == getnbasefunctions(ip)
-        else
-            @test length(celldofs(cell)) == 0
-        end
+    # I smell smoke here
+    # for cell in CellIterator(dh)
+    #     if cellid(cell) <= 3
+    #         @test length(celldofs(cell)) == getnbasefunctions(ip)
+    #     else
+    #         @test length(celldofs(cell)) == 0
+    #     end
+    # end
+    for cell in CellIterator(sdh)
+        @test length(celldofs(cell)) == getnbasefunctions(ip)
+    end
+    for cell in CellIterator(grid)
+        @test length(celldofs(cell)) == 0
     end
     return
 end