diff --git a/CHANGELOG.md b/CHANGELOG.md
index 11add1e764..3ecacc0ac0 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -446,10 +446,18 @@ more discussion).
   `facedof_indices`, `volumedof_indices` (and similar) according to these definitions.
 
 - `Ferrite.getdim` has been changed into `Ferrite.getrefdim` for getting the dimension of the reference shape
-  and `Ferrite.getspatialdim` to get the spatial dimension (of the grid). [#943][github-943]
+  and `Ferrite.getspatialdim` to get the spatial dimension (of the grid). ([#943][github-943])
 
 - `Ferrite.getfielddim(::AbstractDofHandler, args...)` has been renamed to `Ferrite.n_components`.
-  [#943][github-943]
+  ([#943][github-943])
+
+- The constructor for `ExclusiveTopology` only accept an `AbstractGrid` as input,
+  removing the alternative of providing a `Vector{<:AbstractCell}`, as knowing the
+  spatial dimension is required for correct code paths.
+  Furthermore, it uses a new internal data structure, `ArrayOfVectorViews`, to store the neighborhood
+  information more efficiently The datatype for the neighborhood has thus changed to a view of a vector,
+  instead of the now removed `EntityNeighborhood` container. This also applies to `vertex_star_stencils`.
+  ([#974][github-974]).
 
 - `project(::L2Projector, data, qr_rhs)` now expects data to be indexed by the cellid, as opposed to
   the index in the vector of cellids passed to the `L2Projector`. The data may be passed as an
@@ -1037,3 +1045,4 @@ poking into Ferrite internals:
 [github-943]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/943
 [github-949]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/949
 [github-953]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/953
+[github-974]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/974
diff --git a/docs/src/devdocs/index.md b/docs/src/devdocs/index.md
index 9c16b83d7c..bfd680dcc8 100644
--- a/docs/src/devdocs/index.md
+++ b/docs/src/devdocs/index.md
@@ -5,5 +5,5 @@ developing the library.
 
 ```@contents
 Depth = 1
-Pages = ["reference_cells.md", "interpolations.md", "elements.md", "FEValues.md", "dofhandler.md", "performance.md"]
+Pages = ["reference_cells.md", "interpolations.md", "elements.md", "FEValues.md", "dofhandler.md", "performance.md", "special_datastructures.md"]
 ```
diff --git a/docs/src/devdocs/special_datastructures.md b/docs/src/devdocs/special_datastructures.md
new file mode 100644
index 0000000000..e642e5a1be
--- /dev/null
+++ b/docs/src/devdocs/special_datastructures.md
@@ -0,0 +1,16 @@
+# Special data structures
+
+## `ArrayOfVectorViews`
+`ArrayOfVectorViews` is a data structure representing an `Array` of
+vector views (specifically `SubArray{T, 1} where T`). By arranging all
+data (of type `T`) continuously in memory, this will significantly reduce
+the garbage collection time compared to using an `Array{AbstractVector{T}}`. While the data in each view can be mutated, the length of each view is
+fixed after construction.
+This data structure provides two features not provided by `ArraysOfArrays.jl`: Support of matrices and higher order arrays for storing vectors
+of different dimensions and efficient construction when the number of elements in each view is not known in advance.
+
+```@docs
+Ferrite.ArrayOfVectorViews
+Ferrite.ConstructionBuffer
+Ferrite.push_at_index!
+```
diff --git a/docs/src/literate-gallery/topology_optimization.jl b/docs/src/literate-gallery/topology_optimization.jl
index 823756d13e..a03917019f 100644
--- a/docs/src/literate-gallery/topology_optimization.jl
+++ b/docs/src/literate-gallery/topology_optimization.jl
@@ -200,11 +200,11 @@ function cache_neighborhood(dh, topology)
         nbg = zeros(Int,_nfacets)
         i = cellid(element)
         for j in 1:_nfacets
-            nbg_cellid = getcells(getneighborhood(topology, dh.grid, FacetIndex(i,j)))
+            nbg_cellid = getneighborhood(topology, dh.grid, FacetIndex(i,j))
             if(!isempty(nbg_cellid))
-                nbg[j] = first(nbg_cellid) # assuming only one face neighbor per cell
+                nbg[j] = first(nbg_cellid)[1] # assuming only one face neighbor per cell
             else # boundary face
-                nbg[j] = first(getcells(getneighborhood(topology, dh.grid, FacetIndex(i,opp[j]))))
+                nbg[j] = first(getneighborhood(topology, dh.grid, FacetIndex(i,opp[j])))[1]
             end
         end
 
diff --git a/src/CollectionsOfViews.jl b/src/CollectionsOfViews.jl
new file mode 100644
index 0000000000..a85691af06
--- /dev/null
+++ b/src/CollectionsOfViews.jl
@@ -0,0 +1,156 @@
+module CollectionsOfViews
+
+export ArrayOfVectorViews, push_at_index!, ConstructionBuffer
+
+# `AdaptiveRange` and `ConstructionBuffer` are used to efficiently build up an `ArrayOfVectorViews`
+# when the size of each view is unknown.
+struct AdaptiveRange
+    start::Int
+    ncurrent::Int
+    nmax::Int
+end
+
+struct ConstructionBuffer{T, N}
+    indices::Array{AdaptiveRange, N}
+    data::Vector{T}
+    sizehint::Int
+end
+
+"""
+    ConstructionBuffer(data::Vector, dims::NTuple{N, Int}, sizehint)
+
+Create a buffer for creating an [`ArrayOfVectorViews`](@ref), representing an array with `N` axes.
+`sizehint` sets the number of elements in `data` allocated when a new index is added via `push_at_index!`,
+or when the current storage for the index is full, how much many additional elements are reserved for that index.
+Any content in `data` is overwritten, but performance is improved by pre-allocating it to a reasonable size or
+by `sizehint!`ing it.
+"""
+function ConstructionBuffer(data::Vector, dims::NTuple{<:Any, Int}, sizehint::Int)
+    indices = fill(AdaptiveRange(0, 0, 0), dims)
+    return ConstructionBuffer(indices, empty!(data), sizehint)
+end
+
+"""
+    push_at_index!(b::ConstructionBuffer, val, indices::Int...)
+
+`push!` the value `val` to the `Vector` view at the index given by `indices`, typically called
+inside the [`ArrayOfVectorViews`](@ref) constructor do-block. But can also be used when manually
+creating a `ConstructionBuffer`.
+"""
+function push_at_index!(b::ConstructionBuffer, val, indices::Vararg{Int, N}) where {N}
+    r = getindex(b.indices, indices...)
+    n = length(b.data)
+    if r.start == 0
+        # `indices...` not previously added, allocate new space for it at the end of `b.data`
+        resize!(b.data, n + b.sizehint)
+        b.data[n+1] = val
+        setindex!(b.indices, AdaptiveRange(n + 1, 1, b.sizehint), indices...)
+    elseif r.ncurrent == r.nmax
+        # We have used up our space, move data associated with `indices...` to the end of `b.data`
+        resize!(b.data, n + r.nmax + b.sizehint)
+        for i in 1:r.ncurrent
+            b.data[n + i] = b.data[r.start + i - 1]
+        end
+        b.data[n + r.ncurrent + 1] = val
+        setindex!(b.indices, AdaptiveRange(n + 1, r.ncurrent + 1, r.nmax + b.sizehint), indices...)
+    else # We have space in an already allocated section
+        b.data[r.start + r.ncurrent] = val
+        setindex!(b.indices, AdaptiveRange(r.start, r.ncurrent + 1, r.nmax), indices...)
+    end
+    return b
+end
+
+struct ArrayOfVectorViews{T, N} <: AbstractArray{SubArray{T, 1, Vector{T}, Tuple{UnitRange{Int64}}, true}, N}
+    indices::Vector{Int}
+    data::Vector{T}
+    lin_idx::LinearIndices{N, NTuple{N, Base.OneTo{Int}}}
+    function ArrayOfVectorViews{T, N}(indices::Vector{Int}, data::Vector{T}, lin_idx::LinearIndices{N}) where {T, N}
+        return new{T, N}(indices, data, lin_idx)
+    end
+end
+
+# AbstractArray interface (https://docs.julialang.org/en/v1/manual/interfaces/#man-interface-array)
+Base.size(cv::ArrayOfVectorViews) = size(cv.lin_idx)
+@inline function Base.getindex(cv::ArrayOfVectorViews, linear_index::Int)
+    @boundscheck checkbounds(cv.lin_idx, linear_index)
+    return @inbounds view(cv.data, cv.indices[linear_index]:(cv.indices[linear_index+1]-1))
+end
+@inline function Base.getindex(cv::ArrayOfVectorViews, idx...)
+    linear_index = getindex(cv.lin_idx, idx...)
+    return @inbounds getindex(cv, linear_index)
+end
+Base.IndexStyle(::Type{<:ArrayOfVectorViews{<:Any, N}}) where N = Base.IndexStyle(Array{Int, N})
+
+# Constructors
+"""
+    ArrayOfVectorViews(f!::Function, data::Vector{T}, dims::NTuple{N, Int}; sizehint)
+
+Create an `ArrayOfVectorViews` to store many vector views of potentially different sizes,
+emulating an `Array{Vector{T}, N}` with size `dims`. However, it avoids allocating each vector individually
+by storing all data in `data`, and instead of `Vector{T}`, the each element is a `typeof(view(data, 1:2))`.
+
+When the length of each vector is unknown, the `ArrayOfVectorViews` can be created reasonably efficient
+with the following do-block, which creates an intermediate `buffer::ConstructionBuffer` supporting the
+[`push_at_index!`](@ref) function.
+```
+vector_views = ArrayOfVectorViews(data, dims; sizehint) do buffer
+    for (ind, val) in some_data
+        push_at_index!(buffer, val, ind)
+    end
+end
+```
+`sizehint` tells how much space to allocate for the index `ind` if no `val` has been added to that index before,
+or how much more space to allocate in case all previously allocated space for `ind` has been used up.
+"""
+function ArrayOfVectorViews(f!::F, data::Vector, dims::Tuple; sizehint = nothing) where {F <: Function}
+    sizehint === nothing && error("Providing sizehint is mandatory")
+    b = ConstructionBuffer(data, dims, sizehint)
+    f!(b)
+    return ArrayOfVectorViews(b)
+end
+
+"""
+    ArrayOfVectorViews(b::CollectionsOfViews.ConstructionBuffer)
+
+Creates the `ArrayOfVectorViews` directly from the `ConstructionBuffer` that was manually created and filled.
+"""
+function ArrayOfVectorViews(b::ConstructionBuffer{T}) where T
+    indices = Vector{Int}(undef, length(b.indices) + 1)
+    lin_idx = LinearIndices(b.indices)
+    data_length = sum(ar.ncurrent for ar in b.indices; init=0)
+    data = Vector{T}(undef, data_length)
+    data_index = 1
+    for (idx, ar) in pairs(b.indices)
+        copyto!(data, data_index, b.data, ar.start, ar.ncurrent)
+        indices[lin_idx[idx]] = data_index
+        data_index += ar.ncurrent
+    end
+    indices[length(indices)] = data_index
+    # Since user-code in the constructor function has access to `b`, setting dimensions to
+    # zero here allows GC:ing the data in `b` even in cases when the compiler cannot
+    # guarantee that it is unreachable.
+    resize!(b.data, 0); sizehint!(b.data, 0)
+    isa(b.indices, Vector) && (resize!(b.indices, 0); sizehint!(b.indices, 0))
+    return ArrayOfVectorViews(indices, data, lin_idx)
+end
+
+"""
+    ArrayOfVectorViews(indices::Vector{Int}, data::Vector{T}, lin_idx::LinearIndices{N}; checkargs = true)
+
+Creates the `ArrayOfVectorViews` directly where the user is responsible for having the correct input data.
+Checking of the argument dimensions can be elided by setting `checkargs = false`, but incorrect dimensions
+may lead to illegal out of bounds access later.
+
+`data` is indexed by `indices[i]:indices[i+1]`, where `i = lin_idx[idx...]` and `idx...` are the user-provided
+indices to the `ArrayOfVectorViews`.
+"""
+function ArrayOfVectorViews(indices::Vector{Int}, data::Vector{T}, lin_idx::LinearIndices{N}; checkargs = true) where {T, N}
+    if checkargs
+        checkbounds(data, 1:(last(indices) - 1))
+        checkbounds(indices, last(lin_idx) + 1)
+        issorted(indices) || throw(ArgumentError("indices must be weakly increasing"))
+    end
+    return ArrayOfVectorViews{T, N}(indices, data, lin_idx)
+end
+
+end
diff --git a/src/Ferrite.jl b/src/Ferrite.jl
index ca5ef421e5..0776f2f6be 100644
--- a/src/Ferrite.jl
+++ b/src/Ferrite.jl
@@ -25,9 +25,13 @@ using Tensors:
 using ForwardDiff:
     ForwardDiff
 
+include("CollectionsOfViews.jl")
+using .CollectionsOfViews:
+    CollectionsOfViews, ArrayOfVectorViews, push_at_index!, ConstructionBuffer
 
 include("exports.jl")
 
+
 """
     AbstractRefShape{refdim}
 
diff --git a/src/Grid/topology.jl b/src/Grid/topology.jl
index b1caaf2efc..0d314a89eb 100644
--- a/src/Grid/topology.jl
+++ b/src/Grid/topology.jl
@@ -1,4 +1,3 @@
-
 ############
 # Topology #
 ############
@@ -14,65 +13,90 @@ the given entity is included in the returned list as well.
 """
 getneighborhood
 
-struct EntityNeighborhood{T<:Union{BoundaryIndex,CellIndex}}
-    neighbor_info::Vector{T}
-end
-
-EntityNeighborhood(info::T) where T <: BoundaryIndex = EntityNeighborhood([info])
-Base.length(n::EntityNeighborhood) = length(n.neighbor_info)
-Base.getindex(n::EntityNeighborhood,i) = getindex(n.neighbor_info,i)
-Base.firstindex(n::EntityNeighborhood) = 1
-Base.lastindex(n::EntityNeighborhood) = length(n.neighbor_info)
-Base.:(==)(n1::EntityNeighborhood, n2::EntityNeighborhood) = n1.neighbor_info == n2.neighbor_info
-Base.iterate(n::EntityNeighborhood, state=1) = iterate(n.neighbor_info,state)
-
-function Base.show(io::IO, ::MIME"text/plain", n::EntityNeighborhood)
-    if length(n) == 0
-        println(io, "No EntityNeighborhood")
-    elseif length(n) == 1
-        println(io, "$(n.neighbor_info[1])")
-    else
-        println(io, "$(n.neighbor_info...)")
-    end
-end
 
 abstract type AbstractTopology end
 
 """
-    ExclusiveTopology(cells::Vector{C}) where C <: AbstractCell
-    ExclusiveTopology(grid::Grid)
+    ExclusiveTopology(grid::AbstractGrid)
 
-`ExclusiveTopology` saves topological (connectivity/neighborhood) data of the grid. The constructor works with an `AbstractCell`
-vector for all cells that dispatch `vertices`, `edges`, and `faces`.
-The struct saves the highest dimensional neighborhood, i.e. if something is connected by a face and an
- edge only the face neighborhood is saved. The lower dimensional neighborhood is recomputed, if needed.
+The **experimental feature** `ExclusiveTopology` saves topological (connectivity/neighborhood) data of the grid.
+Only the highest dimensional neighborhood is saved. I.e., if something is connected by a face and an
+edge, only the face neighborhood is saved. The lower dimensional neighborhood is recomputed when calling getneighborhood if needed.
 
 # Fields
-- `vertex_to_cell::Vector{Set{Int}}`: global vertex id to all cells containing the vertex
-- `cell_neighbor::Vector{EntityNeighborhood{CellIndex}}`: cellid to all connected cells
-- `face_neighbor::Matrix{EntityNeighborhood,Int}`: `face_neighbor[cellid,local_face_id]` -> neighboring face
-- `vertex_neighbor::Matrix{EntityNeighborhood,Int}`: `vertex_neighbor[cellid,local_vertex_id]` -> neighboring vertex
-- `edge_neighbor::Matrix{EntityNeighborhood,Int}`: `edge_neighbor[cellid_local_vertex_id]` -> neighboring edge
-- `face_skeleton::Union{Vector{FaceIndex}, Nothing}`:
-
-!!! note Currently mixed-dimensional queries do not work at the moment. They will be added back later.
+- `vertex_to_cell::AbstractArray{AbstractVector{Int}, 1}`:           global vertex id to all cells containing the vertex
+- `cell_neighbor::AbstractArray{AbstractVector{Int}, 1}`:            cellid to all connected cells
+- `face_neighbor::AbstractArray{AbstractVector{FaceIndex}, 2}`:      `face_neighbor[cellid,   local_face_id]`   -> neighboring faces
+- `edge_neighbor::AbstractArray{AbstractVector{EdgeIndex}, 2}`:      `edge_neighbor[cellid,   local_edge_id]`   -> neighboring edges
+- `vertex_neighbor::AbstractArray{AbstractVector{VertexIndex}, 2}`:  `vertex_neighbor[cellid, local_vertex_id]` -> neighboring vertices
+- `face_skeleton::Union{Vector{FaceIndex}, Nothing}`:     List of unique faces in the grid given as `FaceIndex`
+- `edge_skeleton::Union{Vector{EdgeIndex}, Nothing}`:     List of unique edges in the grid given as `EdgeIndex`
+- `vertex_skeleton::Union{Vector{VertexIndex}, Nothing}`: List of unique vertices in the grid given as `VertexIndex`
+
+!!! warning "Limitations"
+    The implementation only works with conforming grids, i.e. grids without "hanging nodes". Non-conforming grids will give unexpected results.
+    Grids with embedded cells (different reference dimension compared
+    to the spatial dimension) are not supported, and will error on construction.
+
 """
 mutable struct ExclusiveTopology <: AbstractTopology
-    # maps a global vertex id to all cells containing the vertex
-    vertex_to_cell::Vector{Set{Int}}
-    # index of the vector = cell id ->  all other connected cells
-    cell_neighbor::Vector{EntityNeighborhood{CellIndex}}
-    # face_neighbor[cellid,local_face_id] -> exclusive connected entities (not restricted to one entity)
-    face_face_neighbor::Matrix{EntityNeighborhood{FaceIndex}}
-    # vertex_neighbor[cellid,local_vertex_id] -> exclusive connected entities to the given vertex
-    vertex_vertex_neighbor::Matrix{EntityNeighborhood{VertexIndex}}
-    # edge_neighbor[cellid,local_edge_id] -> exclusive connected entities of the given edge
-    edge_edge_neighbor::Matrix{EntityNeighborhood{EdgeIndex}}
-    # lazy constructed face topology
+    vertex_to_cell::ArrayOfVectorViews{Int, 1}
+    cell_neighbor::ArrayOfVectorViews{Int, 1}
+    # face_face_neighbor[cellid,local_face_id] -> exclusive connected entities (not restricted to one entity)
+    face_face_neighbor::ArrayOfVectorViews{FaceIndex, 2}
+    # edge_edge_neighbor[cellid,local_edge_id] -> exclusive connected entities of the given edge
+    edge_edge_neighbor::ArrayOfVectorViews{EdgeIndex, 2}
+    # vertex_vertex_neighbor[cellid,local_vertex_id] -> exclusive connected entities to the given vertex
+    vertex_vertex_neighbor::ArrayOfVectorViews{VertexIndex, 2}
+
     face_skeleton::Union{Vector{FaceIndex}, Nothing}
     edge_skeleton::Union{Vector{EdgeIndex}, Nothing}
     vertex_skeleton::Union{Vector{VertexIndex}, Nothing}
-    # TODO reintroduce the codimensional connectivity, e.g. 3D edge to 2D face
+end
+
+function ExclusiveTopology(grid::AbstractGrid{sdim}) where sdim
+    if sdim != get_reference_dimension(grid)
+        error("ExclusiveTopology does not support embedded cells (i.e. reference dimensions different from the spatial dimension)")
+    end
+    cells = getcells(grid)
+    nnodes = getnnodes(grid)
+    ncells = length(cells)
+
+    max_vertices, max_edges, max_faces = _max_nentities_per_cell(cells)
+    vertex_to_cell = build_vertex_to_cell(cells; max_vertices, nnodes)
+    cell_neighbor = build_cell_neighbor(grid, cells, vertex_to_cell; ncells)
+
+    # Here we don't use the convenience constructor taking a function,
+    # since we want to do it simultaneously for 3 data-types
+    facedata = sizehint!(FaceIndex[], ncells * max_faces * _getsizehint(grid, FaceIndex))
+    face_face_neighbor_buf = CollectionsOfViews.ConstructionBuffer(facedata, (ncells, max_faces), _getsizehint(grid, FaceIndex))
+    edgedata = sizehint!(EdgeIndex[], ncells * max_edges * _getsizehint(grid, EdgeIndex))
+    edge_edge_neighbor_buf = CollectionsOfViews.ConstructionBuffer(edgedata, (ncells, max_edges), _getsizehint(grid, EdgeIndex))
+    vertdata = sizehint!(VertexIndex[], ncells * max_vertices * _getsizehint(grid, VertexIndex))
+    vertex_vertex_neighbor_buf = CollectionsOfViews.ConstructionBuffer(vertdata, (ncells, max_vertices), _getsizehint(grid, VertexIndex))
+
+    for (cell_id, cell) in enumerate(cells)
+        for neighbor_cell_id in cell_neighbor[cell_id]
+            neighbor_cell = cells[neighbor_cell_id]
+            getrefdim(neighbor_cell) == getrefdim(cell) || error("Not supported")
+            num_shared_vertices = _num_shared_vertices(cell, neighbor_cell)
+            if num_shared_vertices == 1
+                _add_single_vertex_neighbor!(vertex_vertex_neighbor_buf, cell, cell_id, neighbor_cell, neighbor_cell_id)
+            # Shared edge
+            elseif num_shared_vertices == 2
+                _add_single_edge_neighbor!(edge_edge_neighbor_buf, cell, cell_id, neighbor_cell, neighbor_cell_id)
+            # Shared face
+            elseif num_shared_vertices >= 3
+                _add_single_face_neighbor!(face_face_neighbor_buf, cell, cell_id, neighbor_cell, neighbor_cell_id)
+            else
+                error("Found connected elements without shared vertex... Mesh broken?")
+            end
+        end
+    end
+    face_face_neighbor     = ArrayOfVectorViews(face_face_neighbor_buf)
+    edge_edge_neighbor     = ArrayOfVectorViews(edge_edge_neighbor_buf)
+    vertex_vertex_neighbor = ArrayOfVectorViews(vertex_vertex_neighbor_buf)
+    return ExclusiveTopology(vertex_to_cell, cell_neighbor, face_face_neighbor, edge_edge_neighbor, vertex_vertex_neighbor, nothing, nothing, nothing)
 end
 
 function get_facet_facet_neighborhood(t::ExclusiveTopology, g::AbstractGrid)
@@ -86,12 +110,32 @@ function _get_facet_facet_neighborhood(::ExclusiveTopology, #=rdim=#::Val{:mixed
     Access the `vertex_vertex_neighbor`, `edge_edge_neighbor`, or `face_face_neighbor` fields explicitly instead."))
 end
 
-function Base.show(io::IO, ::MIME"text/plain", topology::ExclusiveTopology)
-    println(io, "ExclusiveTopology\n")
-    print(io, "  Vertex neighbors: $(size(topology.vertex_vertex_neighbor))\n")
-    print(io, "  Face neighbors: $(size(topology.face_face_neighbor))\n")
-    println(io, "  Edge neighbors: $(size(topology.edge_edge_neighbor))")
+# Guess of how many neighbors depending on grid dimension and index type.
+# This could be possible to optimize further by studying connectivities of non-uniform
+# grids, see https://github.com/Ferrite-FEM/Ferrite.jl/pull/974#discussion_r1660838649
+function _getsizehint(g::AbstractGrid, ::Type{IDX}) where IDX
+    CT = getcelltype(g)
+    isconcretetype(CT) && return _getsizehint(getrefshape(CT)(), IDX)
+    rdim = get_reference_dimension(g)::Int
+    return _getsizehint(RefSimplex{rdim}(), IDX) # Simplex is "worst case", used as default.
 end
+_getsizehint(::AbstractRefShape,    ::Type{FaceIndex}) = 1 # Always 1 or zero if not mixed rdim
+
+_getsizehint(::AbstractRefShape{1}, ::Type{EdgeIndex}) = 1
+_getsizehint(::AbstractRefShape{2}, ::Type{EdgeIndex}) = 1
+_getsizehint(::AbstractRefShape{3}, ::Type{EdgeIndex}) = 3 # Number for RefTetrahedron
+_getsizehint(::RefHexahedron,       ::Type{EdgeIndex}) = 1 # Optim for RefHexahedron
+
+_getsizehint(::AbstractRefShape{1}, ::Type{VertexIndex}) = 1
+_getsizehint(::AbstractRefShape{2}, ::Type{VertexIndex}) = 3
+_getsizehint(::AbstractRefShape{3}, ::Type{VertexIndex}) = 13
+_getsizehint(::RefHypercube,        ::Type{VertexIndex}) = 1 # Optim for RefHypercube
+
+_getsizehint(::AbstractRefShape{1}, ::Type{CellIndex}) = 2
+_getsizehint(::AbstractRefShape{2}, ::Type{CellIndex}) = 12
+_getsizehint(::AbstractRefShape{3}, ::Type{CellIndex}) = 70
+_getsizehint(::RefQuadrilateral,    ::Type{CellIndex}) = 8
+_getsizehint(::RefHexahedron,       ::Type{CellIndex}) = 26
 
 function _num_shared_vertices(cell_a::C1, cell_b::C2) where {C1, C2}
     num_shared_vertices = 0
@@ -106,164 +150,125 @@ function _num_shared_vertices(cell_a::C1, cell_b::C2) where {C1, C2}
     return num_shared_vertices
 end
 
-function _exclusive_topology_ctor(cells::Vector{C}, vertex_cell_table::Array{Set{Int}}, vertex_table, face_table, edge_table, cell_neighbor_table) where C <: AbstractCell
-    for (cell_id, cell) in enumerate(cells)
-        # Gather all cells which are connected via vertices
-        cell_neighbor_ids = Set{Int}()
-        for vertex ∈ vertices(cell)
-            for vertex_cell_id ∈ vertex_cell_table[vertex]
-                if vertex_cell_id != cell_id
-                    push!(cell_neighbor_ids, vertex_cell_id)
-                end
-            end
-        end
-        cell_neighbor_table[cell_id] = EntityNeighborhood(CellIndex.(collect(cell_neighbor_ids)))
-
-        # Any of the neighbors is now sorted in the respective categories
-        for cell_neighbor_id ∈ cell_neighbor_ids
-            # Buffer neighbor
-            cell_neighbor = cells[cell_neighbor_id]
-            # TODO handle mixed-dimensional case
-            getrefdim(cell_neighbor) == getrefdim(cell) || continue
-
-            num_shared_vertices = _num_shared_vertices(cell, cell_neighbor)
-
-            # Simplest case: Only one vertex is shared => Vertex neighbor
-            if num_shared_vertices == 1
-                for (lvi, vertex) ∈ enumerate(vertices(cell))
-                    for (lvi2, vertex_neighbor) ∈ enumerate(vertices(cell_neighbor))
-                        if vertex_neighbor == vertex
-                            push!(vertex_table[cell_id, lvi].neighbor_info, VertexIndex(cell_neighbor_id, lvi2))
-                            break
-                        end
-                    end
-                end
-            # Shared edge
-            elseif num_shared_vertices == 2
-                _add_single_edge_neighbor!(edge_table, cell, cell_id, cell_neighbor, cell_neighbor_id)
-            # Shared face
-            elseif num_shared_vertices >= 3
-                _add_single_face_neighbor!(face_table, cell, cell_id, cell_neighbor, cell_neighbor_id)
-            else
-                @error "Found connected elements without shared vertex... Mesh broken?"
-            end
-        end
-    end
-end
-
-function ExclusiveTopology(cells::Vector{C}) where C <: AbstractCell
-    # Setup the cell to vertex table
-    cell_vertices_table = vertices.(cells) #needs generic interface for <: AbstractCell
-    vertex_cell_table = Set{Int}[Set{Int}() for _ ∈ 1:maximum(maximum.(cell_vertices_table))]
-
-    # Setup vertex to cell connectivity by flipping the cell to vertex table
-    for (cellid, cell_vertices) in enumerate(cell_vertices_table)
-        for vertex in cell_vertices
-            push!(vertex_cell_table[vertex], cellid)
-        end
-    end
-
-    # Compute correct matrix size
-    celltype = eltype(cells)
-    max_vertices = 0
-    max_faces = 0
-    max_edges = 0
-    if isconcretetype(celltype)
-        max_vertices = nvertices(cells[1])
-        max_faces = nfaces(cells[1])
-        max_edges = nedges(cells[1])
+"Return the highest number of vertices, edges, and faces per cell"
+function _max_nentities_per_cell(cells::Vector{C}) where C
+    if isconcretetype(C)
+        cell = first(cells)
+        return nvertices(cell), nedges(cell), nfaces(cell)
     else
         celltypes = Set(typeof.(cells))
+        max_vertices = 0
+        max_edges = 0
+        max_faces = 0
         for celltype in celltypes
-            celltypeidx = findfirst(x->typeof(x)==celltype,cells)
-            max_vertices = max(max_vertices,nvertices(cells[celltypeidx]))
-            max_faces = max(max_faces, nfaces(cells[celltypeidx]))
+            celltypeidx = findfirst(x -> isa(x, celltype), cells)
+            max_vertices = max(max_vertices, nvertices(cells[celltypeidx]))
             max_edges = max(max_edges, nedges(cells[celltypeidx]))
+            max_faces = max(max_faces, nfaces(cells[celltypeidx]))
         end
+        return max_vertices, max_edges, max_faces
     end
-
-    # Setup matrices
-    vertex_table = Matrix{EntityNeighborhood{VertexIndex}}(undef, length(cells), max_vertices)
-    for j = 1:size(vertex_table,2)
-        for i = 1:size(vertex_table,1)
-            vertex_table[i,j] = EntityNeighborhood{VertexIndex}(VertexIndex[])
-        end
-    end
-    face_table   = Matrix{EntityNeighborhood{FaceIndex}}(undef, length(cells), max_faces)
-    for j = 1:size(face_table,2)
-        for i = 1:size(face_table,1)
-            face_table[i,j] = EntityNeighborhood{FaceIndex}(FaceIndex[])
-        end
-    end
-    edge_table   = Matrix{EntityNeighborhood{EdgeIndex}}(undef, length(cells), max_edges)
-    for j = 1:size(edge_table,2)
-        for i = 1:size(edge_table,1)
-            edge_table[i,j] = EntityNeighborhood{EdgeIndex}(EdgeIndex[])
-        end
-    end
-    cell_neighbor_table = Vector{EntityNeighborhood{CellIndex}}(undef, length(cells))
-
-    _exclusive_topology_ctor(cells, vertex_cell_table, vertex_table, face_table, edge_table, cell_neighbor_table)
-
-    return ExclusiveTopology(vertex_cell_table,cell_neighbor_table,face_table,vertex_table,edge_table,nothing,nothing,nothing)
 end
 
-function _add_single_face_neighbor!(face_table, cell::C1, cell_id, cell_neighbor::C2, cell_neighbor_id) where {C1, C2}
+function _add_single_face_neighbor!(face_table::ConstructionBuffer, cell::AbstractCell, cell_id::Int, cell_neighbor::AbstractCell, cell_neighbor_id::Int)
     for (lfi, face) ∈ enumerate(faces(cell))
         uniqueface = sortface_fast(face)
         for (lfi2, face_neighbor) ∈ enumerate(faces(cell_neighbor))
             uniqueface2 = sortface_fast(face_neighbor)
             if uniqueface == uniqueface2
-                push!(face_table[cell_id, lfi].neighbor_info, FaceIndex(cell_neighbor_id, lfi2))
+                push_at_index!(face_table, FaceIndex(cell_neighbor_id, lfi2), cell_id, lfi)
                 return
             end
         end
     end
 end
 
-function _add_single_edge_neighbor!(edge_table, cell::C1, cell_id, cell_neighbor::C2, cell_neighbor_id) where {C1, C2}
+function _add_single_edge_neighbor!(edge_table::ConstructionBuffer, cell::AbstractCell, cell_id::Int, cell_neighbor::AbstractCell, cell_neighbor_id::Int)
     for (lei, edge) ∈ enumerate(edges(cell))
         uniqueedge = sortedge_fast(edge)
         for (lei2, edge_neighbor) ∈ enumerate(edges(cell_neighbor))
             uniqueedge2 = sortedge_fast(edge_neighbor)
             if uniqueedge == uniqueedge2
-                push!(edge_table[cell_id, lei].neighbor_info, EdgeIndex(cell_neighbor_id, lei2))
+                push_at_index!(edge_table, EdgeIndex(cell_neighbor_id, lei2), cell_id, lei)
                 return
             end
         end
     end
 end
 
+function _add_single_vertex_neighbor!(vertex_table::ConstructionBuffer, cell::AbstractCell, cell_id::Int, cell_neighbor::AbstractCell, cell_neighbor_id::Int)
+    for (lvi, vertex) ∈ enumerate(vertices(cell))
+        for (lvi2, vertex_neighbor) ∈ enumerate(vertices(cell_neighbor))
+            if vertex_neighbor == vertex
+                push_at_index!(vertex_table, VertexIndex(cell_neighbor_id, lvi2), cell_id, lvi)
+                break
+            end
+        end
+    end
+end
+
+function build_vertex_to_cell(cells; max_vertices, nnodes)
+    vertex_to_cell = ArrayOfVectorViews(sizehint!(Int[], max_vertices * nnodes), (nnodes,); sizehint = max_vertices) do cov
+            for (cellid, cell) in enumerate(cells)
+                for vertex in vertices(cell)
+                    push_at_index!(cov, cellid, vertex)
+                end
+            end
+        end
+    return vertex_to_cell
+end
 
-getcells(neighbor::EntityNeighborhood{T}) where T <: BoundaryIndex = first.(neighbor.neighbor_info)
-getcells(neighbor::EntityNeighborhood{CellIndex}) = getproperty.(neighbor.neighbor_info, :idx)
-getcells(neighbors::Vector{T}) where T <: EntityNeighborhood = reduce(vcat, getcells.(neighbors))
-getcells(neighbors::Vector{T}) where T <: BoundaryIndex = getindex.(neighbors,1)
+function build_cell_neighbor(grid, cells, vertex_to_cell; ncells)
+    # In this case, we loop over the cells in order and all neighbors at once.
+    # Then we can create ArrayOfVectorViews directly without the CollectionsOfViews.ConstructionBuffer
+    sizehint = _getsizehint(grid, CellIndex)
+    data = empty!(Vector{Int}(undef, ncells * sizehint))
 
-ExclusiveTopology(grid::AbstractGrid) = ExclusiveTopology(getcells(grid))
+    indices = Vector{Int}(undef, ncells + 1)
+    cell_neighbor_ids = Int[]
+    n = 1
+    for (cell_id, cell) in enumerate(cells)
+        empty!(cell_neighbor_ids)
+        for vertex ∈ vertices(cell)
+            for vertex_cell_id ∈ vertex_to_cell[vertex]
+                if vertex_cell_id != cell_id
+                    vertex_cell_id ∈ cell_neighbor_ids || push!(cell_neighbor_ids, vertex_cell_id)
+                end
+            end
+        end
+        indices[cell_id] = n
+        append!(data, cell_neighbor_ids)
+        n += length(cell_neighbor_ids)
+    end
+    indices[end] = n
+    sizehint!(data, length(data)) # Tell julia that we won't add more data
+    return ArrayOfVectorViews(indices, data, LinearIndices(1:ncells))
+end
 
 function getneighborhood(top::ExclusiveTopology, grid::AbstractGrid, cellidx::CellIndex, include_self=false)
-    patch = getcells(top.cell_neighbor[cellidx.idx])
+    patch = top.cell_neighbor[cellidx.idx]
     if include_self
-        return [patch; cellidx.idx]
+        return view(push!(collect(patch), cellidx.idx), 1:(length(patch) + 1))
     else
         return patch
     end
 end
 
 function getneighborhood(top::ExclusiveTopology, grid::AbstractGrid, faceidx::FaceIndex, include_self=false)
+    neighbors = top.face_face_neighbor[faceidx[1], faceidx[2]]
     if include_self
-        return [top.face_face_neighbor[faceidx[1],faceidx[2]].neighbor_info; faceidx]
+        return view(push!(collect(neighbors), faceidx), 1:(length(neighbors) + 1))
     else
-        return top.face_face_neighbor[faceidx[1],faceidx[2]].neighbor_info
+        return neighbors
     end
 end
 
 function getneighborhood(top::ExclusiveTopology, grid::AbstractGrid{2}, edgeidx::EdgeIndex, include_self=false)
+    neighbors = top.edge_edge_neighbor[edgeidx[1], edgeidx[2]]
     if include_self
-        return [top.edge_edge_neighbor[edgeidx[1],edgeidx[2]].neighbor_info; edgeidx]
+        return view(push!(collect(neighbors), edgeidx), 1:(length(neighbors) + 1))
     else
-        return top.edge_edge_neighbor[edgeidx[1],edgeidx[2]].neighbor_info
+        return neighbors
     end
 end
 
@@ -279,26 +284,27 @@ function getneighborhood(top::ExclusiveTopology, grid::AbstractGrid, vertexidx::
         !include_self && local_vertex == vertexidx && continue
         push!(self_reference_local, local_vertex)
     end
-    return self_reference_local
+    return view(self_reference_local, 1:length(self_reference_local))
 end
 
 function getneighborhood(top::ExclusiveTopology, grid::AbstractGrid{3}, edgeidx::EdgeIndex, include_self=false)
     cellid, local_edgeidx = edgeidx[1], edgeidx[2]
-    cell_edges = edges(getcells(grid,cellid))
+    cell_edges = edges(getcells(grid, cellid))
     nonlocal_edgeid = cell_edges[local_edgeidx]
-    cell_neighbors = getneighborhood(top,grid,CellIndex(cellid))
+    cell_neighbors = getneighborhood(top, grid, CellIndex(cellid))
     self_reference_local = EdgeIndex[]
     for cellid in cell_neighbors
-        local_neighbor_edgeid = findfirst(x->issubset(x,nonlocal_edgeid),edges(getcells(grid,cellid)))
+        local_neighbor_edgeid = findfirst(x -> issubset(x, nonlocal_edgeid), edges(getcells(grid, cellid)))
         local_neighbor_edgeid === nothing && continue
         local_edge = EdgeIndex(cellid,local_neighbor_edgeid)
         push!(self_reference_local, local_edge)
     end
     if include_self
-        return unique([top.edge_edge_neighbor[cellid, local_edgeidx].neighbor_info; self_reference_local; edgeidx])
+        neighbors = unique([top.edge_edge_neighbor[cellid, local_edgeidx]; self_reference_local; edgeidx])
     else
-        return unique([top.edge_edge_neighbor[cellid, local_edgeidx].neighbor_info; self_reference_local])
+        neighbors = unique([top.edge_edge_neighbor[cellid, local_edgeidx]; self_reference_local])
     end
+    return view(neighbors, 1:length(neighbors))
 end
 
 function getneighborhood(top::ExclusiveTopology, grid::AbstractGrid, facetindex::FacetIndex, include_self=false)
@@ -314,41 +320,42 @@ function _getneighborhood(::Val{:mixed}, args...)
 end
 
 """
-    vertex_star_stencils(top::ExclusiveTopology, grid::Grid) -> Vector{Int, EntityNeighborhood{VertexIndex}}()
+    vertex_star_stencils(top::ExclusiveTopology, grid::Grid) -> AbstractVector{AbstractVector{VertexIndex}}
 Computes the stencils induced by the edge connectivity of the vertices.
 """
 function vertex_star_stencils(top::ExclusiveTopology, grid::Grid)
     cells = grid.cells
-    stencil_table = Dict{Int,EntityNeighborhood{VertexIndex}}()
-    # Vertex Connectivity
-    for (global_vertexid,cellset) ∈ enumerate(top.vertex_to_cell)
-        vertex_neighbors_local = VertexIndex[]
-        for cell ∈ cellset
-            neighbor_boundary = edges(cells[cell])
-            neighbor_connected_faces = neighbor_boundary[findall(x->global_vertexid ∈ x, neighbor_boundary)]
-            this_local_vertex = findfirst(i->toglobal(grid, VertexIndex(cell, i)) == global_vertexid, 1:nvertices(cells[cell]))
-            push!(vertex_neighbors_local, VertexIndex(cell, this_local_vertex))
-            other_vertices = findfirst.(x->x!=global_vertexid,neighbor_connected_faces)
-            any(other_vertices .=== nothing) && continue
-            neighbor_vertices_global = getindex.(neighbor_connected_faces, other_vertices)
-            neighbor_vertices_local = [VertexIndex(cell,local_vertex) for local_vertex ∈ findall(x->x ∈ neighbor_vertices_global, vertices(cells[cell]))]
-            append!(vertex_neighbors_local, neighbor_vertices_local)
+    stencil_table = ArrayOfVectorViews(VertexIndex[], (getnnodes(grid),); sizehint = 10) do buf
+        # Vertex Connectivity
+        for (global_vertexid,cellset) ∈ enumerate(top.vertex_to_cell)
+            for cell ∈ cellset
+                neighbor_boundary = edges(cells[cell])
+                neighbor_connected_faces = neighbor_boundary[findall(x->global_vertexid ∈ x, neighbor_boundary)]
+                this_local_vertex = findfirst(i->toglobal(grid, VertexIndex(cell, i)) == global_vertexid, 1:nvertices(cells[cell]))
+                push_at_index!(buf, VertexIndex(cell, this_local_vertex), global_vertexid)
+                other_vertices = findfirst.(x->x!=global_vertexid,neighbor_connected_faces)
+                any(other_vertices .=== nothing) && continue
+                neighbor_vertices_global = getindex.(neighbor_connected_faces, other_vertices)
+                neighbor_vertices_local = [VertexIndex(cell,local_vertex) for local_vertex ∈ findall(x->x ∈ neighbor_vertices_global, vertices(cells[cell]))]
+                for vertex_index in neighbor_vertices_local
+                    push_at_index!(buf, vertex_index, global_vertexid)
+                end
+            end
         end
-        stencil_table[global_vertexid] =  EntityNeighborhood(vertex_neighbors_local)
     end
     return stencil_table
 end
 
 """
-    getstencil(top::Dict{Int, EntityNeighborhood{VertexIndex}}, grid::AbstractGrid, vertex_idx::VertexIndex) -> EntityNeighborhood{VertexIndex}
+    getstencil(top::ArrayOfVectorViews{VertexIndex, 1}, grid::AbstractGrid, vertex_idx::VertexIndex) -> AbstractVector{VertexIndex}
 Get an iterateable over the stencil members for a given local entity.
 """
-function getstencil(top::Dict{Int, EntityNeighborhood{VertexIndex}}, grid::Grid, vertex_idx::VertexIndex)
-    return top[toglobal(grid, vertex_idx)].neighbor_info
+function getstencil(top::ArrayOfVectorViews{VertexIndex, 1}, grid::Grid, vertex_idx::VertexIndex)
+    return top[toglobal(grid, vertex_idx)]
 end
 
 """
-    _create_skeleton(neighborhood::Matrix{EntityNeighborhood{BI}}) where BI <: Union{FaceIndex, EdgeIndex, VertexIndex}
+    _create_skeleton(neighborhood::AbstractMatrix{AbstractVector{BI}}) where BI <: Union{FaceIndex, EdgeIndex, VertexIndex}
 
 Materializes the skeleton from the `neighborhood` information by returning a `Vector{BI}` with `BI`s describing
 the unique entities in the grid.
@@ -356,18 +363,17 @@ the unique entities in the grid.
 *Example:* With `BI=EdgeIndex`, and an edge between cells and 1 and 2, with vertices 2 and 5, could be described by either
 `EdgeIndex(1, 2)` or `EdgeIndex(2, 4)`, but only one of these will be in the vector returned by this function.
 """
-function _create_skeleton(neighborhood::Matrix{EntityNeighborhood{BI}}) where BI <: Union{FaceIndex, EdgeIndex, VertexIndex}
+function _create_skeleton(neighborhood::ArrayOfVectorViews{BI, 2}) where BI <: Union{FaceIndex, EdgeIndex, VertexIndex}
     i = 1
-    skeleton = Vector{BI}(undef, length(neighborhood) - count(neighbors -> !isempty(neighbors) , neighborhood) ÷ 2)
+    skeleton = Vector{BI}(undef, length(neighborhood) - count(neighbors -> !isempty(neighbors) , values(neighborhood)) ÷ 2)
     for (idx, entity) in pairs(neighborhood)
-        isempty(entity.neighbor_info) || entity.neighbor_info[][1] > idx[1] || continue
+        isempty(entity) || entity[][1] > idx[1] || continue
         skeleton[i] = BI(idx[1], idx[2])
         i += 1
     end
     return skeleton
 end
 
-#TODO: For the specific entities the grid input is unused
 """
     vertexskeleton(top::ExclusiveTopology, ::AbstractGrid) -> Vector{VertexIndex}
 
diff --git a/src/Grid/utils.jl b/src/Grid/utils.jl
index fe96df897c..c684409aaf 100644
--- a/src/Grid/utils.jl
+++ b/src/Grid/utils.jl
@@ -160,7 +160,7 @@ function _create_boundaryset(f::Function, grid::AbstractGrid, top::ExclusiveTopo
     function _makeset(ff_nh)
         set = OrderedSet{BI}()
         for (ff_nh_idx, neighborhood) in pairs(ff_nh)
-            # ff_nh_idx::CartesianIndex into Matrix{<:EntityNeighborhood}
+            # ff_nh_idx::CartesianIndex into AbstractMatrix{AbstractVector{BI}}
             isempty(neighborhood) || continue # Skip any facets with neighbors (not on boundary)
             cell_idx  = ff_nh_idx[1]
             facet_nr = ff_nh_idx[2]
diff --git a/src/iterators.jl b/src/iterators.jl
index 6360b76b63..b37721f4b3 100644
--- a/src/iterators.jl
+++ b/src/iterators.jl
@@ -372,7 +372,7 @@ function Base.iterate(ii::InterfaceIterator{<:Any, <:Grid{sdim}}, state...) wher
         if isempty(neighborhood[facet_a[1], facet_a[2]])
             continue
         end
-        neighbors = neighborhood[facet_a[1], facet_a[2]].neighbor_info
+        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)
diff --git a/test/runtests.jl b/test/runtests.jl
index 3a5e35e708..2a04c9dcef 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -35,6 +35,7 @@ end
 include("test_utils.jl")
 
 # Unit tests
+include("test_collectionsofviews.jl")
 include("test_interpolations.jl")
 include("test_cellvalues.jl")
 include("test_facevalues.jl")
diff --git a/test/test_collectionsofviews.jl b/test/test_collectionsofviews.jl
new file mode 100644
index 0000000000..4bca3ff5a2
--- /dev/null
+++ b/test/test_collectionsofviews.jl
@@ -0,0 +1,46 @@
+@testset "ArrayOfVectorViews" begin
+    # Create a vector sorting integers into bins and check
+    test_ints = rand(0:99, 100)
+    # Create for 3 different sizehints
+    aovs = map([20, 1, 100]) do sh
+        Ferrite.ArrayOfVectorViews(Int[], (10,); sizehint=sh) do buf
+            for v in test_ints
+                idx = 1 + v ÷ 10
+                Ferrite.push_at_index!(buf, v, idx)
+            end
+        end
+    end
+    # Check correct values for the first one
+    for (idx, v) in enumerate(aovs[1])
+        interval = (10 * (idx-1)):(10 * idx - 1)
+        @test all(x -> x ∈ interval, v)
+        @test count(x -> x ∈ interval, test_ints) == length(v)
+    end
+    for aov in aovs
+        @test sum(length, aov; init=0) == length(test_ints)
+    end
+    # Check that the result is independent of sizehint
+    for idx in eachindex(aovs[1])
+        for aov in aovs[2:end]
+            @test aovs[1][idx] == aov[idx]
+        end
+    end
+
+    # Create an array with random tuple containing and place the tuples
+    # according to the values. Check for 2d and 3d arrays.
+    for N in 2:3
+        tvals = [ntuple(i->rand(0:9), N) for _ in 1:1000]
+        aov = Ferrite.ArrayOfVectorViews(NTuple{N, Int}[], (5,5,5)[1:N]; sizehint=10) do buf
+            for v in tvals
+                idx = 1 .+ v .÷ 2
+                Ferrite.push_at_index!(buf, v, idx...)
+            end
+        end
+        @test sum(length, aov; init=0) == length(tvals)
+        for (idx, v) in pairs(aov)
+            intervals = map(i -> (2 * (i-1)):(2 * i - 1), idx.I)
+            @test all(x -> all(map((z, r) -> z ∈ r, x, intervals)), v)
+            @test count(x -> all(map((z, r) -> z ∈ r, x, intervals)), tvals) == length(v)
+        end
+    end
+end
diff --git a/test/test_dofs.jl b/test/test_dofs.jl
index f0417f2d5f..7cd3e47011 100644
--- a/test/test_dofs.jl
+++ b/test/test_dofs.jl
@@ -659,8 +659,8 @@ end
             # test inner coupling
             _check_dofs(K, dh, sdh, cell_idx, coupling, coupling_idx, vdim, [cell_idx], false)
             # test cross-element coupling
-            neighborhood = Ferrite.getrefdim(dh.grid.cells[1]) > 1 ? topology.face_face_neighbor : topology.vertex_vertex_neighbor
-            neighbors = neighborhood[cell_idx, :]
+            neighborhood = Ferrite.get_facet_facet_neighborhood(topology, grid)
+            neighbors = [neighborhood[cell_idx, i] for i in 1:size(neighborhood, 2)]
             _check_dofs(K, dh, sdh, cell_idx, interface_coupling, interface_coupling_idx, vdim, [i[1][1] for i in  neighbors[.!isempty.(neighbors)]], true)
         end
     end
diff --git a/test/test_grid_dofhandler_vtk.jl b/test/test_grid_dofhandler_vtk.jl
index 3ae5d15124..fa78faf033 100644
--- a/test/test_grid_dofhandler_vtk.jl
+++ b/test/test_grid_dofhandler_vtk.jl
@@ -283,19 +283,26 @@ end
 end
 
 @testset "Grid topology" begin
+    # Error paths
+    mixed_rdim_grid = Grid([Triangle((1,2,3)), Line((2,3))], [Node((0.0, 0.0)), Node((1.0, 0.0)), Node((1.0, 1.1))])
+    @test_throws ErrorException ExclusiveTopology(mixed_rdim_grid)
+    top_line = ExclusiveTopology(generate_grid(Line, (3,)))
+    @test_throws ArgumentError Ferrite.get_facet_facet_neighborhood(top_line, mixed_rdim_grid)
+    @test_throws ArgumentError Ferrite.getneighborhood(top_line, mixed_rdim_grid, FacetIndex(1,2))
+
 #
 #      (1) (2) (3) (4)
 #       +---+---+---+
 #
     linegrid = generate_grid(Line,(3,))
     linetopo = ExclusiveTopology(linegrid)
-    @test linetopo.vertex_vertex_neighbor[1,2] == Ferrite.EntityNeighborhood(VertexIndex(2,1))
+    @test linetopo.vertex_vertex_neighbor[1,2] == [VertexIndex(2,1)]
     @test getneighborhood(linetopo, linegrid, VertexIndex(1,2)) == [VertexIndex(2,1)]
-    @test linetopo.vertex_vertex_neighbor[2,1] == Ferrite.EntityNeighborhood(VertexIndex(1,2))
+    @test linetopo.vertex_vertex_neighbor[2,1] == [VertexIndex(1,2)]
     @test getneighborhood(linetopo, linegrid, VertexIndex(2,1)) == [VertexIndex(1,2)]
-    @test linetopo.vertex_vertex_neighbor[2,2] == Ferrite.EntityNeighborhood(VertexIndex(3,1))
+    @test linetopo.vertex_vertex_neighbor[2,2] == [VertexIndex(3,1)]
     @test getneighborhood(linetopo, linegrid, VertexIndex(2,2)) == [VertexIndex(3,1)]
-    @test linetopo.vertex_vertex_neighbor[3,1] == Ferrite.EntityNeighborhood(VertexIndex(2,2))
+    @test linetopo.vertex_vertex_neighbor[3,1] == [VertexIndex(2,2)]
     @test getneighborhood(linetopo, linegrid, VertexIndex(3,1)) == [VertexIndex(2,2)]
     @test getneighborhood(linetopo, linegrid, FacetIndex(3,1)) == getneighborhood(linetopo, linegrid, VertexIndex(3,1))
 
@@ -321,36 +328,38 @@ end
     topology = ExclusiveTopology(quadgrid)
     faceskeleton = Ferrite.edgeskeleton(topology, quadgrid)
     #test vertex neighbors maps cellid and local vertex id to neighbor id and neighbor local vertex id
-    @test topology.vertex_vertex_neighbor[1,3] == Ferrite.EntityNeighborhood(VertexIndex(4,1))
-    @test topology.vertex_vertex_neighbor[2,4] == Ferrite.EntityNeighborhood(VertexIndex(3,2))
+    @test topology.vertex_vertex_neighbor[1,3] == [VertexIndex(4,1)]
+    @test topology.vertex_vertex_neighbor[2,4] == [VertexIndex(3,2)]
     @test Set(getneighborhood(topology, quadgrid, VertexIndex(2,4))) == Set([VertexIndex(1,3), VertexIndex(3,2), VertexIndex(4,1)])
-    @test topology.vertex_vertex_neighbor[3,3] == Ferrite.EntityNeighborhood(VertexIndex(6,1))
-    @test topology.vertex_vertex_neighbor[3,2] == Ferrite.EntityNeighborhood(VertexIndex(2,4))
-    @test topology.vertex_vertex_neighbor[4,1] == Ferrite.EntityNeighborhood(VertexIndex(1,3))
-    @test topology.vertex_vertex_neighbor[4,4] == Ferrite.EntityNeighborhood(VertexIndex(5,2))
-    @test topology.vertex_vertex_neighbor[5,2] == Ferrite.EntityNeighborhood(VertexIndex(4,4))
-    @test topology.vertex_vertex_neighbor[6,1] == Ferrite.EntityNeighborhood(VertexIndex(3,3))
+    @test topology.vertex_vertex_neighbor[3,3] == [VertexIndex(6,1)]
+    @test topology.vertex_vertex_neighbor[3,2] == [VertexIndex(2,4)]
+    @test topology.vertex_vertex_neighbor[4,1] == [VertexIndex(1,3)]
+    @test topology.vertex_vertex_neighbor[4,4] == [VertexIndex(5,2)]
+    @test topology.vertex_vertex_neighbor[5,2] == [VertexIndex(4,4)]
+    @test topology.vertex_vertex_neighbor[6,1] == [VertexIndex(3,3)]
     @test isempty(getneighborhood(topology, quadgrid, VertexIndex(2,2)))
     @test length(getneighborhood(topology, quadgrid, VertexIndex(2,4))) == 3
     #test face neighbor maps cellid and local face id to neighbor id and neighbor local face id
-    @test topology.edge_edge_neighbor[1,2] == Ferrite.EntityNeighborhood(EdgeIndex(2,4))
-    @test topology.edge_edge_neighbor[1,3] == Ferrite.EntityNeighborhood(EdgeIndex(3,1))
-    @test topology.edge_edge_neighbor[2,3] == Ferrite.EntityNeighborhood(EdgeIndex(4,1))
-    @test topology.edge_edge_neighbor[2,4] == Ferrite.EntityNeighborhood(EdgeIndex(1,2))
-    @test topology.edge_edge_neighbor[3,1] == Ferrite.EntityNeighborhood(EdgeIndex(1,3))
-    @test topology.edge_edge_neighbor[3,2] == Ferrite.EntityNeighborhood(EdgeIndex(4,4))
-    @test topology.edge_edge_neighbor[3,3] == Ferrite.EntityNeighborhood(EdgeIndex(5,1))
-    @test topology.edge_edge_neighbor[4,1] == Ferrite.EntityNeighborhood(EdgeIndex(2,3))
-    @test topology.edge_edge_neighbor[4,3] == Ferrite.EntityNeighborhood(EdgeIndex(6,1))
-    @test topology.edge_edge_neighbor[4,4] == Ferrite.EntityNeighborhood(EdgeIndex(3,2))
-    @test topology.edge_edge_neighbor[5,1] == Ferrite.EntityNeighborhood(EdgeIndex(3,3))
-    @test topology.edge_edge_neighbor[5,2] == Ferrite.EntityNeighborhood(EdgeIndex(6,4))
-    @test topology.edge_edge_neighbor[5,3] == Ferrite.EntityNeighborhood(Ferrite.BoundaryIndex[])
-    @test topology.edge_edge_neighbor[5,4] == Ferrite.EntityNeighborhood(Ferrite.BoundaryIndex[])
-    @test topology.edge_edge_neighbor[6,1] == Ferrite.EntityNeighborhood(EdgeIndex(4,3))
-    @test topology.edge_edge_neighbor[6,2] == Ferrite.EntityNeighborhood(Ferrite.BoundaryIndex[])
-    @test topology.edge_edge_neighbor[6,3] == Ferrite.EntityNeighborhood(Ferrite.BoundaryIndex[])
-    @test topology.edge_edge_neighbor[6,4] == Ferrite.EntityNeighborhood(EdgeIndex(5,2))
+    @test topology.edge_edge_neighbor[1,2] == [EdgeIndex(2,4)]
+    @test topology.edge_edge_neighbor[1,3] == [EdgeIndex(3,1)]
+    @test topology.edge_edge_neighbor[2,3] == [EdgeIndex(4,1)]
+    @test topology.edge_edge_neighbor[2,4] == [EdgeIndex(1,2)]
+    @test topology.edge_edge_neighbor[3,1] == [EdgeIndex(1,3)]
+    @test topology.edge_edge_neighbor[3,2] == [EdgeIndex(4,4)]
+    @test topology.edge_edge_neighbor[3,3] == [EdgeIndex(5,1)]
+    @test topology.edge_edge_neighbor[4,1] == [EdgeIndex(2,3)]
+    @test topology.edge_edge_neighbor[4,3] == [EdgeIndex(6,1)]
+    @test topology.edge_edge_neighbor[4,4] == [EdgeIndex(3,2)]
+    @test topology.edge_edge_neighbor[5,1] == [EdgeIndex(3,3)]
+    @test topology.edge_edge_neighbor[5,2] == [EdgeIndex(6,4)]
+    @test isempty(topology.edge_edge_neighbor[5,3])
+    @test isempty(topology.edge_edge_neighbor[5,4])
+    @test topology.edge_edge_neighbor[6,1] == [EdgeIndex(4,3)]
+    @test isempty(topology.edge_edge_neighbor[6,2])
+    @test isempty(topology.edge_edge_neighbor[6,3])
+    @test topology.edge_edge_neighbor[6,4] == [EdgeIndex(5,2)]
+    @test getneighborhood(topology, quadgrid, EdgeIndex(6, 4), false) == [EdgeIndex(5,2)]
+    @test Set(getneighborhood(topology, quadgrid, EdgeIndex(6, 4), true)) == Set([EdgeIndex(6, 4), EdgeIndex(5,2)])
 
     @test getneighborhood(topology, quadgrid, EdgeIndex(2,4)) == getneighborhood(topology, quadgrid, FacetIndex(2,4))
 
@@ -383,14 +392,14 @@ end
 #                        (11)
     hexgrid = generate_grid(Hexahedron,(2,2,1))
     topology = ExclusiveTopology(hexgrid)
-    @test topology.edge_edge_neighbor[1,11] == Ferrite.EntityNeighborhood(EdgeIndex(4,9))
+    @test topology.edge_edge_neighbor[1,11] == [EdgeIndex(4,9)]
     @test Set(getneighborhood(topology,hexgrid,EdgeIndex(1,11),true)) == Set([EdgeIndex(4,9),EdgeIndex(2,12),EdgeIndex(3,10),EdgeIndex(1,11)])
     @test Set(getneighborhood(topology,hexgrid,EdgeIndex(1,11),false)) == Set([EdgeIndex(4,9),EdgeIndex(2,12),EdgeIndex(3,10)])
-    @test topology.edge_edge_neighbor[2,12] == Ferrite.EntityNeighborhood(EdgeIndex(3,10))
+    @test topology.edge_edge_neighbor[2,12] == [EdgeIndex(3,10)]
     @test Set(getneighborhood(topology,hexgrid,EdgeIndex(2,12),true)) == Set([EdgeIndex(3,10),EdgeIndex(1,11),EdgeIndex(4,9),EdgeIndex(2,12)])
     @test Set(getneighborhood(topology,hexgrid,EdgeIndex(2,12),false)) == Set([EdgeIndex(3,10),EdgeIndex(1,11),EdgeIndex(4,9)])
-    @test topology.edge_edge_neighbor[3,10] == Ferrite.EntityNeighborhood(EdgeIndex(2,12))
-    @test topology.edge_edge_neighbor[4,9] == Ferrite.EntityNeighborhood(EdgeIndex(1,11))
+    @test topology.edge_edge_neighbor[3,10] == [EdgeIndex(2,12)]
+    @test topology.edge_edge_neighbor[4,9] == [EdgeIndex(1,11)]
     @test getneighborhood(topology,hexgrid,FaceIndex((1,3))) == [FaceIndex((2,5))]
     @test getneighborhood(topology,hexgrid,FaceIndex((1,4))) == [FaceIndex((3,2))]
     @test getneighborhood(topology,hexgrid,FaceIndex((2,4))) == [FaceIndex((4,2))]
@@ -399,6 +408,7 @@ end
     @test getneighborhood(topology,hexgrid,FaceIndex((3,3))) == [FaceIndex((4,5))]
     @test getneighborhood(topology,hexgrid,FaceIndex((4,2))) == [FaceIndex((2,4))]
     @test getneighborhood(topology,hexgrid,FaceIndex((4,5))) == [FaceIndex((3,3))]
+    @test Set(getneighborhood(topology, hexgrid, FaceIndex((4,5)), true)) == Set([FaceIndex((3,3)), FaceIndex((4,5))])
 
     @test getneighborhood(topology, hexgrid, FaceIndex(2,4)) == getneighborhood(topology, hexgrid, FacetIndex(2,4))
 
@@ -429,7 +439,12 @@ end
 # test for multiple vertex_neighbors as in e.g. ele 3, local vertex 3 (middle node)
     trigrid = generate_grid(Triangle,(2,2))
     topology = ExclusiveTopology(trigrid)
-    @test topology.vertex_vertex_neighbor[3,3] == Ferrite.EntityNeighborhood([VertexIndex(5,2),VertexIndex(6,1),VertexIndex(7,1)])
+    tri_vert_nbset = Set([VertexIndex(5,2), VertexIndex(6,1), VertexIndex(7,1)]) # Exclusive neighbors
+    @test Set(topology.vertex_vertex_neighbor[3,3]) == tri_vert_nbset
+    union!(tri_vert_nbset, [VertexIndex(4, 3), VertexIndex(2, 2)]) # Add vertices shared via edges as well
+    @test Set(getneighborhood(topology, trigrid, VertexIndex(3,3), false)) == tri_vert_nbset
+    union!(tri_vert_nbset, [VertexIndex(3,3)]) # Add self
+    @test Set(getneighborhood(topology, trigrid, VertexIndex(3,3), true))  == tri_vert_nbset
 
     quadtrigrid = generate_grid(QuadraticTriangle,(2,2))
     quadtopology = ExclusiveTopology(trigrid)
@@ -461,7 +476,19 @@ end
         FaceIndex(5,1),                 FaceIndex(5,3), FaceIndex(5,4),
         FaceIndex(6,1),                 FaceIndex(6,3),
     ])
-# test mixed grid
+# test grids with mixed celltypes of same refdim
+# (4)---(5)
+#  |     | \
+#  |  1  |2 \
+# (1)---(2)-(3)
+    tet_quad_grid = Grid(
+        [Quadrilateral((1, 2, 5, 4)), Triangle((2, 3, 5))],
+        [Node(Float64.((x, y))) for (x, y) in ((0,0), (1,0), (2,0), (0,1), (1,1))])
+    top = ExclusiveTopology(tet_quad_grid)
+    @test getneighborhood(top, tet_quad_grid, FacetIndex(1, 2)) == [EdgeIndex(2,3)]
+    @test Set(getneighborhood(top, tet_quad_grid, FacetIndex(1, 2), true)) == Set([EdgeIndex(1,2), EdgeIndex(2,3)])
+
+# test grids with mixed refdim
     cells = [
         Hexahedron((1, 2, 3, 4, 5, 6, 7, 8)),
         Hexahedron((11, 13, 14, 12, 15, 16, 17, 18)),
@@ -470,11 +497,12 @@ end
         ]
     nodes = [Node(coord) for coord in zeros(Vec{3,Float64}, 18)]
     grid = Grid(cells, nodes)
-    topology = ExclusiveTopology(grid)
+    @test_throws ErrorException ExclusiveTopology(grid)
+    # topology = ExclusiveTopology(grid)
 
-    @test_throws ArgumentError Ferrite.facetskeleton(topology, grid)
-    # @test topology.face_face_neighbor[3,4] == Ferrite.EntityNeighborhood(EdgeIndex(1,2))
-    # @test topology.edge_edge_neighbor[1,2] == Ferrite.EntityNeighborhood(FaceIndex(3,4))
+    # @test_throws ArgumentError Ferrite.facetskeleton(topology, grid)
+    # @test topology.face_face_neighbor[3,4] == [EdgeIndex(1,2))
+    # @test topology.edge_edge_neighbor[1,2] == [FaceIndex(3,4))
     # # regression that it doesn't error for boundary faces, see https://github.com/Ferrite-FEM/Ferrite.jl/issues/518
     # @test topology.face_face_neighbor[1,6] == topology.face_face_neighbor[1,1] == zero(Ferrite.EntityNeighborhood{FaceIndex})
     # @test topology.edge_edge_neighbor[1,1] == topology.edge_edge_neighbor[1,3] == zero(Ferrite.EntityNeighborhood{FaceIndex})
@@ -489,10 +517,11 @@ end
         ]
     nodes = [Node(coord) for coord in zeros(Vec{2,Float64}, 18)]
     grid = Grid(cells, nodes)
-    topology = ExclusiveTopology(grid)
-    @test_throws ArgumentError Ferrite.facetskeleton(topology, grid)
-    @test_throws ArgumentError getneighborhood(topology, grid, FacetIndex(1,1))
-    @test_throws ArgumentError Ferrite.get_facet_facet_neighborhood(topology, grid)
+    @test_throws ErrorException ExclusiveTopology(grid)
+    # topology = ExclusiveTopology(grid)
+    # @test_throws ArgumentError Ferrite.facetskeleton(topology, grid)
+    # @test_throws ArgumentError getneighborhood(topology, grid, FacetIndex(1,1))
+    # @test_throws ArgumentError Ferrite.get_facet_facet_neighborhood(topology, grid)
 
 #
 #                   +-----+-----+-----+
@@ -516,15 +545,17 @@ end
     @test issubset([4,5,8], patches[7])
     @test issubset([7,4,5,6,9], patches[8])
     @test issubset([8,5,6], patches[9])
+    @test 5 ∉ getneighborhood(topology, quadgrid, CellIndex(5))
+    @test 5 ∈ getneighborhood(topology, quadgrid, CellIndex(5), true)
 
 # test star stencils
     stars = Ferrite.vertex_star_stencils(topology, quadgrid)
     @test Set(Ferrite.getstencil(stars, quadgrid, VertexIndex(1,1))) == Set([VertexIndex(1,2), VertexIndex(1,4), VertexIndex(1,1)])
     @test Set(Ferrite.getstencil(stars, quadgrid, VertexIndex(2,1))) == Set([VertexIndex(1,1), VertexIndex(1,3), VertexIndex(2,2), VertexIndex(2,4), VertexIndex(1,2), VertexIndex(2,1)])
     @test Set(Ferrite.getstencil(stars, quadgrid, VertexIndex(5,4))) == Set([VertexIndex(4,2), VertexIndex(4,4), VertexIndex(5,1), VertexIndex(5,3), VertexIndex(7,1), VertexIndex(7,3), VertexIndex(8,2), VertexIndex(8,4), VertexIndex(4,3), VertexIndex(5,4), VertexIndex(7,2), VertexIndex(8,1)])
-    @test Set(Ferrite.toglobal(quadgrid, Ferrite.getstencil(stars, quadgrid, VertexIndex(1,1)))) == Set([1,2,5])
-    @test Set(Ferrite.toglobal(quadgrid, Ferrite.getstencil(stars, quadgrid, VertexIndex(2,1)))) == Set([2,1,6,3])
-    @test Set(Ferrite.toglobal(quadgrid, Ferrite.getstencil(stars, quadgrid, VertexIndex(5,4)))) == Set([10,6,9,11,14])
+    @test Set(Ferrite.toglobal.((quadgrid,), Ferrite.getstencil(stars, quadgrid, VertexIndex(1,1)))) == Set([1,2,5])
+    @test Set(Ferrite.toglobal.((quadgrid,), Ferrite.getstencil(stars, quadgrid, VertexIndex(2,1)))) == Set([2,1,6,3])
+    @test Set(Ferrite.toglobal.((quadgrid,), Ferrite.getstencil(stars, quadgrid, VertexIndex(5,4)))) == Set([10,6,9,11,14])
 
     face_skeleton = Ferrite.edgeskeleton(topology, quadgrid)
     @test Set(face_skeleton) == Set([EdgeIndex(1,1),EdgeIndex(1,2),EdgeIndex(1,3),EdgeIndex(1,4),