diff --git a/ext/FerriteMetis.jl b/ext/FerriteMetis.jl
index d2a82935ad..4405659c44 100644
--- a/ext/FerriteMetis.jl
+++ b/ext/FerriteMetis.jl
@@ -52,7 +52,7 @@ function Ferrite.compute_renumber_permutation(
         else
             count(couplings[sdhi][i, j] for i in 1:n, j in 1:n if i != j)
         end
-        buffer_length += entries_per_cell * length(sdh.cellset)
+        buffer_length += entries_per_cell * length(getcellset(sdh))
     end
     I = Vector{idx_t}(undef, buffer_length)
     J = Vector{idx_t}(undef, buffer_length)
@@ -60,7 +60,7 @@ function Ferrite.compute_renumber_permutation(
 
     for (sdhi, sdh) in pairs(dh.subdofhandlers)
         coupling === nothing || (coupling_fh = couplings[sdhi])
-        for cc in CellIterator(dh, sdh.cellset)
+        for cc in CellIterator(dh, getcellset(sdh))
             dofs = celldofs(cc)
             for (j, dofj) in pairs(dofs), (i, dofi) in pairs(dofs)
                 dofi == dofj && continue # Metis doesn't want the diagonal
diff --git a/src/Dofs/ConstraintHandler.jl b/src/Dofs/ConstraintHandler.jl
index 8d9b9f1cec..d8a65a5d67 100644
--- a/src/Dofs/ConstraintHandler.jl
+++ b/src/Dofs/ConstraintHandler.jl
@@ -827,7 +827,7 @@ function add!(ch::ConstraintHandler, dbc::Dirichlet)
         dbc.field_name in sdh.field_names || continue
         # Compute the intersection between dbc.set and the cellset of this
         # SubDofHandler and skip if the set is empty
-        filtered_set = filter_dbc_set(get_grid(ch.dh), sdh.cellset, dbc.faces)
+        filtered_set = filter_dbc_set(get_grid(ch.dh), getcellset(sdh), dbc.faces)
         isempty(filtered_set) && continue
         # Fetch information about the field on this SubDofHandler
         field_idx = find_field(sdh, dbc.field_name)
@@ -855,7 +855,7 @@ function add!(ch::ConstraintHandler, dbc::Dirichlet)
         filtered_dbc = Dirichlet(dbc.field_name, filtered_set, dbc.f, components)
         _add!(
             ch, filtered_dbc, filtered_dbc.faces, interpolation, n_comp,
-            field_offset(sdh, field_idx), bcvalues, sdh.cellset,
+            field_offset(sdh, field_idx), bcvalues, getcellset(sdh),
         )
         dbc_added = true
     end
diff --git a/src/Dofs/DofHandler.jl b/src/Dofs/DofHandler.jl
index a653c8a563..a49da3d117 100644
--- a/src/Dofs/DofHandler.jl
+++ b/src/Dofs/DofHandler.jl
@@ -36,7 +36,7 @@ function SubDofHandler(dh::DH, cellset) where {DH <: AbstractDofHandler}
     end
     # Make sure this set is disjoint with all other existing
     for sdh in dh.subdofhandlers
-        if !isdisjoint(cellset, sdh.cellset)
+        if !isdisjoint(cellset, getcellset(sdh))
             error("cellset not disjoint with sets in existing SubDofHandlers")
         end
     end
@@ -46,7 +46,8 @@ function SubDofHandler(dh::DH, cellset) where {DH <: AbstractDofHandler}
     return sdh
 end
 
-@inline getcelltype(sdh::SubDofHandler) = getcelltype(get_grid(sdh.dh), first(sdh.cellset))
+@inline getcelltype(sdh::SubDofHandler) = getcelltype(get_grid(sdh.dh), first(getcellset(sdh)))
+@inline getcellset(sdh::SubDofHandler) = sdh.cellset
 
 function Base.show(io::IO, mime::MIME"text/plain", sdh::SubDofHandler)
     println(io, typeof(sdh))
@@ -147,7 +148,7 @@ function celldofs!(global_dofs::Vector{Int}, dh::DofHandler, i::Int)
     return global_dofs
 end
 function celldofs!(global_dofs::Vector{Int}, sdh::SubDofHandler, i::Int)
-    @assert i in sdh.cellset
+    @assert i in getcellset(sdh)
     return celldofs!(global_dofs, sdh.dh, i)
 end
 
@@ -224,7 +225,7 @@ function add!(sdh::SubDofHandler, name::Symbol, ip::Interpolation)
     end
     
     # Check that interpolation is compatible with cells it it added to
-    refshape_sdh = getrefshape(getcells(sdh.dh.grid, first(sdh.cellset)))
+    refshape_sdh = getrefshape(getcells(sdh.dh.grid, first(getcellset(sdh))))
     if refshape_sdh !== getrefshape(ip)
         error("The refshape of the interpolation $(getrefshape(ip)) is incompatible with the refshape $refshape_sdh of the cells.")
     end
@@ -254,7 +255,7 @@ function add!(dh::DofHandler, name::Symbol, ip::Interpolation)
     elseif length(dh.subdofhandlers) == 1
         # Add to existing SubDofHandler (if it covers all cells)
         sdh = dh.subdofhandlers[1]
-        if length(sdh.cellset) != getncells(get_grid(dh))
+        if length(getcellset(sdh)) != getncells(get_grid(dh))
             error("can not add field to DofHandler with a SubDofHandler for a subregion")
         end
     else # length(dh.subdofhandlers) > 1
@@ -395,7 +396,7 @@ function _close_subdofhandler!(dh::DofHandler{sdim}, sdh::SubDofHandler, sdh_ind
 
     # 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 BitSet(getcellset(sdh))
         @debug println("Creating dofs for cell #$ci")
 
         # TODO: _check_cellset_intersections can be removed in favor of this assertion
diff --git a/src/Dofs/DofRenumbering.jl b/src/Dofs/DofRenumbering.jl
index 704d89a181..a98653c4f8 100644
--- a/src/Dofs/DofRenumbering.jl
+++ b/src/Dofs/DofRenumbering.jl
@@ -196,7 +196,7 @@ function compute_renumber_permutation(dh::DofHandler, _, order::DofOrder.Compone
     for sdh in dh.subdofhandlers
         dof_ranges = [dof_range(sdh, f) for f in eachindex(sdh.field_names)]
         global_idxs = [findfirst(x -> x === f, dh.field_names) for f in sdh.field_names]
-        for cell in CellIterator(dh, sdh.cellset, flags)
+        for cell in CellIterator(dh, getcellset(sdh), flags)
             cdofs = celldofs(cell)
             for (local_idx, global_idx) in pairs(global_idxs)
                 rng = dof_ranges[local_idx]
diff --git a/src/Dofs/apply_analytical.jl b/src/Dofs/apply_analytical.jl
index 847fd58f8c..41c833c846 100644
--- a/src/Dofs/apply_analytical.jl
+++ b/src/Dofs/apply_analytical.jl
@@ -38,10 +38,10 @@ function apply_analytical!(
         ip_fun = getfieldinterpolation(sdh, field_idx)
         field_dim = getfielddim(sdh, field_idx)
         celldofinds = dof_range(sdh, fieldname)
-        set_intersection = if length(cellset) == length(sdh.cellset) == getncells(get_grid(dh))
+        set_intersection = if length(cellset) == length(getcellset(sdh)) == getncells(get_grid(dh))
             BitSet(1:getncells(get_grid(dh)))
         else
-            intersect(BitSet(sdh.cellset), BitSet(cellset))
+            intersect(BitSet(getcellset(sdh)), BitSet(cellset))
         end
         _apply_analytical!(a, dh, celldofinds, field_dim, ip_fun, ip_geo, f, set_intersection)
     end
diff --git a/src/Dofs/sparsity_pattern.jl b/src/Dofs/sparsity_pattern.jl
index 46ef67c3d3..23ba5863a3 100644
--- a/src/Dofs/sparsity_pattern.jl
+++ b/src/Dofs/sparsity_pattern.jl
@@ -179,7 +179,7 @@ function _create_sparsity_pattern(dh::AbstractDofHandler, ch#=::Union{Constraint
     # of entries eliminated by constraints.
     max_buffer_length = ndofs(dh) # diagonal elements
     for (sdh_idx, sdh) in pairs(dh.subdofhandlers)
-        set = sdh.cellset
+        set = getcellset(sdh)
         n = ndofs_per_cell(sdh)
         entries_per_cell = if coupling === nothing
             sym ? div(n * (n + 1), 2) : n^2
@@ -197,7 +197,7 @@ function _create_sparsity_pattern(dh::AbstractDofHandler, ch#=::Union{Constraint
     for (sdh_idx, sdh) in pairs(dh.subdofhandlers)
         coupling === nothing || (coupling_sdh = couplings[sdh_idx])
         # TODO: Remove BitSet construction when SubDofHandler ensures sorted collections
-        set = BitSet(sdh.cellset)
+        set = BitSet(getcellset(sdh))
         n = ndofs_per_cell(sdh)
         resize!(global_dofs, n)
         @inbounds for element_id in set
diff --git a/src/PointEval/PointEvalHandler.jl b/src/PointEval/PointEvalHandler.jl
index 62b2dfdb49..9d99b2bfac 100644
--- a/src/PointEval/PointEvalHandler.jl
+++ b/src/PointEval/PointEvalHandler.jl
@@ -217,7 +217,7 @@ function evaluate_at_points!(out_vals::Vector{T2},
         ip = func_interpolations[sdh_idx]
         if ip !== nothing
             dofrange = dof_range(sdh, fname)
-            cellset = sdh.cellset
+            cellset = getcellset(sdh)
             _evaluate_at_points!(out_vals, dof_vals, ph, dh, ip, cellset, dofrange)
         end
     end
diff --git a/src/iterators.jl b/src/iterators.jl
index 85611406cd..8e09e375e6 100644
--- a/src/iterators.jl
+++ b/src/iterators.jl
@@ -267,7 +267,7 @@ function CellIterator(gridordh::Union{Grid,DofHandler}, flags::UpdateFlags)
     return CellIterator(gridordh, nothing, flags)
 end
 function CellIterator(sdh::SubDofHandler, flags::UpdateFlags=UpdateFlags())
-    CellIterator(CellCache(sdh, flags), sdh.cellset)
+    CellIterator(CellCache(sdh, flags), getcellset(sdh))
 end
 
 @inline _getset(ci::CellIterator) = ci.set
diff --git a/test/test_apply_analytical.jl b/test/test_apply_analytical.jl
index 45140a4109..5bec7a0f5a 100644
--- a/test/test_apply_analytical.jl
+++ b/test/test_apply_analytical.jl
@@ -68,7 +68,7 @@
         dofs = Set{Int}()
         for sdh in dh.subdofhandlers
             if !isnothing(Ferrite._find_field(sdh, field_name))
-                _global_dof_range!(dofs, sdh, field_name, sdh.cellset)
+                _global_dof_range!(dofs, sdh, field_name, getcellset(sdh))
             end
         end
         return sort!(collect(Int, dofs))
diff --git a/test/test_mixeddofhandler.jl b/test/test_mixeddofhandler.jl
index bd02cebb98..2fc11b445d 100644
--- a/test/test_mixeddofhandler.jl
+++ b/test/test_mixeddofhandler.jl
@@ -341,7 +341,7 @@ function test_2_element_heat_eq()
         qr = QuadratureRule{RefQuadrilateral}(2)
         interp = sdh.field_interpolations[1]
         cellvalues = CellValues(qr, interp)
-        doassemble(sdh.cellset, cellvalues, assembler, dh)
+        doassemble(getcellset(sdh), cellvalues, assembler, dh)
     end
 
     update!(ch, 0.0);