diff --git a/docs/src/reference/grid.md b/docs/src/reference/grid.md
index 38d7f7791b..8cf35e93ab 100644
--- a/docs/src/reference/grid.md
+++ b/docs/src/reference/grid.md
@@ -46,6 +46,8 @@ get_cell_coordinates!
 Ferrite.ExclusiveTopology
 Ferrite.getneighborhood
 Ferrite.faceskeleton
+Ferrite.vertex_star_stencils
+Ferrite.getstencil
 ```
 
 ### Grid Sets Utility
diff --git a/src/Dofs/DofHandler.jl b/src/Dofs/DofHandler.jl
index 050ba1e201..f22dbe8c3a 100644
--- a/src/Dofs/DofHandler.jl
+++ b/src/Dofs/DofHandler.jl
@@ -620,6 +620,17 @@ function sortedge(edge::Tuple{Int,Int})
     a < b ? (return (edge, PathOrientationInfo(true))) : (return ((b, a), PathOrientationInfo(false)))
 end
 
+"""
+sortedge_fast(edge::Tuple{Int,Int})
+
+Returns the unique representation of an edge.
+Here the unique representation is the sorted node index tuple.
+"""
+function sortedge_fast(edge::Tuple{Int,Int})
+    a, b = edge
+    a < b ? (return edge) : (return (b, a))
+end
+
 """
     sortface(face::Tuple{Int})
     sortface(face::Tuple{Int,Int})
@@ -633,6 +644,20 @@ so the unique representation is always a tuple length 3.
 """
 sortface(face::Tuple{Int,Int}) = sortedge(face) # Face in 2D is the same as edge in 3D.
 
+
+"""
+    sortface_fast(face::Tuple{Int})
+    sortface_fast(face::Tuple{Int,Int})
+    sortface_fast(face::Tuple{Int,Int,Int})
+    sortface_fast(face::Tuple{Int,Int,Int,Int})
+
+Returns the unique representation of a face.
+Here the unique representation is the sorted node index tuple.
+Note that in 3D we only need indices to uniquely identify a face,
+so the unique representation is always a tuple length 3.
+"""
+sortface_fast(face::Tuple{Int,Int}) = sortedge_fast(face) # Face in 2D is the same as edge in 3D.
+
 """
     !!!NOTE TODO implement me.
 
@@ -669,6 +694,16 @@ function sortface(face::Tuple{Int,Int,Int})
     return (a, b, c), SurfaceOrientationInfo() # TODO fill struct
 end
 
+
+function sortface_fast(face::Tuple{Int,Int,Int})
+    a, b, c = face
+    b, c = minmax(b, c)
+    a, c = minmax(a, c)
+    a, b = minmax(a, b)
+    return (a, b, c)
+end
+
+
 function sortface(face::Tuple{Int,Int,Int,Int})
     a, b, c, d = face
     c, d = minmax(c, d)
@@ -680,7 +715,21 @@ function sortface(face::Tuple{Int,Int,Int,Int})
     return (a, b, c), SurfaceOrientationInfo() # TODO fill struct
 end
 
+
+function sortface_fast(face::Tuple{Int,Int,Int,Int})
+    a, b, c, d = face
+    c, d = minmax(c, d)
+    b, d = minmax(b, d)
+    a, d = minmax(a, d)
+    b, c = minmax(b, c)
+    a, c = minmax(a, c)
+    a, b = minmax(a, b)
+    return (a, b, c)
+end
+
+
 sortface(face::Tuple{Int}) = face, nothing
+sortface_fast(face::Tuple{Int}) = face
 
 """
     find_field(dh::DofHandler, field_name::Symbol)::NTuple{2,Int}
diff --git a/src/Ferrite.jl b/src/Ferrite.jl
index 4cdd8c77ce..533f90f05e 100644
--- a/src/Ferrite.jl
+++ b/src/Ferrite.jl
@@ -89,6 +89,8 @@ include("FEValues/face_integrals.jl")
 
 # Grid
 include("Grid/grid.jl")
+include("Grid/topology.jl")
+include("Grid/utils.jl")
 include("Grid/grid_generators.jl")
 include("Grid/coloring.jl")
 
diff --git a/src/Grid/grid.jl b/src/Grid/grid.jl
index 4e804e317e..d21e0962b5 100644
--- a/src/Grid/grid.jl
+++ b/src/Grid/grid.jl
@@ -238,41 +238,6 @@ struct SerendipityQuadraticHexahedron    <: AbstractCell{RefHexahedron}    nodes
 default_interpolation(::Type{SerendipityQuadraticQuadrilateral}) = Serendipity{RefQuadrilateral, 2}()
 default_interpolation(::Type{SerendipityQuadraticHexahedron})    = Serendipity{RefHexahedron,    2}()
 
-
-############
-# Topology #
-############
-# TODO: Move topology stuff to src/Grid/Topology.jl or something
-
-struct EntityNeighborhood{T<:Union{BoundaryIndex,CellIndex}}
-    neighbor_info::Vector{T}
-end
-
-EntityNeighborhood(info::T) where T <: BoundaryIndex = EntityNeighborhood([info])
-Base.zero(::Type{EntityNeighborhood{T}}) where T = EntityNeighborhood(T[])
-Base.zero(::Type{EntityNeighborhood}) = EntityNeighborhood(BoundaryIndex[])
-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.:+(n1::EntityNeighborhood, n2::EntityNeighborhood)
-    neighbor_info = [n1.neighbor_info; n2.neighbor_info]
-    return EntityNeighborhood(neighbor_info)
-end
-
-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
-
 """
     nvertices_on_face(cell::AbstractCell, local_face_index::Int)
 Specifies for each subtype of AbstractCell how many nodes form a face.
@@ -286,231 +251,12 @@ nvertices_on_edge(cell::AbstractCell, local_edge_index::Int) = length(edges(cell
 
 getdim(::Union{AbstractCell{refshape},Type{<:AbstractCell{refshape}}}) where {refdim, refshape <: AbstractRefShape{refdim}} = refdim
 
-abstract type AbstractTopology end
-
-"""
-    ExclusiveTopology(cells::Vector{C}) where C <: AbstractCell
-`ExclusiveTopology` saves topological (connectivity) data of the grid. The constructor works with an `AbstractCell`
-vector for all cells that dispatch `vertices`, `faces` and in 3D `edges`.
-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.
-
-# Fields
-- `vertex_to_cell::Dict{Int,Vector{Int}}`: global vertex id to all cells containing the vertex
-- `cell_neighbor::Vector{EntityNeighborhood{CellIndex}}`: cellid to all connected cells
-- `face_neighbor::SparseMatrixCSC{EntityNeighborhood,Int}`: `face_neighbor[cellid,local_face_id]` -> neighboring face
-- `vertex_neighbor::SparseMatrixCSC{EntityNeighborhood,Int}`: `vertex_neighbor[cellid,local_vertex_id]` -> neighboring vertex
-- `edge_neighbor::SparseMatrixCSC{EntityNeighborhood,Int}`: `edge_neighbor[cellid_local_vertex_id]` -> neighboring edge
-- `vertex_vertex_neighbor::Dict{Int,EntityNeighborhood{VertexIndex}}`: global vertex id -> all connected vertices by edge or face
-- `face_skeleton::Vector{FaceIndex}`: list of unique faces in the grid 
-"""
-struct ExclusiveTopology <: AbstractTopology
-    # maps a global vertex id to all cells containing the vertex
-    vertex_to_cell::Dict{Int,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_neighbor::SparseMatrixCSC{EntityNeighborhood,Int}
-    # vertex_neighbor[cellid,local_vertex_id] -> exclusive connected entities to the given vertex
-    vertex_neighbor::SparseMatrixCSC{EntityNeighborhood,Int}
-    # edge_neighbor[cellid,local_edge_id] -> exclusive connected entities of the given edge
-    edge_neighbor::SparseMatrixCSC{EntityNeighborhood,Int}
-    # maps global vertex id to all directly (by edge or face) connected vertices (no diagonal connection considered)
-    vertex_vertex_neighbor::Dict{Int,EntityNeighborhood{VertexIndex}}
-    # list of unique faces in the grid given as FaceIndex
-    face_skeleton::Vector{FaceIndex}
-end
-
-function ExclusiveTopology(cells::Vector{C}) where C <: AbstractCell
-    cell_vertices_table = vertices.(cells) #needs generic interface for <: AbstractCell
-    vertex_cell_table = Dict{Int,Set{Int}}() 
-
-    for (cellid, cell_vertices) in enumerate(cell_vertices_table)
-       for vertex in cell_vertices
-            if haskey(vertex_cell_table, vertex)
-                push!(vertex_cell_table[vertex], cellid)
-            else
-                vertex_cell_table[vertex] = Set([cellid])
-            end
-        end 
-    end
-
-    I_face = Int[]; J_face = Int[]; V_face = EntityNeighborhood[]
-    I_edge = Int[]; J_edge = Int[]; V_edge = EntityNeighborhood[]
-    I_vertex = Int[]; J_vertex = Int[]; V_vertex = EntityNeighborhood[]   
-    cell_neighbor_table = Vector{EntityNeighborhood{CellIndex}}(undef, length(cells)) 
-
-    for (cellid, cell) in enumerate(cells)
-        cell_neighbors = reduce(union!, [Set{Int}(vertex_cell_table[vertex]) for vertex ∈ vertices(cell) if vertex_cell_table[vertex] != cellid])
-        cell_neighbor_table[cellid] = EntityNeighborhood(CellIndex.(collect(cell_neighbors)))
-
-        face_neighbors = Set{Int}()
-        for (face_idx,face) ∈ enumerate(faces(cell))
-            neighbor_candidates = Set{Int}(c for c ∈ vertex_cell_table[face[1]] if c != cellid)
-            for face_vertex ∈ face[2:end]
-                intersect!(neighbor_candidates, vertex_cell_table[face_vertex])
-            end
-            union!(face_neighbors, neighbor_candidates)
-        end
-
-        if getdim(cell) > 2
-            edge_neighbors = Set{Int}()
-            for (edge_idx,edge) ∈ enumerate(edges(cell))
-                neighbor_candidates = Set{Int}(c for c ∈ vertex_cell_table[edge[1]] if c != cellid)
-                for edge_vertex ∈ edge[2:end]
-                    edge_neighbor = vertex_cell_table[edge_vertex]
-                    if edge_neighbor != cellid && edge_neighbor ∉ face_neighbors
-                        intersect!(neighbor_candidates, edge_neighbor)
-                    end
-                end
-                union!(edge_neighbors, neighbor_candidates)
-            end
-        end
-
-        for neighbor_cellid in cell_neighbors
-            cell_local_ids = findall(x->x in cell_vertices_table[neighbor_cellid], cell_vertices_table[cellid])
-            # vertex neighbor
-            if length(cell_local_ids) == 1
-                neighbor_local_ids = findall(x->x in cell_vertices_table[cellid], cell_vertices_table[neighbor_cellid])
-                _vertex_neighbor!(V_vertex, I_vertex, J_vertex, cellid, cell, neighbor_local_ids, neighbor_cellid, cells[neighbor_cellid])
-            # face neighbor
-            elseif neighbor_cellid ∈ face_neighbors
-                neighbor_local_ids = findall(x->x in cell_vertices_table[cellid], cell_vertices_table[neighbor_cellid])
-                _face_neighbor!(V_face, I_face, J_face, cellid, cell, neighbor_local_ids, neighbor_cellid, cells[neighbor_cellid]) 
-            # edge neighbor
-            elseif getdim(cell) > 2 && neighbor_cellid ∈ edge_neighbors
-                neighbor_local_ids = findall(x->x in cell_vertices_table[cellid], cell_vertices_table[neighbor_cellid])
-                _edge_neighbor!(V_edge, I_edge, J_edge, cellid, cell, neighbor_local_ids, neighbor_cellid, cells[neighbor_cellid])
-            end
-        end   
-    end
-
-    celltype = eltype(cells)
-    if isconcretetype(celltype)
-        dim = getdim(cells[1])
-        _nvertices = nvertices(cells[1])
-        push!(V_vertex,zero(EntityNeighborhood{VertexIndex}))
-        push!(I_vertex,1); push!(J_vertex,_nvertices)
-        if dim > 1
-            _nfaces = nfaces(cells[1])
-            push!(V_face,zero(EntityNeighborhood{FaceIndex}))
-            push!(I_face,1); push!(J_face,_nfaces)
-        end
-        if dim > 2
-            _nedges = nedges(cells[1])
-            push!(V_edge,zero(EntityNeighborhood{EdgeIndex}))
-            push!(I_edge,1); push!(J_edge,_nedges)
-        end
-    else
-        celltypes = typeof.(cells) 
-        for celltype in celltypes
-            celltypeidx = findfirst(x->typeof(x)==celltype,cells)
-            dim = getdim(cells[celltypeidx])
-            _nvertices = nvertices(cells[celltypeidx])
-            push!(V_vertex,zero(EntityNeighborhood{VertexIndex}))
-            push!(I_vertex,celltypeidx); push!(J_vertex,_nvertices)
-            if dim > 1
-                _nfaces = nfaces(cells[celltypeidx])
-                push!(V_face,zero(EntityNeighborhood{FaceIndex}))
-                push!(I_face,celltypeidx); push!(J_face,_nfaces)
-            end
-            if dim > 2
-                _nedges = nedges(cells[celltypeidx])
-                push!(V_edge,zero(EntityNeighborhood{EdgeIndex}))
-                push!(I_edge,celltypeidx); push!(J_edge,_nedges)
-            end
-        end
-    end
-    face_neighbor = sparse(I_face,J_face,V_face)
-    vertex_neighbor = sparse(I_vertex,J_vertex,V_vertex) 
-    edge_neighbor = sparse(I_edge,J_edge,V_edge)
-
-    vertex_vertex_table = Dict{Int,EntityNeighborhood}()
-    vertex_vertex_global = Dict{Int,Vector{Int}}()
-    # Vertex Connectivity
-    for global_vertexid in keys(vertex_cell_table)
-        #Cellset that contains given vertex 
-        cellset = vertex_cell_table[global_vertexid]
-        vertex_neighbors_local = VertexIndex[]
-        vertex_neighbors_global = Int[]
-        for cell in cellset
-            neighbor_boundary = getdim(cells[cell]) > 2 ? collect(edges(cells[cell])) : collect(faces(cells[cell])) #get lowest dimension boundary
-            neighbor_connected_faces = neighbor_boundary[findall(x->global_vertexid in x, neighbor_boundary)]
-            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 in findall(x->x in neighbor_vertices_global, vertices(cells[cell]))]
-            append!(vertex_neighbors_local, neighbor_vertices_local)
-            append!(vertex_neighbors_global, neighbor_vertices_global)
-        end
-        vertex_vertex_table[global_vertexid] =  EntityNeighborhood(vertex_neighbors_local)
-        vertex_vertex_global[global_vertexid] = vertex_neighbors_global
-    end 
-
-    # Face Skeleton
-    face_skeleton_global = Set{NTuple}()
-    face_skeleton_local = Vector{FaceIndex}()
-    fs_length = length(face_skeleton_global)
-    for (cellid,cell) in enumerate(cells)
-        for (local_face_id,face) in enumerate(faces(cell))
-            push!(face_skeleton_global, first(sortface(face)))
-            fs_length_new = length(face_skeleton_global)
-            if fs_length != fs_length_new
-                push!(face_skeleton_local, FaceIndex(cellid,local_face_id)) 
-                fs_length = fs_length_new
-            end
-        end
-    end
-    return ExclusiveTopology(vertex_cell_table,cell_neighbor_table,face_neighbor,vertex_neighbor,edge_neighbor,vertex_vertex_table,face_skeleton_local)
-end
-
-function _vertex_neighbor!(V_vertex, I_vertex, J_vertex, cellid, cell, neighbor, neighborid, neighbor_cell)
-    vertex_neighbor = VertexIndex((neighborid, neighbor[1]))
-    cell_vertex_id = findfirst(x->x==neighbor_cell.nodes[neighbor[1]], cell.nodes)
-    push!(V_vertex,EntityNeighborhood(vertex_neighbor))
-    push!(I_vertex,cellid)
-    push!(J_vertex,cell_vertex_id)
-end
-
-function _edge_neighbor!(V_edge, I_edge, J_edge, cellid, cell, neighbor, neighborid, neighbor_cell)
-    neighbor_edge = neighbor_cell.nodes[neighbor]
-    if getdim(neighbor_cell) < 3
-        neighbor_edge_id = findfirst(x->issubset(x,neighbor_edge), faces(neighbor_cell))
-        edge_neighbor = FaceIndex((neighborid, neighbor_edge_id))
-    else
-        neighbor_edge_id = findfirst(x->issubset(x,neighbor_edge), edges(neighbor_cell))
-        edge_neighbor = EdgeIndex((neighborid, neighbor_edge_id))
-    end
-    cell_edge_id = findfirst(x->issubset(x,neighbor_edge),edges(cell))
-    push!(V_edge, EntityNeighborhood(edge_neighbor))
-    push!(I_edge, cellid)
-    push!(J_edge, cell_edge_id)
-end
-
-function _face_neighbor!(V_face, I_face, J_face, cellid, cell, neighbor, neighborid, neighbor_cell)
-    neighbor_face = neighbor_cell.nodes[neighbor]
-    if getdim(neighbor_cell) == getdim(cell)
-        neighbor_face_id = findfirst(x->issubset(x,neighbor_face), faces(neighbor_cell))
-        face_neighbor = FaceIndex((neighborid, neighbor_face_id))
-    else
-        neighbor_face_id = findfirst(x->issubset(x,neighbor_face), edges(neighbor_cell))
-        face_neighbor = EdgeIndex((neighborid, neighbor_face_id))
-    end
-    cell_face_id = findfirst(x->issubset(x,neighbor_face),faces(cell))
-    push!(V_face, EntityNeighborhood(face_neighbor))
-    push!(I_face, cellid)
-    push!(J_face, cell_face_id)
-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)
 
+######################
+### Mesh interface ###
+######################
 abstract type AbstractGrid{dim} end
 
-ExclusiveTopology(grid::AbstractGrid) = ExclusiveTopology(getcells(grid))
-
 """
     Grid{dim, C<:AbstractCell, T<:Real} <: AbstractGrid}
 
@@ -562,79 +308,6 @@ Get the datatype for a single point in the grid.
 """
 get_coordinate_type(grid::Grid{dim,C,T}) where {dim,C,T} = Vec{dim,T} # Node is baked into the mesh type.
 
-"""
-    getneighborhood(top::ExclusiveTopology, grid::AbstractGrid, cellidx::CellIndex, include_self=false)
-    getneighborhood(top::ExclusiveTopology, grid::AbstractGrid, faceidx::FaceIndex, include_self=false)
-    getneighborhood(top::ExclusiveTopology, grid::AbstractGrid, vertexidx::VertexIndex, include_self=false)
-    getneighborhood(top::ExclusiveTopology, grid::AbstractGrid, edgeidx::EdgeIndex, include_self=false)
-
-Returns all directly connected entities of the same type, i.e. calling the function with a `VertexIndex` will return
-a list of directly connected vertices (connected via face/edge). If `include_self` is true, the given `*Index` is included 
-in the returned list.
-
-!!! warning
-    This feature is highly experimental and very likely subjected to interface changes in the future.
-"""
-function getneighborhood(top::ExclusiveTopology, grid::AbstractGrid, cellidx::CellIndex, include_self=false)
-    patch = getcells(top.cell_neighbor[cellidx.idx])
-    if include_self
-        return [patch; cellidx.idx]
-    else 
-        return patch
-    end
-end
-
-function getneighborhood(top::ExclusiveTopology, grid::AbstractGrid, faceidx::FaceIndex, include_self=false)
-    if include_self 
-        return [top.face_neighbor[faceidx[1],faceidx[2]].neighbor_info; faceidx]
-    else
-        return top.face_neighbor[faceidx[1],faceidx[2]].neighbor_info
-    end
-end
-
-function getneighborhood(top::ExclusiveTopology, grid::AbstractGrid, vertexidx::VertexIndex, include_self=false)
-    cellid, local_vertexid = vertexidx[1], vertexidx[2]
-    cell_vertices = vertices(getcells(grid,cellid))
-    global_vertexid = cell_vertices[local_vertexid]
-    if include_self
-        vertex_to_cell = top.vertex_to_cell[global_vertexid]
-        self_reference_local = Vector{VertexIndex}(undef,length(vertex_to_cell))
-        for (i,cellid) in enumerate(vertex_to_cell)
-            local_vertex = VertexIndex(cellid,findfirst(x->x==global_vertexid,vertices(getcells(grid,cellid))))
-            self_reference_local[i] = local_vertex
-        end
-        return [top.vertex_vertex_neighbor[global_vertexid].neighbor_info; self_reference_local]
-    else
-        return top.vertex_vertex_neighbor[global_vertexid].neighbor_info
-    end
-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))
-    nonlocal_edgeid = cell_edges[local_edgeidx] 
-    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 === nothing && continue
-        local_edge = EdgeIndex(cellid,local_neighbor_edgeid)
-        push!(self_reference_local, local_edge)
-    end
-    if include_self  
-        return unique([top.edge_neighbor[cellid, local_edgeidx].neighbor_info; self_reference_local; edgeidx])
-    else
-        return unique([top.edge_neighbor[cellid, local_edgeidx].neighbor_info; self_reference_local])
-    end
-end
-
-"""
-    faceskeleton(grid) -> Vector{FaceIndex}
-Returns an iterateable face skeleton. The skeleton consists of `FaceIndex` that can be used to `reinit` 
-`FaceValues`.
-"""
-faceskeleton(top::ExclusiveTopology, grid::AbstractGrid) =  top.face_skeleton
-
 """
     toglobal(grid::AbstractGrid, vertexidx::VertexIndex) -> Int
     toglobal(grid::AbstractGrid, vertexidx::Vector{VertexIndex}) -> Vector{Int}
@@ -945,349 +618,6 @@ for (func,             entity_f, subentity_f, entity_t,   subentity_t) in (
     end
 end
 
-"""
-    getfaceinstances(grid::AbstractGrid, topology::ExclusiveTopology, face::FaceIndex)
-
-Returns all the faces as `Set{FaceIndex}` that share all their vertices with a given face
-represented as `FaceIndex`. The returned set includes the input face.
-
-```julia-repl
-julia> using Ferrite; using Ferrite: getfaceinstances
-
-julia> grid = generate_grid(Tetrahedron, (2,1,1));
-
-julia> topology = ExclusiveTopology(grid);
-
-julia> getfaceinstances(grid, topology, FaceIndex(4,2))
-Set{FaceIndex} with 2 elements:
-  FaceIndex((6, 4))
-  FaceIndex((4, 2))
-```
-"""
-function getfaceinstances end
-
-"""
-    getedgeinstances(grid::AbstractGrid, topology::ExclusiveTopology, edge::EdgeIndex)
-
-Returns all the edges as `Set{EdgeIndex}` that share all their vertices with a given edge
-represented as `EdgeIndex`.
-The returned set includes the input edge.
-
-```julia-repl
-julia> using Ferrite; using Ferrite: getedgeinstances
-
-julia> grid = generate_grid(Tetrahedron, (2,1,1));
-
-julia> topology = ExclusiveTopology(grid);
-
-julia> getedgeinstances(grid, topology, EdgeIndex(4,2))
-Set{EdgeIndex} with 3 elements:
-  EdgeIndex((4, 2))
-  EdgeIndex((9, 6))
-  EdgeIndex((7, 6))
-```
-"""
-function getedgeinstances end
-
-"""
-    getvertexinstances(grid::AbstractGrid, topology::ExclusiveTopology, vertex::EdgeIndex)
-
-Returns all the vertices as `Set{::VertexIndex}` that use a given vertex represented as
-`VertexIndex` in all cells.
-The returned set includes the input vertex.
-
-```julia-repl
-julia> using Ferrite; using Ferrite: getvertexinstances
-
-julia> grid = generate_grid(Tetrahedron,(2,1,1));
-
-julia> topology = ExclusiveTopology(grid);
-
-julia> getvertexinstances(grid, topology, VertexIndex(4,2))
-Set{VertexIndex} with 8 elements:
-  VertexIndex((7, 4))
-  VertexIndex((10, 4))
-  VertexIndex((12, 4))
-  VertexIndex((6, 3))
-  VertexIndex((4, 2))
-  VertexIndex((9, 4))
-  VertexIndex((11, 4))
-  VertexIndex((8, 4))
-```
-"""
-function getvertexinstances end
-
-for (func,                entity_f,  entity_t) in (
-    (:getvertexinstances, :vertices, :VertexIndex),
-    (:getedgeinstances,   :edges,    :EdgeIndex),
-    (:getfaceinstances,   :faces,    :FaceIndex),
-)
-    @eval begin
-      function $(func)(grid::AbstractGrid, topology::ExclusiveTopology, entity::$(entity_t))
-        _set = Set{$(entity_t)}()
-        cells = getcells(grid)
-        cell = cells[entity[1]]
-        verts = $(entity_f)(cell)[entity[2]]
-        # Since we are looking for an entity that share *all* vertices, the first one can be
-        # used here to query potiential neighbor cells
-        for cell_idx in topology.vertex_to_cell[verts[1]] # Since all vertices should be shared, the first one can be used here
-            cell_entities = $(entity_f)(cells[cell_idx])
-            for (entity_idx, cell_entity) in pairs(cell_entities)
-                if all(x -> x in verts, cell_entity)
-                    push!(_set, $(entity_t)((cell_idx, entity_idx)))
-                end
-            end
-        end
-        return _set
-      end
-    end
-end
-
-"""
-    filterfaces(grid::AbstractGrid, faces::Set{FaceIndex}, f::Function; all::Bool=true)
-
-Returns the faces in `faces` that satisfy `f` as a `Set{FaceIndex}`.
-`all=true` implies that `f(x)` must return `true` for all nodal coordinates `x` on the face
-if the face should be added to the set, otherwise it suffices that `f(x)` returns `true` for
-one node.
-
-```julia-repl
-julia> using Ferrite; using Ferrite: filterfaces
-
-julia> grid = generate_grid(Tetrahedron, (2,2,2));
-
-julia> topology = ExclusiveTopology(grid);
-
-julia> addboundaryfaceset!(grid, topology, "b", x -> true);
-
-julia> filterfaces(grid, grid.facesets["b"], x -> x[3] ≈ -1)
-Set{FaceIndex} with 8 elements:
-  FaceIndex((7, 1))
-  FaceIndex((3, 1))
-  FaceIndex((21, 1))
-  FaceIndex((13, 1))
-  FaceIndex((19, 1))
-  FaceIndex((15, 1))
-  FaceIndex((1, 1))
-  FaceIndex((9, 1))
-```
-"""
-function filterfaces end
-
-"""
-    filteredges(grid::AbstractGrid, edges::Set{EdgeIndex}, f::Function; all::Bool=true)
-
-Returns the edges in `edges` that satisfy `f` as a `Set{EdgeIndex}`.
-`all=true` implies that `f(x)` must return `true` for all nodal coordinates `x` on the face
-if the face should be added to the set, otherwise it suffices that `f(x)` returns `true` for
-one node.
-
-```julia-repl
-julia> using Ferrite; using Ferrite: filteredges
-
-julia> grid = generate_grid(Tetrahedron, (1,1,1));
-
-julia> topology = ExclusiveTopology(grid);
-
-julia> addboundaryedgeset!(grid, topology, "b", x -> true);
-
-julia> filteredges(grid, grid.edgesets["b"], x -> x[3] ≈ -1)
-Set{EdgeIndex} with 8 elements:
-  EdgeIndex((1, 2))
-  EdgeIndex((3, 2))
-  EdgeIndex((4, 3))
-  EdgeIndex((1, 3))
-  EdgeIndex((3, 3))
-  EdgeIndex((1, 1))
-  EdgeIndex((3, 1))
-  EdgeIndex((2, 3))
-```
-"""
-function filteredges end
-
-"""
-    filtervertices(grid::AbstractGrid, vertices::Set{VertexIndex}, f::Function; all::Bool=true)
-
-Returns the vertices in `vertices` that satisfy `f` as a `Set{VertexIndex}`.
-`all=true` implies that `f(x)` must return `true` for all nodal coordinates `x` on the face
-if the face should be added to the set, otherwise it suffices that `f(x)` returns `true` for
-one node.
-
-```julia-repl
-julia> using Ferrite; using Ferrite: filtervertices
-
-julia> grid = generate_grid(Tetrahedron,(1,1,1));
-
-julia> topology = ExclusiveTopology(grid);
-
-julia> addboundaryvertexset!(grid, topology, "b", x -> true);
-
-julia> filtervertices(grid, grid.vertexsets["b"], x -> x[3] ≈ -1)
-Set{VertexIndex} with 12 elements:
-  VertexIndex((2, 3))
-  VertexIndex((4, 3))
-  VertexIndex((4, 1))
-  VertexIndex((3, 3))
-  VertexIndex((3, 2))
-  VertexIndex((1, 1))
-  VertexIndex((2, 1))
-  VertexIndex((3, 1))
-  VertexIndex((1, 3))
-  VertexIndex((5, 1))
-  VertexIndex((1, 2))
-  VertexIndex((6, 1))
-```
-"""
-function filtervertices end
-
-for (func,            entity_f,  entity_t) in (
-    (:filtervertices, :vertices, :VertexIndex),
-    (:filteredges,    :edges,    :EdgeIndex),
-    (:filterfaces,    :faces,    :FaceIndex),
-)
-    @eval begin
-        function $(func)(grid::AbstractGrid, set::Set{$(entity_t)}, f::Function; all::Bool=true)
-            _set = Set{$(entity_t)}()
-            cells = getcells(grid)
-            for entity in set # entities can be edges/vertices in the face/edge
-                cell = cells[entity[1]]
-                cell_entities = $(entity_f)(cell)
-                pass = all
-                for node_idx in cell_entities[entity[2]] # using cell entities common with boundary face
-                    v = f(grid.nodes[node_idx].x)
-                    all ? (!v && (pass = false; break)) : (v && (pass = true; break))
-                end
-                pass && push!(_set, entity)
-            end
-            return _set
-        end
-    end
-end
-
-"""
-    addboundaryfaceset!(grid::AbstractGrid, topology::ExclusiveTopology, name::String, f::Function; all::Bool=true)
-
-Adds a boundary faceset to the grid with key `name`.
-A faceset maps a `String` key to a `Set` of tuples corresponding to `(global_cell_id,
-local_face_id)`. Facesets are used to initialize `Dirichlet` structs, that are needed to
-specify the boundary for the `ConstraintHandler`. `all=true` implies that `f(x)` must return
-`true` for all nodal coordinates `x` on the face if the face should be added to the set,
-otherwise it suffices that `f(x)` returns `true` for one node.
-
-```julia-repl
-julia> using Ferrite
-
-julia> grid = generate_grid(Tetrahedron, (1,1,1));
-
-julia> topology = ExclusiveTopology(grid);
-
-julia> addboundaryfaceset!(grid, topology, "b", x -> true);
-
-julia> grid.facesets["b"]
-Set{FaceIndex} with 12 elements:
-  FaceIndex((3, 1))
-  FaceIndex((4, 3))
-  FaceIndex((3, 3))
-  FaceIndex((4, 1))
-  FaceIndex((5, 1))
-  FaceIndex((2, 2))
-  FaceIndex((1, 4))
-  FaceIndex((2, 1))
-  FaceIndex((6, 1))
-  FaceIndex((6, 3))
-  FaceIndex((5, 3))
-  FaceIndex((1, 1))
-```
-"""
-function addboundaryfaceset! end
-
-"""
-    addboundaryedgeset!(grid::AbstractGrid, topology::ExclusiveTopology, name::String, f::Function; all::Bool=true)
-
-Adds a boundary edgeset to the grid with key `name`.
-An edgeset maps a `String` key to a `Set` of tuples corresponding to `(global_cell_id,
-local_edge_id)`. `all=true` implies that `f(x)` must return `true` for all nodal coordinates
-`x` on the face if the face should be added to the set, otherwise it suffices that `f(x)`
-returns `true` for one node.
-
-```julia-repl
-julia> using Ferrite
-
-julia> grid = generate_grid(Tetrahedron, (1,1,1));
-
-julia> topology = ExclusiveTopology(grid);
-
-julia> addboundaryedgeset!(grid, topology, "b", x -> true);
-
-julia> grid.edgesets["b"]
-Set{EdgeIndex} with 30 elements:
-  EdgeIndex((6, 6))
-  EdgeIndex((2, 1))
-  EdgeIndex((5, 3))
-  .
-  .
-  .
-  EdgeIndex((2, 5))
-  EdgeIndex((1, 4))
-```
-"""
-function addboundaryedgeset! end
-
-"""
-    addboundaryvertexset!(grid::AbstractGrid, topology::ExclusiveTopology, name::String, f::Function; all::Bool=true)
-
-Adds a boundary vertexset to the grid with key `name`.
-A vertexset maps a `String` key to a `Set` of tuples corresponding to `(global_cell_id,
-local_vertex_id)`. `all=true` implies that `f(x)` must return `true` for all nodal
-coordinates `x` on the face if the face should be added to the set, otherwise it suffices
-that `f(x)` returns `true` for one node.
-
-```julia-repl
-julia> using Ferrite
-
-julia> grid = generate_grid(Tetrahedron, (1,1,1));
-
-julia> topology = ExclusiveTopology(grid);
-
-julia> addboundaryvertexset!(grid, topology, "b", x -> true);
-
-julia> grid.vertexsets["b"]
-Set{VertexIndex} with 24 elements:
-  VertexIndex((2, 3))
-  VertexIndex((5, 2))
-  VertexIndex((4, 1))
-  .
-  .
-  .
-  VertexIndex((1, 4))
-  VertexIndex((4, 4))
-```
-"""
-function addboundaryvertexset! end
-
-for (func,                   entity_f,            entity_t,     filter_f,        instance_f,          destination) in (
-    (:addboundaryfaceset!,   :((_, x)->Set([x])), :FaceIndex,   :filterfaces,    :getfaceinstances,   :(grid.facesets)),
-    (:addboundaryedgeset!,   :getfaceedges,       :EdgeIndex,   :filteredges,    :getedgeinstances,   :(grid.edgesets)),
-    (:addboundaryvertexset!, :getfacevertices,    :VertexIndex, :filtervertices, :getvertexinstances, :(grid.vertexsets)),
-)
-    @eval begin
-        function $(func)(grid::AbstractGrid, topology::ExclusiveTopology, name::String, f::Function; all::Bool=true)
-            _check_setname($(destination), name)
-            _set = Set{$(entity_t)}()
-            for (face_idx, neighborhood) in pairs(topology.face_neighbor)
-                isempty(neighborhood) || continue # Skip any faces with neighbors (not on boundary)
-                entities =  $(entity_f)(grid, FaceIndex((face_idx[1], face_idx[2])))
-                for entity in $(filter_f)(grid, entities, f; all=all)
-                    union!(_set, $(instance_f)(grid, topology, entity))
-                end
-            end
-            _warn_emptyset(_set, name)
-            $(destination)[name] = _set
-            return grid
-        end
-    end
-end
-
 """
     addnodeset!(grid::AbstractGrid, name::String, nodeid::Union{Vector{Int},Set{Int}})
     addnodeset!(grid::AbstractGrid, name::String, f::Function)    
diff --git a/src/Grid/topology.jl b/src/Grid/topology.jl
new file mode 100644
index 0000000000..ed5f1e26ac
--- /dev/null
+++ b/src/Grid/topology.jl
@@ -0,0 +1,362 @@
+
+############
+# Topology #
+############
+
+"""
+    getneighborhood(topology, grid::AbstractGrid, cellidx::CellIndex, include_self=false)
+    getneighborhood(topology, grid::AbstractGrid, faceidx::FaceIndex, include_self=false)
+    getneighborhood(topology, grid::AbstractGrid, vertexidx::VertexIndex, include_self=false)
+    getneighborhood(topology, grid::AbstractGrid, edgeidx::EdgeIndex, include_self=false)
+
+Returns all connected entities of the same type as defined by the respective topology. If `include_self` is true,
+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` saves topological (connectivity/neighborhood) data of the grid. The constructor works with an `AbstractCell`
+vector for all cells that dispatch `vertices`, `faces` and in 3D `edges`.
+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.
+
+# 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.
+"""
+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
+    face_skeleton::Union{Vector{FaceIndex}, Nothing}
+    # TODO reintroduce the codimensional connectivity, e.g. 3D edge to 2D face
+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))")
+end
+
+function _num_shared_vertices(cell_a::C1, cell_b::C2) where {C1, C2}
+    num_shared_vertices = 0
+    for vertex ∈ vertices(cell_a)
+        for vertex_neighbor ∈ vertices(cell_b)
+            if vertex_neighbor == vertex
+                num_shared_vertices += 1
+                continue
+            end
+        end
+    end
+    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
+            getdim(cell_neighbor) == getdim(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 path
+            elseif num_shared_vertices == 2
+                if getdim(cell) == 2
+                    _add_single_face_neighbor!(face_table, cell, cell_id, cell_neighbor, cell_neighbor_id)
+                elseif getdim(cell) == 3
+                    _add_single_edge_neighbor!(edge_table, cell, cell_id, cell_neighbor, cell_neighbor_id)
+                else
+                    @error "Case not implemented."
+                end
+            # Shared surface
+            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)
+        dim = getdim(cells[1])
+
+        max_vertices = nvertices(cells[1])
+        dim > 1 && (max_faces = nfaces(cells[1]))
+        dim > 2 && (max_edges = nedges(cells[1]))
+    else
+        celltypes = Set(typeof.(cells))
+        for celltype in celltypes
+            celltypeidx = findfirst(x->typeof(x)==celltype,cells)
+            dim = getdim(cells[celltypeidx])
+
+            max_vertices = max(max_vertices,nvertices(cells[celltypeidx]))
+            dim > 1 && (max_faces = max(max_faces, nfaces(cells[celltypeidx])))
+            dim > 2 && (max_edges = max(max_edges, nedges(cells[celltypeidx])))
+        end
+    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)
+end
+
+function _add_single_face_neighbor!(face_table, cell::C1, cell_id, cell_neighbor::C2, cell_neighbor_id) where {C1, C2}
+    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))
+                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}
+    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))
+                return
+            end
+        end
+    end
+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)
+
+ExclusiveTopology(grid::AbstractGrid) = ExclusiveTopology(getcells(grid))
+
+function getneighborhood(top::ExclusiveTopology, grid::AbstractGrid, cellidx::CellIndex, include_self=false)
+    patch = getcells(top.cell_neighbor[cellidx.idx])
+    if include_self
+        return [patch; cellidx.idx]
+    else
+        return patch
+    end
+end
+
+function getneighborhood(top::ExclusiveTopology, grid::AbstractGrid, faceidx::FaceIndex, include_self=false)
+    if include_self
+        return [top.face_face_neighbor[faceidx[1],faceidx[2]].neighbor_info; faceidx]
+    else
+        return top.face_face_neighbor[faceidx[1],faceidx[2]].neighbor_info
+    end
+end
+
+function getneighborhood(top::ExclusiveTopology, grid::AbstractGrid, vertexidx::VertexIndex, include_self=false)
+    cellid, local_vertexid = vertexidx[1], vertexidx[2]
+    cell_vertices = vertices(getcells(grid,cellid))
+    global_vertexid = cell_vertices[local_vertexid]
+    if include_self
+        vertex_to_cell = top.vertex_to_cell[global_vertexid]
+        self_reference_local = Vector{VertexIndex}(undef,length(vertex_to_cell))
+        for (i,cellid) in enumerate(vertex_to_cell)
+            local_vertex = VertexIndex(cellid,findfirst(x->x==global_vertexid,vertices(getcells(grid,cellid)))::Int)
+            self_reference_local[i] = local_vertex
+        end
+        return [top.vertex_vertex_neighbor[global_vertexid].neighbor_info; self_reference_local]
+    else
+        return top.vertex_vertex_neighbor[global_vertexid].neighbor_info
+    end
+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))
+    nonlocal_edgeid = cell_edges[local_edgeidx]
+    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)))::Int
+        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])
+    else
+        return unique([top.edge_edge_neighbor[cellid, local_edgeidx].neighbor_info; self_reference_local])
+    end
+end
+
+"""
+    vertex_star_stencils(top::ExclusiveTopology, grid::Grid) -> Vector{Int, EntityNeighborhood{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 = getdim(cells[cell]) > 2 ? collect(edges(cells[cell])) : collect(faces(cells[cell])) #get lowest dimension boundary
+            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)
+        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}
+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
+end
+
+"""
+    _faceskeleton(topology::ExclusiveTopology, grid::Grid) -> Iterable{FaceIndex}
+Creates an iterateable face skeleton. The skeleton consists of `FaceIndex` that can be used to `reinit`
+`FaceValues`.
+"""
+function _faceskeleton(top::ExclusiveTopology, grid::Grid)
+    face_skeleton_global = Set{NTuple}()
+    face_skeleton_local = Vector{FaceIndex}()
+    fs_length = length(face_skeleton_global)
+    # TODO use topology to speed up :)
+    for (cellid,cell) ∈ enumerate(grid.cells)
+        for (local_face_id,face) ∈ enumerate(faces(cell))
+            push!(face_skeleton_global, sortface_fast(face))
+            fs_length_new = length(face_skeleton_global)
+            if fs_length != fs_length_new
+                push!(face_skeleton_local, FaceIndex(cellid,local_face_id))
+                fs_length = fs_length_new
+            end
+        end
+    end
+    return face_skeleton_local
+end
+
+"""
+    face_skeleton(top::ExclusiveTopology, grid::Grid) -> Vector{FaceIndex}
+Creates an iterateable face skeleton. The skeleton consists of `FaceIndex` that can be used to `reinit`
+`FaceValues`.
+"""
+function faceskeleton(top::ExclusiveTopology, grid::Grid)
+    if top.face_skeleton === nothing
+        top.face_skeleton = _faceskeleton(top, grid)
+    end
+    return top.face_skeleton
+end
diff --git a/src/Grid/utils.jl b/src/Grid/utils.jl
new file mode 100644
index 0000000000..61c76e04a7
--- /dev/null
+++ b/src/Grid/utils.jl
@@ -0,0 +1,342 @@
+"""
+getfaceinstances(grid::AbstractGrid, topology::ExclusiveTopology, face::FaceIndex)
+
+Returns all the faces as `Set{FaceIndex}` that share all their vertices with a given face
+represented as `FaceIndex`. The returned set includes the input face.
+
+```julia-repl
+julia> using Ferrite; using Ferrite: getfaceinstances
+
+julia> grid = generate_grid(Tetrahedron, (2,1,1));
+
+julia> topology = ExclusiveTopology(grid);
+
+julia> getfaceinstances(grid, topology, FaceIndex(4,2))
+Set{FaceIndex} with 2 elements:
+FaceIndex((6, 4))
+FaceIndex((4, 2))
+```
+"""
+function getfaceinstances end
+
+"""
+getedgeinstances(grid::AbstractGrid, topology::ExclusiveTopology, edge::EdgeIndex)
+
+Returns all the edges as `Set{EdgeIndex}` that share all their vertices with a given edge
+represented as `EdgeIndex`.
+The returned set includes the input edge.
+
+```julia-repl
+julia> using Ferrite; using Ferrite: getedgeinstances
+
+julia> grid = generate_grid(Tetrahedron, (2,1,1));
+
+julia> topology = ExclusiveTopology(grid);
+
+julia> getedgeinstances(grid, topology, EdgeIndex(4,2))
+Set{EdgeIndex} with 3 elements:
+EdgeIndex((4, 2))
+EdgeIndex((9, 6))
+EdgeIndex((7, 6))
+```
+"""
+function getedgeinstances end
+
+"""
+getvertexinstances(grid::AbstractGrid, topology::ExclusiveTopology, vertex::EdgeIndex)
+
+Returns all the vertices as `Set{::VertexIndex}` that use a given vertex represented as
+`VertexIndex` in all cells.
+The returned set includes the input vertex.
+
+```julia-repl
+julia> using Ferrite; using Ferrite: getvertexinstances
+
+julia> grid = generate_grid(Tetrahedron,(2,1,1));
+
+julia> topology = ExclusiveTopology(grid);
+
+julia> getvertexinstances(grid, topology, VertexIndex(4,2))
+Set{VertexIndex} with 8 elements:
+VertexIndex((7, 4))
+VertexIndex((10, 4))
+VertexIndex((12, 4))
+VertexIndex((6, 3))
+VertexIndex((4, 2))
+VertexIndex((9, 4))
+VertexIndex((11, 4))
+VertexIndex((8, 4))
+```
+"""
+function getvertexinstances end
+
+for (func,                entity_f,  entity_t) in (
+(:getvertexinstances, :vertices, :VertexIndex),
+(:getedgeinstances,   :edges,    :EdgeIndex),
+(:getfaceinstances,   :faces,    :FaceIndex),
+)
+@eval begin
+  function $(func)(grid::AbstractGrid, topology::ExclusiveTopology, entity::$(entity_t))
+    _set = Set{$(entity_t)}()
+    cells = getcells(grid)
+    cell = cells[entity[1]]
+    verts = $(entity_f)(cell)[entity[2]]
+    # Since we are looking for an entity that share *all* vertices, the first one can be
+    # used here to query potiential neighbor cells
+    for cell_idx in topology.vertex_to_cell[verts[1]] # Since all vertices should be shared, the first one can be used here
+        cell_entities = $(entity_f)(cells[cell_idx])
+        for (entity_idx, cell_entity) in pairs(cell_entities)
+            if all(x -> x in verts, cell_entity)
+                push!(_set, $(entity_t)((cell_idx, entity_idx)))
+            end
+        end
+    end
+    return _set
+  end
+end
+end
+
+"""
+filterfaces(grid::AbstractGrid, faces::Set{FaceIndex}, f::Function; all::Bool=true)
+
+Returns the faces in `faces` that satisfy `f` as a `Set{FaceIndex}`.
+`all=true` implies that `f(x)` must return `true` for all nodal coordinates `x` on the face
+if the face should be added to the set, otherwise it suffices that `f(x)` returns `true` for
+one node.
+
+```julia-repl
+julia> using Ferrite; using Ferrite: filterfaces
+
+julia> grid = generate_grid(Tetrahedron, (2,2,2));
+
+julia> topology = ExclusiveTopology(grid);
+
+julia> addboundaryfaceset!(grid, topology, "b", x -> true);
+
+julia> filterfaces(grid, grid.facesets["b"], x -> x[3] ≈ -1)
+Set{FaceIndex} with 8 elements:
+FaceIndex((7, 1))
+FaceIndex((3, 1))
+FaceIndex((21, 1))
+FaceIndex((13, 1))
+FaceIndex((19, 1))
+FaceIndex((15, 1))
+FaceIndex((1, 1))
+FaceIndex((9, 1))
+```
+"""
+function filterfaces end
+
+"""
+filteredges(grid::AbstractGrid, edges::Set{EdgeIndex}, f::Function; all::Bool=true)
+
+Returns the edges in `edges` that satisfy `f` as a `Set{EdgeIndex}`.
+`all=true` implies that `f(x)` must return `true` for all nodal coordinates `x` on the face
+if the face should be added to the set, otherwise it suffices that `f(x)` returns `true` for
+one node.
+
+```julia-repl
+julia> using Ferrite; using Ferrite: filteredges
+
+julia> grid = generate_grid(Tetrahedron, (1,1,1));
+
+julia> topology = ExclusiveTopology(grid);
+
+julia> addboundaryedgeset!(grid, topology, "b", x -> true);
+
+julia> filteredges(grid, grid.edgesets["b"], x -> x[3] ≈ -1)
+Set{EdgeIndex} with 8 elements:
+EdgeIndex((1, 2))
+EdgeIndex((3, 2))
+EdgeIndex((4, 3))
+EdgeIndex((1, 3))
+EdgeIndex((3, 3))
+EdgeIndex((1, 1))
+EdgeIndex((3, 1))
+EdgeIndex((2, 3))
+```
+"""
+function filteredges end
+
+"""
+filtervertices(grid::AbstractGrid, vertices::Set{VertexIndex}, f::Function; all::Bool=true)
+
+Returns the vertices in `vertices` that satisfy `f` as a `Set{VertexIndex}`.
+`all=true` implies that `f(x)` must return `true` for all nodal coordinates `x` on the face
+if the face should be added to the set, otherwise it suffices that `f(x)` returns `true` for
+one node.
+
+```julia-repl
+julia> using Ferrite; using Ferrite: filtervertices
+
+julia> grid = generate_grid(Tetrahedron,(1,1,1));
+
+julia> topology = ExclusiveTopology(grid);
+
+julia> addboundaryvertexset!(grid, topology, "b", x -> true);
+
+julia> filtervertices(grid, grid.vertexsets["b"], x -> x[3] ≈ -1)
+Set{VertexIndex} with 12 elements:
+VertexIndex((2, 3))
+VertexIndex((4, 3))
+VertexIndex((4, 1))
+VertexIndex((3, 3))
+VertexIndex((3, 2))
+VertexIndex((1, 1))
+VertexIndex((2, 1))
+VertexIndex((3, 1))
+VertexIndex((1, 3))
+VertexIndex((5, 1))
+VertexIndex((1, 2))
+VertexIndex((6, 1))
+```
+"""
+function filtervertices end
+
+for (func,            entity_f,  entity_t) in (
+(:filtervertices, :vertices, :VertexIndex),
+(:filteredges,    :edges,    :EdgeIndex),
+(:filterfaces,    :faces,    :FaceIndex),
+)
+@eval begin
+    function $(func)(grid::AbstractGrid, set::Set{$(entity_t)}, f::Function; all::Bool=true)
+        _set = Set{$(entity_t)}()
+        cells = getcells(grid)
+        for entity in set # entities can be edges/vertices in the face/edge
+            cell = cells[entity[1]]
+            cell_entities = $(entity_f)(cell)
+            pass = all
+            for node_idx in cell_entities[entity[2]] # using cell entities common with boundary face
+                v = f(grid.nodes[node_idx].x)
+                all ? (!v && (pass = false; break)) : (v && (pass = true; break))
+            end
+            pass && push!(_set, entity)
+        end
+        return _set
+    end
+end
+end
+
+"""
+addboundaryfaceset!(grid::AbstractGrid, topology::ExclusiveTopology, name::String, f::Function; all::Bool=true)
+
+Adds a boundary faceset to the grid with key `name`.
+A faceset maps a `String` key to a `Set` of tuples corresponding to `(global_cell_id,
+local_face_id)`. Facesets are used to initialize `Dirichlet` structs, that are needed to
+specify the boundary for the `ConstraintHandler`. `all=true` implies that `f(x)` must return
+`true` for all nodal coordinates `x` on the face if the face should be added to the set,
+otherwise it suffices that `f(x)` returns `true` for one node.
+
+```julia-repl
+julia> using Ferrite
+
+julia> grid = generate_grid(Tetrahedron, (1,1,1));
+
+julia> topology = ExclusiveTopology(grid);
+
+julia> addboundaryfaceset!(grid, topology, "b", x -> true);
+
+julia> grid.facesets["b"]
+Set{FaceIndex} with 12 elements:
+FaceIndex((3, 1))
+FaceIndex((4, 3))
+FaceIndex((3, 3))
+FaceIndex((4, 1))
+FaceIndex((5, 1))
+FaceIndex((2, 2))
+FaceIndex((1, 4))
+FaceIndex((2, 1))
+FaceIndex((6, 1))
+FaceIndex((6, 3))
+FaceIndex((5, 3))
+FaceIndex((1, 1))
+```
+"""
+function addboundaryfaceset! end
+
+"""
+addboundaryedgeset!(grid::AbstractGrid, topology::ExclusiveTopology, name::String, f::Function; all::Bool=true)
+
+Adds a boundary edgeset to the grid with key `name`.
+An edgeset maps a `String` key to a `Set` of tuples corresponding to `(global_cell_id,
+local_edge_id)`. `all=true` implies that `f(x)` must return `true` for all nodal coordinates
+`x` on the face if the face should be added to the set, otherwise it suffices that `f(x)`
+returns `true` for one node.
+
+```julia-repl
+julia> using Ferrite
+
+julia> grid = generate_grid(Tetrahedron, (1,1,1));
+
+julia> topology = ExclusiveTopology(grid);
+
+julia> addboundaryedgeset!(grid, topology, "b", x -> true);
+
+julia> grid.edgesets["b"]
+Set{EdgeIndex} with 30 elements:
+EdgeIndex((6, 6))
+EdgeIndex((2, 1))
+EdgeIndex((5, 3))
+.
+.
+.
+EdgeIndex((2, 5))
+EdgeIndex((1, 4))
+```
+"""
+function addboundaryedgeset! end
+
+"""
+addboundaryvertexset!(grid::AbstractGrid, topology::ExclusiveTopology, name::String, f::Function; all::Bool=true)
+
+Adds a boundary vertexset to the grid with key `name`.
+A vertexset maps a `String` key to a `Set` of tuples corresponding to `(global_cell_id,
+local_vertex_id)`. `all=true` implies that `f(x)` must return `true` for all nodal
+coordinates `x` on the face if the face should be added to the set, otherwise it suffices
+that `f(x)` returns `true` for one node.
+
+```julia-repl
+julia> using Ferrite
+
+julia> grid = generate_grid(Tetrahedron, (1,1,1));
+
+julia> topology = ExclusiveTopology(grid);
+
+julia> addboundaryvertexset!(grid, topology, "b", x -> true);
+
+julia> grid.vertexsets["b"]
+Set{VertexIndex} with 24 elements:
+VertexIndex((2, 3))
+VertexIndex((5, 2))
+VertexIndex((4, 1))
+.
+.
+.
+VertexIndex((1, 4))
+VertexIndex((4, 4))
+```
+"""
+function addboundaryvertexset! end
+
+for (func,                   entity_f,            entity_t,     filter_f,        instance_f,          destination) in (
+(:addboundaryfaceset!,   :((_, x)->Set([x])), :FaceIndex,   :filterfaces,    :getfaceinstances,   :(grid.facesets)),
+(:addboundaryedgeset!,   :getfaceedges,       :EdgeIndex,   :filteredges,    :getedgeinstances,   :(grid.edgesets)),
+(:addboundaryvertexset!, :getfacevertices,    :VertexIndex, :filtervertices, :getvertexinstances, :(grid.vertexsets)),
+)
+@eval begin
+    function $(func)(grid::AbstractGrid, topology::ExclusiveTopology, name::String, f::Function; all::Bool=true)
+        _check_setname($(destination), name)
+        _set = Set{$(entity_t)}()
+        for (face_idx, neighborhood) in pairs(topology.face_face_neighbor)
+            isempty(neighborhood) || continue # Skip any faces with neighbors (not on boundary)
+            entities =  $(entity_f)(grid, FaceIndex((face_idx[1], face_idx[2])))
+            for entity in $(filter_f)(grid, entities, f; all=all)
+                union!(_set, $(instance_f)(grid, topology, entity))
+            end
+        end
+        _warn_emptyset(_set, name)
+        $(destination)[name] = _set
+        return grid
+    end
+end
+end
diff --git a/src/exports.jl b/src/exports.jl
index 1ea17ce9d9..d5c4a9536a 100644
--- a/src/exports.jl
+++ b/src/exports.jl
@@ -67,6 +67,8 @@ export
     ExclusiveTopology,
     getneighborhood,
     faceskeleton,
+    vertex_star_stencils,
+    getstencil,
     getcells,
     getncells,
     getnodes,
diff --git a/test/test_grid_addboundaryset.jl b/test/test_grid_addboundaryset.jl
index 6c69fbb78f..494eac6af2 100644
--- a/test/test_grid_addboundaryset.jl
+++ b/test/test_grid_addboundaryset.jl
@@ -165,7 +165,7 @@
         topology = ExclusiveTopology(grid);
         addfaceset!(grid, "all", x->true)
         addvertexset!(grid, "all", x->true)
-        @test ∪([Ferrite.getfaceinstances(grid, topology,face) for face in topology.face_skeleton]...)  == grid.facesets["all"]
+        @test ∪([Ferrite.getfaceinstances(grid, topology,face) for face in Ferrite.faceskeleton(topology, grid)]...)  == grid.facesets["all"]
     end
     @testset "addboundaryset" for cell_type in [
         # Line, # topology construction error
diff --git a/test/test_grid_dofhandler_vtk.jl b/test/test_grid_dofhandler_vtk.jl
index 4eecde12ac..5d5e95420e 100644
--- a/test/test_grid_dofhandler_vtk.jl
+++ b/test/test_grid_dofhandler_vtk.jl
@@ -209,11 +209,12 @@ end
 #
     linegrid = generate_grid(Line,(3,))
     linetopo = ExclusiveTopology(linegrid)
-    @test linetopo.vertex_neighbor[1,2] == Ferrite.EntityNeighborhood(VertexIndex(2,1))
-    @test linetopo.vertex_neighbor[2,1] == Ferrite.EntityNeighborhood(VertexIndex(1,2))
-    @test linetopo.vertex_neighbor[2,2] == Ferrite.EntityNeighborhood(VertexIndex(3,1))
-    @test linetopo.vertex_neighbor[3,1] == Ferrite.EntityNeighborhood(VertexIndex(2,2))
-    @test length(linetopo.face_skeleton) == 4
+    @test linetopo.vertex_vertex_neighbor[1,2] == Ferrite.EntityNeighborhood(VertexIndex(2,1))
+    @test linetopo.vertex_vertex_neighbor[2,1] == Ferrite.EntityNeighborhood(VertexIndex(1,2))
+    @test linetopo.vertex_vertex_neighbor[2,2] == Ferrite.EntityNeighborhood(VertexIndex(3,1))
+    @test linetopo.vertex_vertex_neighbor[3,1] == Ferrite.EntityNeighborhood(VertexIndex(2,2))
+    linefaceskeleton = Ferrite.faceskeleton(linetopo, linegrid)
+    @test length(linefaceskeleton) == 4
 
 #                           (11)
 #                   (10)+-----+-----+(12)
@@ -227,33 +228,33 @@ end
     quadgrid = generate_grid(Quadrilateral,(2,3))
     topology = ExclusiveTopology(quadgrid)
     #test vertex neighbors maps cellid and local vertex id to neighbor id and neighbor local vertex id
-    @test topology.vertex_neighbor[1,3] == Ferrite.EntityNeighborhood(VertexIndex(4,1))
-    @test topology.vertex_neighbor[2,4] == Ferrite.EntityNeighborhood(VertexIndex(3,2))
-    @test topology.vertex_neighbor[3,3] == Ferrite.EntityNeighborhood(VertexIndex(6,1))
-    @test topology.vertex_neighbor[3,2] == Ferrite.EntityNeighborhood(VertexIndex(2,4))
-    @test topology.vertex_neighbor[4,1] == Ferrite.EntityNeighborhood(VertexIndex(1,3))
-    @test topology.vertex_neighbor[4,4] == Ferrite.EntityNeighborhood(VertexIndex(5,2))
-    @test topology.vertex_neighbor[5,2] == Ferrite.EntityNeighborhood(VertexIndex(4,4))
-    @test topology.vertex_neighbor[6,1] == Ferrite.EntityNeighborhood(VertexIndex(3,3))
+    @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[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 face neighbor maps cellid and local face id to neighbor id and neighbor local face id
-    @test topology.face_neighbor[1,2] == Ferrite.EntityNeighborhood(FaceIndex(2,4))
-    @test topology.face_neighbor[1,3] == Ferrite.EntityNeighborhood(FaceIndex(3,1))
-    @test topology.face_neighbor[2,3] == Ferrite.EntityNeighborhood(FaceIndex(4,1))
-    @test topology.face_neighbor[2,4] == Ferrite.EntityNeighborhood(FaceIndex(1,2))
-    @test topology.face_neighbor[3,1] == Ferrite.EntityNeighborhood(FaceIndex(1,3))
-    @test topology.face_neighbor[3,2] == Ferrite.EntityNeighborhood(FaceIndex(4,4))
-    @test topology.face_neighbor[3,3] == Ferrite.EntityNeighborhood(FaceIndex(5,1))
-    @test topology.face_neighbor[4,1] == Ferrite.EntityNeighborhood(FaceIndex(2,3))
-    @test topology.face_neighbor[4,3] == Ferrite.EntityNeighborhood(FaceIndex(6,1))
-    @test topology.face_neighbor[4,4] == Ferrite.EntityNeighborhood(FaceIndex(3,2))
-    @test topology.face_neighbor[5,1] == Ferrite.EntityNeighborhood(FaceIndex(3,3))
-    @test topology.face_neighbor[5,2] == Ferrite.EntityNeighborhood(FaceIndex(6,4))
-    @test topology.face_neighbor[5,3] == Ferrite.EntityNeighborhood(Ferrite.BoundaryIndex[])
-    @test topology.face_neighbor[5,4] == Ferrite.EntityNeighborhood(Ferrite.BoundaryIndex[])
-    @test topology.face_neighbor[6,1] == Ferrite.EntityNeighborhood(FaceIndex(4,3))
-    @test topology.face_neighbor[6,2] == Ferrite.EntityNeighborhood(Ferrite.BoundaryIndex[])
-    @test topology.face_neighbor[6,3] == Ferrite.EntityNeighborhood(Ferrite.BoundaryIndex[])
-    @test topology.face_neighbor[6,4] == Ferrite.EntityNeighborhood(FaceIndex(5,2))
+    @test topology.face_face_neighbor[1,2] == Ferrite.EntityNeighborhood(FaceIndex(2,4))
+    @test topology.face_face_neighbor[1,3] == Ferrite.EntityNeighborhood(FaceIndex(3,1))
+    @test topology.face_face_neighbor[2,3] == Ferrite.EntityNeighborhood(FaceIndex(4,1))
+    @test topology.face_face_neighbor[2,4] == Ferrite.EntityNeighborhood(FaceIndex(1,2))
+    @test topology.face_face_neighbor[3,1] == Ferrite.EntityNeighborhood(FaceIndex(1,3))
+    @test topology.face_face_neighbor[3,2] == Ferrite.EntityNeighborhood(FaceIndex(4,4))
+    @test topology.face_face_neighbor[3,3] == Ferrite.EntityNeighborhood(FaceIndex(5,1))
+    @test topology.face_face_neighbor[4,1] == Ferrite.EntityNeighborhood(FaceIndex(2,3))
+    @test topology.face_face_neighbor[4,3] == Ferrite.EntityNeighborhood(FaceIndex(6,1))
+    @test topology.face_face_neighbor[4,4] == Ferrite.EntityNeighborhood(FaceIndex(3,2))
+    @test topology.face_face_neighbor[5,1] == Ferrite.EntityNeighborhood(FaceIndex(3,3))
+    @test topology.face_face_neighbor[5,2] == Ferrite.EntityNeighborhood(FaceIndex(6,4))
+    @test topology.face_face_neighbor[5,3] == Ferrite.EntityNeighborhood(Ferrite.BoundaryIndex[])
+    @test topology.face_face_neighbor[5,4] == Ferrite.EntityNeighborhood(Ferrite.BoundaryIndex[])
+    @test topology.face_face_neighbor[6,1] == Ferrite.EntityNeighborhood(FaceIndex(4,3))
+    @test topology.face_face_neighbor[6,2] == Ferrite.EntityNeighborhood(Ferrite.BoundaryIndex[])
+    @test topology.face_face_neighbor[6,3] == Ferrite.EntityNeighborhood(Ferrite.BoundaryIndex[])
+    @test topology.face_face_neighbor[6,4] == Ferrite.EntityNeighborhood(FaceIndex(5,2))
 #                         (8)
 #                (7) +-----+-----+(9)
 #                    |  3  |  4  |
@@ -270,12 +271,12 @@ end
 #                        (11)
     hexgrid = generate_grid(Hexahedron,(2,2,1))
     topology = ExclusiveTopology(hexgrid)
-    @test topology.edge_neighbor[1,11] == Ferrite.EntityNeighborhood(EdgeIndex(4,9))
+    @test topology.edge_edge_neighbor[1,11] == Ferrite.EntityNeighborhood(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 topology.edge_neighbor[2,12] == Ferrite.EntityNeighborhood(EdgeIndex(3,10))
+    @test topology.edge_edge_neighbor[2,12] == Ferrite.EntityNeighborhood(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 topology.edge_neighbor[3,10] == Ferrite.EntityNeighborhood(EdgeIndex(2,12))
-    @test topology.edge_neighbor[4,9] == Ferrite.EntityNeighborhood(EdgeIndex(1,11))
+    @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 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))]
@@ -284,11 +285,12 @@ 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))]
+
+    # regression for https://github.com/Ferrite-FEM/Ferrite.jl/issues/518
     serendipitygrid = generate_grid(SerendipityQuadraticHexahedron,(2,2,1))
     stopology = ExclusiveTopology(serendipitygrid)
-    # regression for https://github.com/Ferrite-FEM/Ferrite.jl/issues/518
-    @test all(stopology.face_neighbor .== topology.face_neighbor)
-    @test all(stopology.vertex_neighbor .== topology.vertex_neighbor)
+    @test all(stopology.face_face_neighbor .== topology.face_face_neighbor)
+    @test all(stopology.vertex_vertex_neighbor .== topology.vertex_vertex_neighbor)
 
 #                   +-----+-----+
 #                   |\  6 |\  8 |
@@ -302,32 +304,32 @@ 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_neighbor[3,3] == Ferrite.EntityNeighborhood([VertexIndex(5,2),VertexIndex(6,1),VertexIndex(7,1)])
+    @test topology.vertex_vertex_neighbor[3,3] == Ferrite.EntityNeighborhood([VertexIndex(5,2),VertexIndex(6,1),VertexIndex(7,1)])
+
     quadtrigrid = generate_grid(QuadraticTriangle,(2,2))
     quadtopology = ExclusiveTopology(trigrid)
     # add more regression for https://github.com/Ferrite-FEM/Ferrite.jl/issues/518
-    @test all(quadtopology.face_neighbor .== topology.face_neighbor)
-    @test all(quadtopology.vertex_neighbor .== topology.vertex_neighbor)
+    @test all(quadtopology.face_face_neighbor .== topology.face_face_neighbor)
+    @test all(quadtopology.vertex_vertex_neighbor .== topology.vertex_vertex_neighbor)
 
 # test mixed grid
-    cells = [
-        Hexahedron((1, 2, 3, 4, 5, 6, 7, 8)),
-        Hexahedron((11, 13, 14, 12, 15, 16, 17, 18)),
-        Quadrilateral((2, 9, 10, 3)),
-        Quadrilateral((9, 11, 12, 10)),
-        ]
-    nodes = [Node(coord) for coord in zeros(Vec{2,Float64}, 18)]
-    grid = Grid(cells, nodes)
-    topology = ExclusiveTopology(grid)
-
-    @test topology.face_neighbor[3,4] == Ferrite.EntityNeighborhood(EdgeIndex(1,2))
-    @test topology.edge_neighbor[1,2] == Ferrite.EntityNeighborhood(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_neighbor[1,6] == topology.face_neighbor[1,1] == zero(Ferrite.EntityNeighborhood{FaceIndex})
-    @test topology.edge_neighbor[1,1] == topology.edge_neighbor[1,3] == zero(Ferrite.EntityNeighborhood{FaceIndex})
-    @test topology.face_neighbor[3,1] == topology.face_neighbor[3,3] == zero(Ferrite.EntityNeighborhood{FaceIndex})
-    @test topology.face_neighbor[4,1] == topology.face_neighbor[4,3] == zero(Ferrite.EntityNeighborhood{FaceIndex})
-
+    # cells = [
+    #     Hexahedron((1, 2, 3, 4, 5, 6, 7, 8)),
+    #     Hexahedron((11, 13, 14, 12, 15, 16, 17, 18)),
+    #     Quadrilateral((2, 9, 10, 3)),
+    #     Quadrilateral((9, 11, 12, 10)),
+    #     ]
+    # nodes = [Node(coord) for coord in zeros(Vec{2,Float64}, 18)]
+    # grid = Grid(cells, nodes)
+    # topology = ExclusiveTopology(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))
+    # # 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})
+    # @test topology.face_face_neighbor[3,1] == topology.face_face_neighbor[3,3] == zero(Ferrite.EntityNeighborhood{FaceIndex})
+    # @test topology.face_face_neighbor[4,1] == topology.face_face_neighbor[4,3] == zero(Ferrite.EntityNeighborhood{FaceIndex})
 
 #
 #                   +-----+-----+-----+
@@ -352,31 +354,32 @@ end
     @test issubset([7,4,5,6,9], patches[8])
     @test issubset([8,5,6], patches[9])
 
-    @test Set(Ferrite.getneighborhood(topology, quadgrid, VertexIndex(1,1))) == Set([VertexIndex(1,2), VertexIndex(1,4)])
-    @test Set(Ferrite.getneighborhood(topology, quadgrid, VertexIndex(2,1))) == Set([VertexIndex(1,1), VertexIndex(1,3), VertexIndex(2,2), VertexIndex(2,4)])
-    @test Set(Ferrite.getneighborhood(topology, 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)])
-    @test Set(Ferrite.getneighborhood(topology, quadgrid, VertexIndex(1,1),true)) == Set([VertexIndex(1,2), VertexIndex(1,4), VertexIndex(1,1)])
-    @test Set(Ferrite.getneighborhood(topology, quadgrid, VertexIndex(2,1),true)) == Set([VertexIndex(1,1), VertexIndex(1,3), VertexIndex(2,2), VertexIndex(2,4), VertexIndex(1,2), VertexIndex(2,1)])
-    @test Set(Ferrite.getneighborhood(topology, quadgrid, VertexIndex(5,4),true)) == 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.getneighborhood(topology, quadgrid, VertexIndex(1,1)))) == Set([2,5])
-    @test Set(Ferrite.toglobal(quadgrid, Ferrite.getneighborhood(topology, quadgrid, VertexIndex(2,1)))) == Set([1,6,3])
-    @test Set(Ferrite.toglobal(quadgrid, Ferrite.getneighborhood(topology, quadgrid, VertexIndex(5,4)))) == Set([6,9,11,14])
-
-    @test Set(topology.face_skeleton) == Set([FaceIndex(1,1),FaceIndex(1,2),FaceIndex(1,3),FaceIndex(1,4),
+# 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])
+
+    face_skeleton = Ferrite.faceskeleton(topology, quadgrid)
+    @test Set(face_skeleton) == Set([FaceIndex(1,1),FaceIndex(1,2),FaceIndex(1,3),FaceIndex(1,4),
                                           FaceIndex(2,1),FaceIndex(2,2),FaceIndex(2,3),
                                           FaceIndex(3,1),FaceIndex(3,2),FaceIndex(3,3),
                                           FaceIndex(4,2),FaceIndex(4,3),FaceIndex(4,4),
                                           FaceIndex(5,2),FaceIndex(5,3),FaceIndex(6,2),FaceIndex(6,3),
                                           FaceIndex(7,2),FaceIndex(7,3),FaceIndex(7,4),
                                           FaceIndex(8,2),FaceIndex(8,3),FaceIndex(9,2),FaceIndex(9,3)])
-    @test length(topology.face_skeleton) == 4*3 + 3*4
+    @test length(face_skeleton) == 4*3 + 3*4
 
     quadratic_quadgrid = generate_grid(QuadraticQuadrilateral,(3,3))
     quadgrid_topology = ExclusiveTopology(quadratic_quadgrid)
-    @test quadgrid_topology.face_skeleton == topology.face_skeleton
+    quadface_skeleton = Ferrite.faceskeleton(topology, quadgrid)
+    @test quadface_skeleton == face_skeleton
     # add more regression for https://github.com/Ferrite-FEM/Ferrite.jl/issues/518
-    @test all(quadgrid_topology.face_neighbor .== topology.face_neighbor)
-    @test all(quadgrid_topology.vertex_neighbor .== topology.vertex_neighbor)
+    @test all(quadgrid_topology.face_face_neighbor .== topology.face_face_neighbor)
+    @test all(quadgrid_topology.vertex_vertex_neighbor .== topology.vertex_vertex_neighbor)
     quadratic_patches = Vector{Int}[Ferrite.getneighborhood(quadgrid_topology, quadratic_quadgrid, CellIndex(i)) for i in 1:getncells(quadratic_quadgrid)]
     @test all(patches .== quadratic_patches)
 
@@ -394,7 +397,6 @@ end
         Ferrite.reinit!(fv, coords, faceid)
     end
     reinit!(fv::FaceValues, faceid::FaceIndex, grid) = reinit!(fv,faceid[1],faceid[2],grid) # wrapper for reinit!(fv,cellid,faceid,grid)
-    face_neighbors_ele5 = nonzeros(topology.face_neighbor[5,:])
     ip = Lagrange{RefQuadrilateral, 1}()^2
     qr_face = FaceQuadratureRule{RefQuadrilateral}(2)
     fv_ele = FaceValues(qr_face, ip)
@@ -409,7 +411,7 @@ end
             dΩ = getdetJdV(fv_ele, q_point)
             normal_5 = getnormal(fv_ele, q_point)
             u_5_n = function_value(fv_ele, q_point, u_ele5) ⋅ normal_5
-            for neighbor_entity in face_neighbors_ele5[ele_faceid].neighbor_info # only one entity can be changed to get rid of the for loop
+            for neighbor_entity in getneighborhood(topology, quadgrid, FaceIndex(5, ele_faceid)) # only one entity can be changed to get rid of the for loop
                 reinit!(fv_neighbor, neighbor_entity, quadgrid)
                 normal_neighbor = getnormal(fv_neighbor, q_point)
                 u_neighbor = function_value(fv_neighbor, q_point, u_neighbors) ⋅ normal_neighbor