diff --git a/Project.toml b/Project.toml
index d38c648f5d..69476ec25b 100644
--- a/Project.toml
+++ b/Project.toml
@@ -6,6 +6,7 @@ version = "0.3.14"
 EnumX = "4e289a0a-7415-4d19-859d-a7e5c4648b56"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce"
+OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
 Preferences = "21216c6a-2e73-6563-6e65-726566657250"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
@@ -26,6 +27,7 @@ BlockArrays = "0.16"
 EnumX = "1"
 Metis = "1.3"
 NearestNeighbors = "0.4"
+OrderedCollections = "1"
 Preferences = "1"
 Reexport = "1"
 Tensors = "1.14"
diff --git a/docs/Manifest.toml b/docs/Manifest.toml
index 2c9f66289c..ca66e397c7 100644
--- a/docs/Manifest.toml
+++ b/docs/Manifest.toml
@@ -429,7 +429,7 @@ uuid = "29a986be-02c6-4525-aec4-84b980013641"
 version = "2.0.0"
 
 [[deps.Ferrite]]
-deps = ["EnumX", "LinearAlgebra", "NearestNeighbors", "Preferences", "Reexport", "SparseArrays", "StaticArrays", "Tensors", "WriteVTK"]
+deps = ["EnumX", "LinearAlgebra", "NearestNeighbors", "OrderedCollections", "Preferences", "Reexport", "SparseArrays", "StaticArrays", "Tensors", "WriteVTK"]
 path = ".."
 uuid = "c061ca5d-56c9-439f-9c0e-210fe06d3992"
 version = "0.3.14"
@@ -450,9 +450,9 @@ version = "1.0.1"
 
 [[deps.FerriteMeshParser]]
 deps = ["Ferrite"]
-git-tree-sha1 = "8b948577bc4066e9c8693438fd511309c7383761"
+git-tree-sha1 = "47a1215bc76454082e265cb5193cffc9700bccb7"
 uuid = "0f8c756f-80dd-4a75-85c6-b0a5ab9d4620"
-version = "0.1.7"
+version = "0.1.8"
 
 [[deps.FileIO]]
 deps = ["Pkg", "Requires", "UUIDs"]
diff --git a/src/Dofs/ConstraintHandler.jl b/src/Dofs/ConstraintHandler.jl
index d398bfeae0..b5f13c7772 100644
--- a/src/Dofs/ConstraintHandler.jl
+++ b/src/Dofs/ConstraintHandler.jl
@@ -33,14 +33,14 @@ which applies the condition via [`apply!`](@ref) and/or [`apply_zero!`](@ref).
 """
 struct Dirichlet # <: Constraint
     f::Function # f(x) or f(x,t) -> value(s)
-    facets::Union{Set{Int},Set{FacetIndex},Set{FaceIndex},Set{EdgeIndex},Set{VertexIndex}}
+    facets::Union{OrderedSet{Int},OrderedSet{FacetIndex},OrderedSet{FaceIndex},OrderedSet{EdgeIndex},OrderedSet{VertexIndex}}
     field_name::Symbol
     components::Vector{Int} # components of the field
     local_facet_dofs::Vector{Int}
     local_facet_dofs_offset::Vector{Int}
 end
-function Dirichlet(field_name::Symbol, facets::Set, f::Function, components=nothing)
-    return Dirichlet(f, facets, field_name, __to_components(components), Int[], Int[])
+function Dirichlet(field_name::Symbol, facets::AbstractVecOrSet, f::Function, components=nothing)
+    return Dirichlet(f, convert_to_orderedset(facets), field_name, __to_components(components), Int[], Int[])
 end
 
 # components=nothing is default and means that all components should be constrained
@@ -293,7 +293,7 @@ function add_prescribed_dof!(ch::ConstraintHandler, constrained_dof::Int, inhomo
 end
 
 # Dirichlet on (facet|face|edge|vertex)set
-function _add!(ch::ConstraintHandler, dbc::Dirichlet, bcfacets::Set{Index}, interpolation::Interpolation, field_dim::Int, offset::Int, bcvalue::BCValues, _) where {Index<:BoundaryIndex}
+function _add!(ch::ConstraintHandler, dbc::Dirichlet, bcfacets::AbstractVecOrSet{Index}, interpolation::Interpolation, field_dim::Int, offset::Int, bcvalue::BCValues, _) where {Index<:BoundaryIndex}
     local_facet_dofs, local_facet_dofs_offset =
         _local_facet_dofs_for_bc(interpolation, field_dim, dbc.components, offset, dirichlet_boundarydof_indices(eltype(bcfacets)))
     copy!(dbc.local_facet_dofs, local_facet_dofs)
@@ -335,7 +335,7 @@ function _local_facet_dofs_for_bc(interpolation, field_dim, components, offset,
     return local_facet_dofs, local_facet_dofs_offset
 end
 
-function _add!(ch::ConstraintHandler, dbc::Dirichlet, bcnodes::Set{Int}, interpolation::Interpolation, field_dim::Int, offset::Int, bcvalue::BCValues, cellset::Set{Int}=Set{Int}(1:getncells(get_grid(ch.dh))))
+function _add!(ch::ConstraintHandler, dbc::Dirichlet, bcnodes::AbstractVecOrSet{Int}, interpolation::Interpolation, field_dim::Int, offset::Int, bcvalue::BCValues, cellset::AbstractVecOrSet{Int}=OrderedSet{Int}(1:getncells(get_grid(ch.dh))))
     grid = get_grid(ch.dh)
     if interpolation !== default_interpolation(getcelltype(grid, first(cellset)))
         @warn("adding constraint to nodeset is not recommended for sub/super-parametric approximations.")
@@ -427,7 +427,7 @@ function update!(ch::ConstraintHandler, time::Real=0.0)
 end
 
 # for facets, vertices, faces and edges
-function _update!(inhomogeneities::Vector{T}, f::Function, boundary_entities::Set{<:BoundaryIndex}, field::Symbol, local_facet_dofs::Vector{Int}, local_facet_dofs_offset::Vector{Int},
+function _update!(inhomogeneities::Vector{T}, f::Function, boundary_entities::AbstractVecOrSet{<:BoundaryIndex}, field::Symbol, local_facet_dofs::Vector{Int}, local_facet_dofs_offset::Vector{Int},
                   components::Vector{Int}, dh::AbstractDofHandler, boundaryvalues::BCValues,
                   dofmapping::Dict{Int,Int}, dofcoefficients::Vector{Union{Nothing,DofCoefficients{T}}}, time::Real) where {T}
 
@@ -464,7 +464,7 @@ function _update!(inhomogeneities::Vector{T}, f::Function, boundary_entities::Se
 end
 
 # for nodes
-function _update!(inhomogeneities::Vector{T}, f::Function, ::Set{Int}, field::Symbol, nodeidxs::Vector{Int}, globaldofs::Vector{Int},
+function _update!(inhomogeneities::Vector{T}, f::Function, ::AbstractVecOrSet{Int}, field::Symbol, nodeidxs::Vector{Int}, globaldofs::Vector{Int},
                   components::Vector{Int}, dh::AbstractDofHandler, facetvalues::BCValues,
                   dofmapping::Dict{Int,Int}, dofcoefficients::Vector{Union{Nothing,DofCoefficients{T}}}, time::Real) where T
     counter = 1
@@ -895,7 +895,7 @@ end
 
 function filter_dbc_set(grid::AbstractGrid, fhset::AbstractSet{Int}, dbcset::AbstractSet{Int})
     ret = empty(dbcset)
-    nodes_in_fhset = Set{Int}()
+    nodes_in_fhset = OrderedSet{Int}()
     for cc in CellIterator(grid, fhset, UpdateFlags(; nodes=true, coords=false))
         union!(nodes_in_fhset, cc.nodes)
     end
@@ -1096,7 +1096,7 @@ function _add!(ch::ConstraintHandler, pdbc::PeriodicDirichlet, interpolation::In
         points = construct_cornerish(min_x, max_x)
         tree = KDTree(Tx[get_node_coordinate(grid, i) for i in all_node_idxs_v])
         idxs, _ = NearestNeighbors.nn(tree, points)
-        corner_set = Set{Int}(all_node_idxs_v[i] for i in idxs)
+        corner_set = OrderedSet{Int}(all_node_idxs_v[i] for i in idxs)
 
         dbc = Dirichlet(pdbc.field_name, corner_set,
             pdbc.func === nothing ? (x, _) -> pdbc.components * eltype(x)(0) : pdbc.func,
@@ -1123,7 +1123,7 @@ function _add!(ch::ConstraintHandler, pdbc::PeriodicDirichlet, interpolation::In
     if pdbc.func !== nothing
         # Create another temp constraint handler if we need to compute inhomogeneities
         chtmp2 = ConstraintHandler(ch.dh)
-        all_facets = Set{FacetIndex}()
+        all_facets = OrderedSet{FaceIndex}()
         union!(all_facets, (x.mirror for x in facet_map))
         union!(all_facets, (x.image for x in facet_map))
         dbc_all = Dirichlet(pdbc.field_name, all_facets, pdbc.func, pdbc.components)
@@ -1292,7 +1292,7 @@ dictionary which maps each mirror facet to a image facet. The result can then be
 [`PeriodicDirichlet`](@ref).
 
 `mset` and `iset` can be given as a `String` (an existing facet set in the grid) or as a
-`Set{FacetIndex}` directly.
+`AbstractSet{FacetIndex}` directly.
 
 By default this function looks for a matching facet in the directions of the coordinate
 system. For other types of periodicities the `transform` function can be used. The
@@ -1304,12 +1304,12 @@ between a image-facet and mirror-facet, for them to be considered matched.
 
 See also: [`collect_periodic_facets!`](@ref), [`PeriodicDirichlet`](@ref).
 """
-function collect_periodic_facets(grid::Grid, mset::Union{Set{FacetIndex},String}, iset::Union{Set{FacetIndex},String}, transform::Union{Function,Nothing}=nothing; tol::Float64=1e-12)
+function collect_periodic_facets(grid::Grid, mset::Union{AbstractSet{FacetIndex},String}, iset::Union{AbstractSet{FacetIndex},String}, transform::Union{Function,Nothing}=nothing; tol::Float64=1e-12)
     return collect_periodic_facets!(PeriodicFacetPair[], grid, mset, iset, transform; tol)
 end
 
 """
-    collect_periodic_facets(grid::Grid, all_facets::Union{Set{FacetIndex},String,Nothing}=nothing; tol=1e-12)
+    collect_periodic_facets(grid::Grid, all_facets::Union{AbstractSet{FacetIndex},String,Nothing}=nothing; tol=1e-12)
 
 Split all facets in `all_facets` into image and mirror sets. For each matching pair, the facet
 located further along the vector `(1, 1, 1)` becomes the image facet.
@@ -1319,7 +1319,7 @@ have a neighbor) is used.
 
 See also: [`collect_periodic_facets!`](@ref), [`PeriodicDirichlet`](@ref).
 """
-function collect_periodic_facets(grid::Grid, all_facets::Union{Set{FacetIndex},String,Nothing}=nothing; tol::Float64=1e-12)
+function collect_periodic_facets(grid::Grid, all_facets::Union{AbstractSet{FacetIndex},String,Nothing}=nothing; tol::Float64=1e-12)
     return collect_periodic_facets!(PeriodicFacetPair[], grid, all_facets; tol)
 end
 
@@ -1329,7 +1329,7 @@ end
 
 Same as [`collect_periodic_facets`](@ref) but adds all matches to the existing `facet_map`.
 """
-function collect_periodic_facets!(facet_map::Vector{PeriodicFacetPair}, grid::Grid, mset::Union{Set{FacetIndex},String}, iset::Union{Set{FacetIndex},String}, transform::Union{Function,Nothing}=nothing; tol::Float64=1e-12)
+function collect_periodic_facets!(facet_map::Vector{PeriodicFacetPair}, grid::Grid, mset::Union{AbstractSet{FacetIndex},String}, iset::Union{AbstractSet{FacetIndex},String}, transform::Union{Function,Nothing}=nothing; tol::Float64=1e-12)
     mset = __to_facetset(grid, mset)
     iset = __to_facetset(grid, iset)
     if transform === nothing
@@ -1342,7 +1342,7 @@ function collect_periodic_facets!(facet_map::Vector{PeriodicFacetPair}, grid::Gr
     return facet_map
 end
 
-function collect_periodic_facets!(facet_map::Vector{PeriodicFacetPair}, grid::Grid, facetset::Union{Set{FacetIndex},String,Nothing}; tol::Float64=1e-12)
+function collect_periodic_facets!(facet_map::Vector{PeriodicFacetPair}, grid::Grid, facetset::Union{AbstractSet{FacetIndex},String,Nothing}; tol::Float64=1e-12)
     facetset = facetset === nothing ? __collect_boundary_facets(grid) : copy(__to_facetset(grid, facetset))
     if mod(length(facetset), 2) != 0
         error("uneven number of facets")
@@ -1350,7 +1350,7 @@ function collect_periodic_facets!(facet_map::Vector{PeriodicFacetPair}, grid::Gr
     return __collect_periodic_facets_bruteforce!(facet_map, grid, facetset, facetset, #=known_order=#false, tol)
 end
 
-__to_facetset(_, set::Set{FacetIndex}) = set
+__to_facetset(_, set::AbstractSet{FacetIndex}) = set
 __to_facetset(grid, set::String) = getfacetset(grid, set)
 function __collect_boundary_facets(grid::Grid)
     candidates = Dict{Tuple, FacetIndex}()
@@ -1364,7 +1364,7 @@ function __collect_boundary_facets(grid::Grid)
             end
         end
     end
-    return Set{FacetIndex}(values(candidates))
+    return OrderedSet{FacetIndex}(values(candidates))
 end
 
 function __collect_periodic_facets_tree!(facet_map::Vector{PeriodicFacetPair}, grid::Grid, mset::Vector{FacetIndex}, iset::Vector{FacetIndex}, transformation::F, tol::Float64) where F <: Function
@@ -1409,7 +1409,7 @@ function __collect_periodic_facets_tree!(facet_map::Vector{PeriodicFacetPair}, g
 end
 
 # This method empties mset and iset
-function __collect_periodic_facets_bruteforce!(facet_map::Vector{PeriodicFacetPair}, grid::Grid, mset::Set{FacetIndex}, iset::Set{FacetIndex}, known_order::Bool, tol::Float64)
+function __collect_periodic_facets_bruteforce!(facet_map::Vector{PeriodicFacetPair}, grid::Grid, mset::AbstractSet{FacetIndex}, iset::AbstractSet{FacetIndex}, known_order::Bool, tol::Float64)
     if length(mset) != length(iset)
         error("different facets in mirror and image")
     end
diff --git a/src/Dofs/DofHandler.jl b/src/Dofs/DofHandler.jl
index 12ec6f2e05..9c6eabf65c 100644
--- a/src/Dofs/DofHandler.jl
+++ b/src/Dofs/DofHandler.jl
@@ -15,7 +15,7 @@ get_grid(dh::AbstractDofHandler)
 struct SubDofHandler{DH} <: AbstractDofHandler
     # From constructor
     dh::DH
-    cellset::Set{Int}
+    cellset::OrderedSet{Int}
     # Populated in add!
     field_names::Vector{Symbol}
     field_interpolations::Vector{Interpolation}
@@ -26,7 +26,7 @@ struct SubDofHandler{DH} <: AbstractDofHandler
 end
 
 """
-    SubDofHandler(dh::AbstractDofHandler, cellset::Set{Int})
+    SubDofHandler(dh::AbstractDofHandler, cellset::AbstractVecOrSet{Int})
 
 Create an `sdh::SubDofHandler` from the parent `dh`, pertaining to the 
 cells in `cellset`. This allows you to add fields to parts of the domain, or using 
@@ -53,7 +53,7 @@ add!(sdh_quad, :u, ip_quad)
 close!(dh) # Finalize by closing the parent 
 ```
 """
-function SubDofHandler(dh::DH, cellset) where {DH <: AbstractDofHandler}
+function SubDofHandler(dh::DH, cellset::AbstractVecOrSet{Int}) where {DH <: AbstractDofHandler}
     # TODO: Should be an inner constructor.
     isclosed(dh) && error("DofHandler already closed")
     # Compute the celltype and make sure all elements have the same one
@@ -68,7 +68,7 @@ function SubDofHandler(dh::DH, cellset) where {DH <: AbstractDofHandler}
         end
     end
     # Construct and insert into the parent dh
-    sdh = SubDofHandler{typeof(dh)}(dh, cellset, Symbol[], Interpolation[], Int[], ScalarWrapper(-1))
+    sdh = SubDofHandler{typeof(dh)}(dh, convert_to_orderedset(cellset), Symbol[], Interpolation[], Int[], ScalarWrapper(-1))
     push!(dh.subdofhandlers, sdh)
     return sdh
 end
@@ -298,7 +298,7 @@ function add!(dh::DofHandler, name::Symbol, ip::Interpolation)
     @assert isconcretetype(celltype)
     if isempty(dh.subdofhandlers)
         # Create a new SubDofHandler for all cells
-        sdh = SubDofHandler(dh, Set(1:getncells(get_grid(dh))))
+        sdh = SubDofHandler(dh, OrderedSet(1:getncells(get_grid(dh))))
     elseif length(dh.subdofhandlers) == 1
         # Add to existing SubDofHandler (if it covers all cells)
         sdh = dh.subdofhandlers[1]
@@ -437,8 +437,7 @@ function _close_subdofhandler!(dh::DofHandler{sdim}, sdh::SubDofHandler, sdh_ind
     global_fidxs = Int[findfirst(gname -> gname === lname, dh.field_names) for lname in sdh.field_names]
 
     # loop over all the cells, and distribute dofs for all the fields
-    # TODO: Remove BitSet construction when SubDofHandler ensures sorted collections
-    for ci in BitSet(sdh.cellset)
+    for ci in sdh.cellset
         @debug println("Creating dofs for cell #$ci")
 
         # TODO: _check_cellset_intersections can be removed in favor of this assertion
diff --git a/src/Ferrite.jl b/src/Ferrite.jl
index 4ba6282a87..1c95b91ff3 100644
--- a/src/Ferrite.jl
+++ b/src/Ferrite.jl
@@ -13,6 +13,8 @@ using LinearAlgebra:
     pinv, tr
 using NearestNeighbors:
     NearestNeighbors, KDTree, knn
+using OrderedCollections:
+    OrderedSet
 using SparseArrays:
     SparseArrays, SparseMatrixCSC, nonzeros, nzrange, rowvals, sparse, spzeros
 using StaticArrays:
@@ -92,6 +94,9 @@ struct FacetIndex <: BoundaryIndex
     idx::Tuple{Int,Int} # cell and side
 end
 
+const AbstractVecOrSet{T} = Union{AbstractSet{T}, AbstractVector{T}}
+const IntegerCollection = AbstractVecOrSet{<:Integer}
+
 include("utils.jl")
 
 # Matrix/Vector utilities
diff --git a/src/Grid/grid.jl b/src/Grid/grid.jl
index 7b1ede7eef..a085191cf1 100644
--- a/src/Grid/grid.jl
+++ b/src/Grid/grid.jl
@@ -288,43 +288,51 @@ There are multiple helper structures to apply boundary conditions or define subd
 # Fields
 - `cells::Vector{C}`: stores all cells of the grid
 - `nodes::Vector{Node{dim,T}}`: stores the `dim` dimensional nodes of the grid
-- `cellsets::Dict{String,Set{Int}}`: maps a `String` key to a `Set` of cell ids
-- `nodesets::Dict{String,Set{Int}}`: maps a `String` key to a `Set` of global node ids
-- `facetsets::Dict{String,Set{FacetIndex}}`: maps a `String` to a `Set` of `Set{FacetIndex} (global_cell_id, local_facet_id)`
-- `vertexsets::Dict{String,Set{VertexIndex}}`: maps a `String` key to a `Set` of local vertex ids
+- `cellsets::Dict{String,OrderedSet{Int}}`: maps a `String` key to an `OrderedSet` of cell ids
+- `nodesets::Dict{String,OrderedSet{Int}}`: maps a `String` key to an `OrderedSet` of global node ids
+- `facetsets::Dict{String,OrderedSet{FacetIndex}}`: maps a `String` to an `OrderedSet` of `FacetIndex`
+- `vertexsets::Dict{String,OrderedSet{VertexIndex}}`: maps a `String` key to an `OrderedSet` of `VertexIndex`
 - `boundary_matrix::SparseMatrixCSC{Bool,Int}`: optional, only needed by `onboundary` to check if a cell is on the boundary, see, e.g. Helmholtz example
 """
 mutable struct Grid{dim,C<:AbstractCell,T<:Real} <: AbstractGrid{dim}
     cells::Vector{C}
     nodes::Vector{Node{dim,T}}
     # Sets
-    cellsets::Dict{String,Set{Int}}
-    nodesets::Dict{String,Set{Int}}
-    facetsets::Dict{String,Set{FacetIndex}}
-    vertexsets::Dict{String,Set{VertexIndex}}
+    cellsets::Dict{String,OrderedSet{Int}}
+    nodesets::Dict{String,OrderedSet{Int}}
+    facetsets::Dict{String,OrderedSet{FacetIndex}}
+    vertexsets::Dict{String,OrderedSet{VertexIndex}}
     # Boundary matrix (faces per cell × cell)
     boundary_matrix::SparseMatrixCSC{Bool,Int} # TODO: Deprecate!
 end
 
 function Grid(cells::Vector{C},
               nodes::Vector{Node{dim,T}};
-              cellsets::Dict{String,Set{Int}}=Dict{String,Set{Int}}(),
-              nodesets::Dict{String,Set{Int}}=Dict{String,Set{Int}}(),
-              facetsets::Dict{String,Set{FacetIndex}}=Dict{String,Set{FacetIndex}}(),
-              facesets = nothing,
-              vertexsets::Dict{String,Set{VertexIndex}}=Dict{String,Set{VertexIndex}}(),
+              cellsets::Dict{String, <:AbstractVecOrSet{Int}}=Dict{String,OrderedSet{Int}}(),
+              nodesets::Dict{String, <:AbstractVecOrSet{Int}}=Dict{String,OrderedSet{Int}}(),
+              facetsets::Dict{String, <:AbstractVecOrSet{FacetIndex}}=Dict{String,OrderedSet{FacetIndex}}(),
+              facesets=nothing, # deprecated
+              vertexsets::Dict{String, <:AbstractVecOrSet{VertexIndex}}=Dict{String,OrderedSet{VertexIndex}}(),
               boundary_matrix::SparseMatrixCSC{Bool,Int}=spzeros(Bool, 0, 0)) where {dim,C,T}
     if facesets !== nothing 
         if isempty(facetsets)
             @warn "facesets in Grid is deprecated, use facetsets instead" maxlog=1
             for (key, set) in facesets
-                facetsets[key] = Set(FacetIndex(cellnr, facenr) for (cellnr, facenr) in set)
+                facetsets[key] = OrderedSet(FacetIndex(cellnr, facenr) for (cellnr, facenr) in set)
             end
         else
             error("facesets are deprecated, use only facetsets")
         end
     end
-    return Grid(cells, nodes, cellsets, nodesets, facetsets, vertexsets, boundary_matrix)
+    return Grid(
+        cells,
+        nodes,
+        convert_to_orderedsets(cellsets),
+        convert_to_orderedsets(nodesets),
+        convert_to_orderedsets(facetsets),
+        convert_to_orderedsets(vertexsets),
+        boundary_matrix
+    )
 end
 
 ##########################
@@ -416,7 +424,7 @@ end
 """
     getcellset(grid::AbstractGrid, setname::String)
 
-Returns all cells as cellid in a `Set` of a given `setname`.
+Returns all cells as cellid in the set with name `setname`.
 """
 @inline getcellset(grid::AbstractGrid, setname::String) = grid.cellsets[setname]
 """
@@ -429,7 +437,7 @@ Returns all cellsets of the `grid`.
 """
     getnodeset(grid::AbstractGrid, setname::String)
 
-Returns all nodes as nodeid in a `Set` of a given `setname`.
+Returns all nodes as nodeid in the set with name `setname`.
 """
 @inline getnodeset(grid::AbstractGrid, setname::String) = grid.nodesets[setname]
 """
@@ -442,7 +450,7 @@ Returns all nodesets of the `grid`.
 """
     getfacetset(grid::AbstractGrid, setname::String)
 
-Returns all facets as `FacetIndex` in a `Set` of a given `setname`.
+Returns all faces as `FacetIndex` in the set with name `setname`.
 """
 @inline getfacetset(grid::AbstractGrid, setname::String) = grid.facetsets[setname]
 """
@@ -456,7 +464,7 @@ Returns all facet sets of the `grid`.
 """
     getvertexset(grid::AbstractGrid, setname::String)
 
-Returns all vertices as `VertexIndex` in a `Set` of a given `setname`.
+Returns all vertices as `VertexIndex` in the set with name `setname`.
 """
 @inline getvertexset(grid::AbstractGrid, setname::String) = grid.vertexsets[setname]
 """
@@ -535,7 +543,7 @@ function Base.show(io::IO, ::MIME"text/plain", grid::Grid)
     if isconcretetype(eltype(grid.cells))
         typestrs = [repr(eltype(grid.cells))]
     else
-        typestrs = sort!(repr.(Set(typeof(x) for x in grid.cells)))
+        typestrs = sort!(repr.(OrderedSet(typeof(x) for x in grid.cells)))
     end
     join(io, typestrs, '/')
     print(io, " cells and $(getnnodes(grid)) nodes")
diff --git a/src/Grid/utils.jl b/src/Grid/utils.jl
index ed02b8e893..4fc83a6df7 100644
--- a/src/Grid/utils.jl
+++ b/src/Grid/utils.jl
@@ -4,7 +4,7 @@ _check_setname(dict, name) = haskey(dict, name) && throw(ArgumentError("there al
 _warn_emptyset(set, name) = length(set) == 0 && @warn("no entities added to the set with name: $name")
 
 """
-    addcellset!(grid::AbstractGrid, name::String, cellid::Union{Set{Int}, Vector{Int}})
+    addcellset!(grid::AbstractGrid, name::String, cellid::AbstractVecOrSet{Int})
     addcellset!(grid::AbstractGrid, name::String, f::function; all::Bool=true)
 
 Adds a cellset to the grid with key `name`.
@@ -18,7 +18,7 @@ addcellset!(grid, "left", Set((1,3))) #add cells with id 1 and 3 to cellset left
 addcellset!(grid, "right", x -> norm(x[1]) < 2.0 ) #add cell to cellset right, if x[1] of each cell's node is smaller than 2.0
 ```
 """
-function addcellset!(grid::AbstractGrid, name::String, cellid::Union{Set{Int},Vector{Int}})
+function addcellset!(grid::AbstractGrid, name::String, cellid::AbstractVecOrSet{Int})
     _addset!(grid, name, cellid, getcellsets(grid))
 end
 
@@ -27,24 +27,24 @@ function addcellset!(grid::AbstractGrid, name::String, f::Function; all::Bool=tr
 end
 
 """
-    addnodeset!(grid::AbstractGrid, name::String, nodeid::Union{Vector{Int}, Set{Int}})
+    addnodeset!(grid::AbstractGrid, name::String, nodeid::AbstractVecOrSet{Int})
     addnodeset!(grid::AbstractGrid, name::String, f::Function)    
 
-Adds a `nodeset::Set{Int}` to the `grid`'s `nodesets` with key `name`. Has the same interface as `addcellset`. 
+Adds a `nodeset::OrderedSet{Int}` to the `grid`'s `nodesets` with key `name`. Has the same interface as `addcellset`. 
 However, instead of mapping a cell id to the `String` key, a set of node ids is returned.
 """
-addnodeset!(grid::AbstractGrid, name::String, nodeid::Union{Vector{Int}, Set{Int}}) = 
+addnodeset!(grid::AbstractGrid, name::String, nodeid::AbstractVecOrSet{Int}) = 
     _addset!(grid, name, nodeid, getnodesets(grid))
 
 addnodeset!(grid::AbstractGrid, name::String, f::Function) = 
     _addset!(grid, name, create_nodeset(grid, f), getnodesets(grid))
 
 """
-    addfacetset!(grid::AbstractGrid, name::String, faceid::Union{Set{FacetIndex},Vector{FacetIndex}})
+    addfacetset!(grid::AbstractGrid, name::String, faceid::AbstractVecOrSet{FacetIndex})
     addfacetset!(grid::AbstractGrid, name::String, f::Function; all::Bool=true) 
 
 Adds a facetset to the grid with key `name`.
-A facetset maps a `String` key to a `Set` of tuples corresponding to `(global_cell_id, local_facet_id)`.
+A facetset maps a `String` key to a `OrderedSet` of tuples corresponding to `(global_cell_id, local_facet_id)`.
 Facetsets can be used to initialize `Dirichlet` boundary conditions for the `ConstraintHandler`.
 `all=true` implies that `f(x)` must return `true` for all nodal coordinates `x` on the facet if the facet
 should be added to the set, otherwise it suffices that `f(x)` returns `true` for one node. 
@@ -54,18 +54,18 @@ addfacetset!(grid, "right", Set((FacetIndex(2,2), FacetIndex(4,2)))) #see grid m
 addfacetset!(grid, "clamped", x -> norm(x[1]) ≈ 0.0) #see incompressible elasticity example for reference
 ```
 """
-addfacetset!(grid::AbstractGrid, name::String, set::Union{Set{FacetIndex},Vector{FacetIndex}}) = 
+addfacetset!(grid::AbstractGrid, name::String, set::AbstractVecOrSet{FacetIndex}) = 
     _addset!(grid, name, set, getfacetsets(grid))
 
 addfacetset!(grid::AbstractGrid, name::String, f::Function; all::Bool=true) = 
     _addset!(grid, name, create_facetset(grid, f; all=all), getfacetsets(grid))
 
 """
-    addvertexset!(grid::AbstractGrid, name::String, faceid::Union{Set{FaceIndex},Vector{FaceIndex}})
+    addvertexset!(grid::AbstractGrid, name::String, faceid::AbstractVecOrSet{FaceIndex})
     addvertexset!(grid::AbstractGrid, name::String, f::Function) 
 
 Adds a 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)`.
+A vertexset maps a `String` key to a `OrderedSet` of tuples corresponding to `(global_cell_id, local_vertex_id)`.
 Vertexsets can be used to initialize `Dirichlet` boundary conditions for the `ConstraintHandler`.
 
 ```julia
@@ -73,15 +73,15 @@ addvertexset!(grid, "right", Set((VertexIndex(2,2), VertexIndex(4,2))))
 addvertexset!(grid, "clamped", x -> norm(x[1]) ≈ 0.0)
 ```
 """
-addvertexset!(grid::AbstractGrid, name::String, set::Union{Set{VertexIndex},Vector{VertexIndex}}) = 
+addvertexset!(grid::AbstractGrid, name::String, set::AbstractVecOrSet{VertexIndex}) = 
     _addset!(grid, name, set, getvertexsets(grid))
 
 addvertexset!(grid::AbstractGrid, name::String, f::Function) = 
     _addset!(grid, name, create_vertexset(grid, f; all=true), getvertexsets(grid))
 
-function _addset!(grid::AbstractGrid, name::String, _set, dict::Dict)
+function _addset!(grid::AbstractGrid, name::String, _set::AbstractVecOrSet, dict::Dict)
     _check_setname(dict, name)
-    set = Set(_set)
+    set = convert_to_orderedset(_set)
     _warn_emptyset(set, name)
     dict[name] = set
     grid
@@ -91,7 +91,7 @@ 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,
+A vertexset maps a `String` key to an `OrderedSet` 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.
@@ -105,7 +105,7 @@ end
 addboundaryfacetset!(grid::AbstractGrid, topology::ExclusiveTopology, name::String, f::Function; all::Bool=true)
 
 Adds a boundary facetset to the grid with key `name`.
-A facetset maps a `String` key to a `Set` of tuples corresponding to `(global_cell_id,
+A facetset maps a `String` key to a `OrderedSet` of tuples corresponding to `(global_cell_id,
 local_facet_id)`. Facetsets 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 facet if the facet should be added to the set,
@@ -117,7 +117,9 @@ function addboundaryfacetset!(grid::AbstractGrid, top::ExclusiveTopology, name::
 end
 
 function _create_set(f::Function, grid::AbstractGrid, ::Type{BI}; all=true) where {BI <: BoundaryIndex}
-    set = Set{BI}()
+    set = OrderedSet{BI}()
+    # Since we loop over the cells in order the resulting set will be sorted
+    # lexicographically based on the (cell_idx, entity_idx) tuple
     for (cell_idx, cell) in enumerate(getcells(grid))
         for (entity_idx, entity) in enumerate(boundaryfunction(BI)(cell))
             pass = all
@@ -133,7 +135,7 @@ end
 
 # Given a boundary index, for example EdgeIndex(2, 1), add this to `set`, as well as any other `EdgeIndex` in the grid 
 # pointing to the same edge (i.e. indices belong to neighboring cells)
-function push_entity_instances!(set::Set{BI}, grid::AbstractGrid, top::ExclusiveTopology, entity::BI) where {BI <: BoundaryIndex}
+function push_entity_instances!(set::OrderedSet{BI}, grid::AbstractGrid, top::ExclusiveTopology, entity::BI) where {BI <: BoundaryIndex}
     push!(set, entity) # Add the given entity
     cell = getcells(grid, entity[1])
     verts = boundaryfunction(BI)(cell)[entity[2]]
@@ -148,7 +150,7 @@ function push_entity_instances!(set::Set{BI}, grid::AbstractGrid, top::Exclusive
     return set
 end
 
-# Create a `Set{BI}` whose entities are a subset of facets which do not have neighbors, i.e. that are on the boundary. 
+# Create a `OrderedSet{BI}` whose entities are a subset of facets which do not have neighbors, i.e. that are on the boundary. 
 # Note that this may include entities from cells incident to the facet, e.g. 
 # ____  consider the case of a vertex boundary set, with f(x) only true on the right side. Then also the VertexIndex
 # |\ |  belong to the lower left cell, in its lower right corner, is on the boundary, so this should be added too.  
@@ -156,7 +158,7 @@ end
 function _create_boundaryset(f::Function, grid::AbstractGrid, top::ExclusiveTopology, ::Type{BI}; all = true) where {BI <: BoundaryIndex}
     # Function barrier as get_facet_facet_neighborhood is not always type stable
     function _makeset(ff_nh)
-        set = Set{BI}()
+        set = OrderedSet{BI}()
         for (ff_nh_idx, neighborhood) in pairs(ff_nh)
             # ff_nh_idx::CartesianIndex into Matrix{<:EntityNeighborhood}
             isempty(neighborhood) || continue # Skip any facets with neighbors (not on boundary)
@@ -177,11 +179,12 @@ function _create_boundaryset(f::Function, grid::AbstractGrid, top::ExclusiveTopo
         end
         return set
     end
-    return _makeset(get_facet_facet_neighborhood(top, grid))::Set{BI}
+    return _makeset(get_facet_facet_neighborhood(top, grid))::OrderedSet{BI}
 end
 
 function create_cellset(grid::AbstractGrid, f::Function; all::Bool=true)
-    cells = Set{Int}()
+    cells = OrderedSet{Int}()
+    # Since we loop over the cells in order the resulting set will be sorted
     for (i, cell) in enumerate(getcells(grid))
         pass = all
         for node_idx in get_node_ids(cell)
@@ -193,7 +196,7 @@ function create_cellset(grid::AbstractGrid, f::Function; all::Bool=true)
     return cells 
 end
 function create_nodeset(grid::AbstractGrid, f::Function)
-    nodes = Set{Int}()
+    nodes = OrderedSet{Int}()
     for (i, n) in pairs(getnodes(grid))
         f(get_node_coordinate(n)) && push!(nodes, i)
     end
diff --git a/src/L2_projection.jl b/src/L2_projection.jl
index d9bc37c939..98f6cd01df 100644
--- a/src/L2_projection.jl
+++ b/src/L2_projection.jl
@@ -6,7 +6,7 @@ struct L2Projector <: AbstractProjector
     geom_ip::Interpolation
     M_cholesky #::SuiteSparse.CHOLMOD.Factor{Float64}
     dh::DofHandler
-    set::Vector{Int}
+    set::OrderedSet{Int}
 end
 
 """
@@ -36,7 +36,7 @@ function L2Projector(
         func_ip::Interpolation,
         grid::AbstractGrid;
         qr_lhs::QuadratureRule = _mass_qr(func_ip),
-        set = 1:getncells(grid),
+        set = OrderedSet(1:getncells(grid)),
         geom_ip::Interpolation = default_interpolation(getcelltype(grid, first(set))),
     )
 
@@ -46,19 +46,20 @@ function L2Projector(
     end
 
     _check_same_celltype(grid, set)
+    _set = convert_to_orderedset(set)
 
     fe_values_mass = CellValues(qr_lhs, func_ip, geom_ip)
 
     # Create an internal scalar valued field. This is enough since the projection is done on a component basis, hence a scalar field.
     dh = DofHandler(grid)
-    sdh = SubDofHandler(dh, Set(set))
+    sdh = SubDofHandler(dh, _set)
     add!(sdh, :_, func_ip) # we need to create the field, but the interpolation is not used here
     close!(dh)
 
-    M = _assemble_L2_matrix(fe_values_mass, set, dh)  # the "mass" matrix
+    M = _assemble_L2_matrix(fe_values_mass, _set, dh)  # the "mass" matrix
     M_cholesky = cholesky(M)
 
-    return L2Projector(func_ip, geom_ip, M_cholesky, dh, collect(set))
+    return L2Projector(func_ip, geom_ip, M_cholesky, dh, _set)
 end
 
 # Quadrature sufficient for integrating a mass matrix
diff --git a/src/PointEvalHandler.jl b/src/PointEvalHandler.jl
index 7e27417d5a..7b7e2e14c3 100644
--- a/src/PointEvalHandler.jl
+++ b/src/PointEvalHandler.jl
@@ -234,7 +234,7 @@ function _evaluate_at_points!(
     ph::PointEvalHandler,
     dh::AbstractDofHandler,
     pv::PointValues,
-    cellset::Union{Nothing, Set{Int}},
+    cellset::Union{Nothing, AbstractSet{Int}},
     dofrange::AbstractRange{Int},
     ) where {T2,T}
 
diff --git a/src/iterators.jl b/src/iterators.jl
index ea2c834c3c..4d47d84133 100644
--- a/src/iterators.jl
+++ b/src/iterators.jl
@@ -225,9 +225,6 @@ getcoordinates(ic::InterfaceCache) = (getcoordinates(ic.a), getcoordinates(ic.b)
 ####################
 
 ## CellIterator ##
-
-const IntegerCollection = Union{Set{<:Integer}, AbstractVector{<:Integer}}
-
 """
     CellIterator(grid::Grid, cellset=1:getncells(grid))
     CellIterator(dh::AbstractDofHandler, cellset=1:getncells(dh))
@@ -289,7 +286,7 @@ FaceIterator(args...) = error("FaceIterator is deprecated, use FacetIterator ins
 
 # Leaving flags undocumented as for CellIterator
 """
-    FacetIterator(gridordh::Union{Grid,AbstractDofHandler}, faceset::Set{FacetIndex})
+    FacetIterator(gridordh::Union{Grid,AbstractDofHandler}, faceset::AbstractVecOrSet{FaceIndex})
 
 Create a `FacetIterator` to conveniently iterate over the faces in `faceset`. The elements of
 the iterator are [`FacetCache`](@ref)s which are properly `reinit!`ialized. See
@@ -311,7 +308,7 @@ end
 """
 struct FacetIterator{FC<:FacetCache}
     fc::FC
-    set::Set{FacetIndex}
+    set::OrderedSet{FacetIndex}
 end
 
 function FacetIterator(gridordh::Union{Grid,AbstractDofHandler},
@@ -410,7 +407,7 @@ function _check_same_celltype(grid::AbstractGrid, cellset::IntegerCollection)
     end
 end
 
-function _check_same_celltype(grid::AbstractGrid, facetset::Set{<:BoundaryIndex})
+function _check_same_celltype(grid::AbstractGrid, facetset::AbstractVecOrSet{FaceIndex})
     isconcretetype(getcelltype(grid)) && return nothing # Short circuit check
     celltype = getcelltype(grid, first(facetset)[1])
     if !all(getcelltype(grid, facet[1]) == celltype for facet in facetset)
diff --git a/src/utils.jl b/src/utils.jl
index 09e57b0277..c77469b357 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -42,3 +42,11 @@ end
 @inline Base.getindex(s::ScalarWrapper) = s.x
 @inline Base.setindex!(s::ScalarWrapper, v) = s.x = v
 Base.copy(s::ScalarWrapper{T}) where {T} = ScalarWrapper{T}(copy(s.x))
+
+convert_to_orderedset(set::AbstractVector{T}) where T = OrderedSet{T}(set)
+convert_to_orderedset(set::AbstractSet{T}) where T = convert(OrderedSet{T}, set)
+
+function convert_to_orderedsets(namedsets::Dict{String, <: AbstractVecOrSet{T}}) where T 
+    return Dict{String,OrderedSet{T}}(k => convert_to_orderedset(v) for (k,v) in namedsets)
+end
+convert_to_orderedsets(namedsets::Dict{String, <: OrderedSet}) = namedsets
diff --git a/test/runtests.jl b/test/runtests.jl
index f864f31b9a..cfbfdd43fe 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -7,6 +7,7 @@ import SHA
 using Random
 using LinearAlgebra
 using SparseArrays
+using OrderedCollections
 
 const HAS_EXTENSIONS = isdefined(Base, :get_extension)
 
diff --git a/test/test_mixeddofhandler.jl b/test/test_mixeddofhandler.jl
index 8e2d548e52..6324dfec4b 100644
--- a/test/test_mixeddofhandler.jl
+++ b/test/test_mixeddofhandler.jl
@@ -32,7 +32,7 @@ function test_1d_bar_beam()
     sdh1 = SubDofHandler(dh, Set(3))
     add!(sdh1, :u, ip^2)
     add!(sdh1, :θ, ip)
-    sdh2 = SubDofHandler(dh, Set((1,2)))
+    sdh2 = SubDofHandler(dh, OrderedSet((1,2)))
     add!(sdh2, :u, ip^2)
     close!(dh)
     @test ndofs(dh) == 8
@@ -386,7 +386,7 @@ function test_element_order()
 
     dh = DofHandler(grid);
     # Note the jump in cell numbers
-    sdh_tri = SubDofHandler(dh, Set((1,3)))
+    sdh_tri = SubDofHandler(dh, OrderedSet((1,3)))
     add!(sdh_tri, :u, Lagrange{RefTriangle,1}()^2)
     sdh_quad = SubDofHandler(dh, Set(2))
     add!(sdh_quad, :u, Lagrange{RefQuadrilateral,1}()^2)
@@ -447,8 +447,8 @@ function test_evaluate_at_grid_nodes()
     Triangle((3,4,6)),
     Triangle((3,6,5))]
     mesh = Grid(cells, nodes)
-    addcellset!(mesh, "quads", Set{Int}((1,)))
-    addcellset!(mesh, "tris", Set{Int}((2, 3)))
+    addcellset!(mesh, "quads", Set((1,)))
+    addcellset!(mesh, "tris", OrderedSet((2, 3)))
 
     ip_quad = Lagrange{RefQuadrilateral,1}()
     ip_tri = Lagrange{RefTriangle,1}()
@@ -554,8 +554,8 @@ function test_separate_fields_on_separate_domains()
             Triangle((3,4,5)),
             Triangle((4,6,5))]
     mesh = Grid(cells, nodes)
-    addcellset!(mesh, "quads", Set{Int}((1,)))
-    addcellset!(mesh, "tris", Set{Int}((2, 3)))
+    addcellset!(mesh, "quads", Set((1,)))
+    addcellset!(mesh, "tris", OrderedSet((2, 3)))
 
     ip_tri = Lagrange{RefTriangle,1}()
     ip_quad = Lagrange{RefQuadrilateral,1}()
@@ -578,8 +578,8 @@ end
 
 function test_unique_cellsets()
     grid = generate_grid(Quadrilateral, (2, 1))
-    set_u = Set(1:2)
-    set_v = Set(1:1)
+    set_u = OrderedSet(1:2)
+    set_v = OrderedSet(1:1)
 
     ip = Lagrange{RefQuadrilateral,1}()