diff --git a/docs/make.jl b/docs/make.jl
index 8f24bcf3b2..9410423598 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -14,7 +14,8 @@ GENERATEDEXAMPLES = [joinpath("examples", f) for f in (
     "threaded_assembly.md",
     "plasticity.md",
     "transient_heat_equation.md",
-    "landau.md"
+    "landau.md",
+    "linear_shell.md"
     )]
 
 # Build documentation.
diff --git a/docs/src/literate/linear_shell.jl b/docs/src/literate/linear_shell.jl
new file mode 100644
index 0000000000..e7d6ee0a99
--- /dev/null
+++ b/docs/src/literate/linear_shell.jl
@@ -0,0 +1,314 @@
+# # Linear shell
+#
+# ![](linear_shell.png)
+#-
+# ## Introduction
+#
+# In this example we show how shell elements can be analyzed in Ferrite.jl. The shell implemented here comes from the book 
+# "The finite elment method - Linear static and dynamic finite element analysis" by Hughes (1987), and a brief description of it is 
+# given at the end of this tutorial.  The first part of the tutorial explains how to set up the problem.
+
+# ## Setting up the problem
+using Ferrite
+using ForwardDiff
+
+function main() #wrap everything in a function...
+
+# First we generate a flat rectangular mesh. There is currently no built-in function for generating
+# shell meshes in Ferrite, so we have to create our own simple mesh generator (see the  
+# function `generate_shell_grid` further down in this file).
+#+
+nels = (10,10)
+size = (10.0, 10.0)
+grid = generate_shell_grid(nels, size)
+
+# Here we define the bi-linear interpolation used for the geometrical description of the shell.
+# We also create two quadrature rules for the in-plane and out-of-plane directions. Note that we use 
+# under integration for the inplane integration, to avoid shear locking. 
+#+
+ip = Lagrange{2,RefCube,1}()
+qr_inplane = QuadratureRule{2,RefCube}(1)
+qr_ooplane = QuadratureRule{1,RefCube}(2)
+cv = CellScalarValues(qr_inplane, ip)
+
+# Next we distribute displacement dofs,`:u = (x,y,z)` and rotational dofs, `:θ = (θ₁,  θ₂)`.
+#+
+dh = DofHandler(grid)
+push!(dh, :u, 3, ip)
+push!(dh, :θ, 2, ip)
+close!(dh)
+
+# In order to apply our boundary conditions, we first need to create some edge- and vertex-sets. This 
+# is done with `addedgeset!` and `addvertexset!` (similar to `addfaceset!`)
+#+
+addedgeset!(grid, "left",  (x) -> x[1] ≈ 0.0)
+addedgeset!(grid, "right", (x) -> x[1] ≈ size[1])
+addvertexset!(grid, "corner", (x) -> x[1] ≈ 0.0 && x[2] ≈ 0.0 && x[3] ≈ 0.0)
+
+# Here we define the boundary conditions. On the left edge, we lock the displacements in the x- and z- directions, and all the rotations.
+#+
+ch = ConstraintHandler(dh)
+add!(ch,  Dirichlet(:u, getedgeset(grid, "left"), (x, t) -> (0.0, 0.0), [1,3])  )
+add!(ch,  Dirichlet(:θ, getedgeset(grid, "left"), (x, t) -> (0.0, 0.0), [1,2])  )
+
+# On the right edge, we also lock the displacements in the x- and z- directions, but apply a precribed roation.
+#+
+add!(ch,  Dirichlet(:u, getedgeset(grid, "right"), (x, t) -> (0.0, 0.0), [1,3])  )
+add!(ch,  Dirichlet(:θ, getedgeset(grid, "right"), (x, t) -> (0.0, pi/10), [1,2])  )
+
+# In order to not get rigid body motion, we lock the y-displacement in one fo the corners.
+#+
+add!(ch,  Dirichlet(:θ, getvertexset(grid, "corner"), (x, t) -> (0.0), [2])  )
+
+close!(ch)
+update!(ch, 0.0)
+
+# Next we define relevant data for the shell, such as shear correction factor and stiffness matrix for the material. 
+# In this linear shell, plane stress is assumed, ie $\\sigma_{zz} = 0 $. Therefor, the stiffness matrix is 5x5 (opposed to the normal 6x6).
+#+
+κ = 5/6 # Shear correction factor
+E = 210.0
+ν = 0.3
+a = (1-ν)/2
+C = E/(1-ν^2) * [1 ν 0   0   0;
+                ν 1 0   0   0;
+                0 0 a*κ 0   0;
+                0 0 0   a*κ 0;
+                0 0 0   0   a*κ]
+
+
+data = (thickness = 1.0, C = C); #Named tuple
+
+# We now assemble the problem in standard finite element fashion
+#+
+nnodes = getnbasefunctions(ip)
+ndofs_shell = ndofs_per_cell(dh)
+
+K = create_sparsity_pattern(dh)
+f = zeros(Float64, ndofs(dh))
+
+ke = zeros(ndofs_shell, ndofs_shell)
+fe = zeros(ndofs_shell)
+
+celldofs = zeros(Int, ndofs_shell)
+cellcoords = zeros(Vec{3,Float64}, nnodes)
+
+assembler = start_assemble(K, f)
+for cellid in 1:getncells(grid)
+    fill!(ke, 0.0)
+
+    celldofs!(celldofs, dh, cellid)
+    getcoordinates!(cellcoords, grid, cellid)
+
+    #Call the element routine
+    integrate_shell!(ke, cv, qr_ooplane, cellcoords, data)
+
+    assemble!(assembler, celldofs, fe, ke)
+end
+
+# Apply BC and solve.
+#+
+apply!(K, f, ch)
+a = K\f
+
+# Output results.
+#+
+vtk_grid("linear_shell", dh) do vtk
+    vtk_point_data(vtk, dh, a)
+end
+
+end; #end main functions
+
+# Below is the function that creates the shell mesh. It simply generates a 2d-quadrature mesh, and appends
+# a third coordinate (z-direction) to the node-positions. 
+function generate_shell_grid(nels, size)
+    _grid = generate_grid(Quadrilateral, nels, Vec((0.0,0.0)), Vec(size))
+    nodes = [(n.x[1], n.x[2], 0.0) |> Vec{3} |> Node  for n in _grid.nodes]
+    cells = [Quadrilateral3D(cell.nodes) for cell in _grid.cells]
+
+    grid = Grid(cells, nodes)
+
+    return grid
+end;
+
+# ## The shell element
+#
+# The shell presented here comes from the book "The finite elment method - Linear static and dynamic finite element analysis" by Hughes (1987).
+# The shell is a so called degenerate shell element, meaning it is based on a continuum element.
+# A brief describtion of the shell is given here.
+
+#md # !!! note
+#md #     This element might experience various locking phenomenas, and should only be seen as a proof of concept.
+
+# ##### Fiber coordinate system 
+# The element uses two coordinate systems. The first coordianate system, called the fiber system, is created for each
+# element node, and is used as a reference frame for the rotations. The function below implements an algorthim that return the
+# fiber directions, $\boldsymbol{e}^{f}_{a1}$, $\boldsymbol{e}^{f}_{a2}$ and $\boldsymbol{e}^{f}_{a3}$, at each node $a$.
+function fiber_coordsys(Ps::Vector{Vec{3,Float64}})
+
+    ef1 = Vec{3,Float64}[]
+    ef2 = Vec{3,Float64}[]
+    ef3 = Vec{3,Float64}[]
+    for P in Ps
+        a = abs.(P)
+        j = 1
+        if a[1] > a[3]; a[3] = a[1]; j = 2; end
+        if a[2] > a[3]; j = 3; end
+       
+        e3 = P
+        e2 = Tensors.cross(P, basevec(Vec{3}, j))
+        e2 /= norm(e2)
+        e1 = Tensors.cross(e2, P)
+
+        push!(ef1, e1)
+        push!(ef2, e2)
+        push!(ef3, e3)
+    end
+    return ef1, ef2, ef3
+
+end;
+
+# ##### Lamina coordinate system 
+# The second coordinate system is the so called Lamina Coordinate system. It is
+# created for each integration point, and is defined to be tangent to the
+# mid-surface. It is in this system that we enforce that plane stress assumption,
+# i.e. $\sigma_{zz} = 0$. The function below returns the rotation matrix, $\boldsymbol{q}$, for this coordinate system.
+function lamina_coordsys(dNdξ, ζ, x, p, h)
+
+    e1 = zero(Vec{3})
+    e2 = zero(Vec{3})
+
+    for i in 1:length(dNdξ)
+        e1 += dNdξ[i][1] * x[i] + 0.5*h*ζ * dNdξ[i][1] * p[i]
+        e2 += dNdξ[i][2] * x[i] + 0.5*h*ζ * dNdξ[i][1] * p[i]
+    end
+
+    e1 /= norm(e1)
+    e2 /= norm(e2)
+
+    ez = Tensors.cross(e1,e2)
+    ez /= norm(ez)
+
+    a = 0.5*(e1 + e2)
+    a /= norm(a)
+
+    b = Tensors.cross(ez,a)
+    b /= norm(b)
+
+    ex = sqrt(2)/2 * (a - b)
+    ey = sqrt(2)/2 * (a + b)
+
+    return Tensor{2,3}(hcat(ex,ey,ez))
+end;
+
+
+# ##### Geometrical description
+# A material point in the shell is defined as
+# ```math
+# \boldsymbol x(\xi, \eta, \zeta) = \sum_{a=1}^{N_{\text{nodes}}} N_a(\xi, \eta) \boldsymbol{\bar{x}}_{a} + ζ \frac{h}{2} \boldsymbol{\bar{p}_a}
+# ``` 
+# where $\boldsymbol{\bar{x}}_{a}$ are nodal positions on the mid-surface, and $\boldsymbol{\bar{p}_a}$ is an vector that defines the fiber direction
+# on the reference surface. $N_a$ arethe shape functions.
+#
+# Based on the defintion of the position vector, we create an function for obtaining the Jacobian-matrix,
+# ```math
+# J_{ij} = \frac{\partial x_i}{\partial \xi_j},
+# ```
+function getjacobian(q, N, dNdξ, ζ, X, p, h)
+
+    J = zeros(3,3)
+    for a in 1:length(N)
+        for i in 1:3, j in 1:3
+            _dNdξ = (j==3) ? 0.0 : dNdξ[a][j]
+            _dζdξ = (j==3) ? 1.0 : 0.0
+            _N = N[a]
+
+            J[i,j] += _dNdξ * X[a][i]  +  (_dNdξ*ζ + _N*_dζdξ) * h/2 * p[a][i]
+        end
+    end
+
+    return (q' * J) |> Tensor{2,3,Float64}
+end;
+
+# ##### Strains
+# Small deformation is assumed,
+# ```math
+# \varepsilon_{ij}= \frac{1}{2}(\frac{\partial u_{i}}{\partial x_j} + \frac{\partial u_{j}}{\partial x_i})
+# ```
+# The displacement field is calculated as:
+# ```math
+# \boldsymbol u = \sum_{a=1}^{N_{\text{nodes}}} N_a \bar{\boldsymbol u}_{a} + 
+#  N_a ζ\frac{h}{2}(\theta_{a2} \boldsymbol e^{f}_{a1} - \theta_{a1} \boldsymbol e^{f}_{a2})
+#
+# ```
+# The gradient of the displacement (in the lamina coordinate system), then becomes:
+# ```math
+# \frac{\partial u_{i}}{\partial x_j} = \sum_{m=1}^3 q_{im} \sum_{a=1}^{N_{\text{nodes}}} \frac{\partial N_a}{\partial x_j} \bar{u}_{am} + 
+#  \frac{\partial(N_a ζ)}{\partial x_j} \frac{h}{2} (\theta_{a2} e^{f}_{am1} - \theta_{a1} e^{f}_{am2})
+# ```
+function strain(dofvec::Vector{T}, N, dNdx, ζ, dζdx, q, ef1, ef2, h) where T
+
+    u = reinterpret(Vec{3,T}, dofvec[1:12])
+    θ = reinterpret(Vec{2,T}, dofvec[13:20])
+
+    dudx = zeros(T, 3, 3)
+    for m in 1:3, j in 1:3
+        for a in 1:length(N)
+            dudx[m,j] += dNdx[a][j] * u[a][m] + h/2 * (dNdx[a][j]*ζ + N[a]*dζdx[j]) * (θ[a][2]*ef1[a][m] - θ[a][1]*ef2[a][m])
+        end
+    end
+
+    dudx = q*dudx
+    ε = [dudx[1,1], dudx[2,2], dudx[1,2]+dudx[2,1], dudx[2,3]+dudx[3,2], dudx[1,3]+dudx[3,1]]
+    return ε
+
+end;
+
+# ##### Main element routine
+# Below is the main routine that calculates the stiffness matrix of the shell element. 
+# Since it is a so called degenerate shell element, the code is similar to that for an standard continuum element.
+function integrate_shell!(ke, cv, qr_ooplane, X, data)
+    nnodes = getnbasefunctions(cv)
+    ndofs = nnodes*5
+    h = data.thickness
+
+    #Create the directors in each node.
+    #Note: For a more general case, the directors should
+    #be input parameters for the element routine.
+    p = zeros(Vec{3}, nnodes)
+    for i in 1:nnodes
+        a = Vec{3}((0.0, 0.0, 1.0))
+        p[i] = a/norm(a)
+    end
+    
+    ef1, ef2, ef3 = fiber_coordsys(p)
+
+    for iqp in 1:getnquadpoints(cv)
+
+        dNdξ = cv.dNdξ[:,iqp]
+        N = cv.N[:,iqp]
+
+        for oqp in 1:length(qr_ooplane.weights)
+
+            ζ = qr_ooplane.points[oqp][1]
+           
+            q = lamina_coordsys(dNdξ, ζ, X, p, h)
+            J = getjacobian(q, N, dNdξ, ζ, X, p, h)
+            Jinv = inv(J)  
+           
+            dζdx = Vec{3}((0.0, 0.0, 1.0)) ⋅ Jinv
+            dNdx = [Vec{3}((dNdξ[i][1], dNdξ[i][2], 0.0)) ⋅ Jinv for i in 1:nnodes]
+
+            
+            #For simplicity, use automatic differentiation to construct the B-matrix from the strain.
+            B = ForwardDiff.jacobian(
+                (a) -> strain(a, N, dNdx, ζ, dζdx, q, ef1, ef2, h), zeros(Float64, ndofs) )
+
+            dV = det(J) * cv.qr_weights[iqp] * qr_ooplane.weights[oqp]
+            ke .+= B'*data.C*B * dV
+        end
+    end
+end;
+
+# Run everything:
+main()
\ No newline at end of file
diff --git a/docs/src/literate/linear_shell.png b/docs/src/literate/linear_shell.png
new file mode 100644
index 0000000000..5693696f42
Binary files /dev/null and b/docs/src/literate/linear_shell.png differ
diff --git a/src/Dofs/ConstraintHandler.jl b/src/Dofs/ConstraintHandler.jl
index f053f8c5d2..5f4382ea4c 100644
--- a/src/Dofs/ConstraintHandler.jl
+++ b/src/Dofs/ConstraintHandler.jl
@@ -28,16 +28,16 @@ which applies the condition via `apply!`.
 """
 struct Dirichlet # <: Constraint
     f::Function # f(x,t) -> value
-    faces::Union{Set{Int},Set{Tuple{Int,Int}}}
+    faces::Union{Set{Int},Set{FaceIndex},Set{EdgeIndex},Set{VertexIndex}}
     field_name::Symbol
     components::Vector{Int} # components of the field
     local_face_dofs::Vector{Int}
     local_face_dofs_offset::Vector{Int}
 end
-function Dirichlet(field_name::Symbol, faces::Union{Set{Int},Set{Tuple{Int,Int}}}, f::Function, component::Int=1)
+function Dirichlet(field_name::Symbol, faces::Union{T}, f::Function, component::Int=1) where T
     Dirichlet(field_name, copy(faces), f, [component])
 end
-function Dirichlet(field_name::Symbol, faces::Union{Set{Int},Set{Tuple{Int,Int}}}, f::Function, components::AbstractVector{Int})
+function Dirichlet(field_name::Symbol, faces::Set{T}, f::Function, components::AbstractVector{Int}) where T
     unique(components) == components || error("components not unique: $components")
     # issorted(components) || error("components not sorted: $components")
     return Dirichlet(f, copy(faces), field_name, Vector(components), Int[], Int[])
@@ -164,21 +164,32 @@ Add a `Dirichlet boundary` condition to the `ConstraintHandler`.
 """
 function add!(ch::ConstraintHandler, dbc::Dirichlet)
     dbc_check(ch, dbc)
+    celltype = getcelltype(ch.dh.grid)
+    @assert isconcretetype(celltype)
+
     field_idx = find_field(ch.dh, dbc.field_name)
     # Extract stuff for the field
     interpolation = getfieldinterpolation(ch.dh, field_idx)#ch.dh.field_interpolations[field_idx]
     field_dim = getfielddim(ch.dh, field_idx)#ch.dh.field_dims[field_idx] # TODO: I think we don't need to extract these here ...
-    bcvalue = getbcvalue(ch.dh, field_idx)
+    
+    
+    if eltype(dbc.faces)==Int #Special case when dbc.faces is a nodeset
+        bcvalue = BCValues(interpolation, default_interpolation(celltype), FaceIndex) #Not used by node bcs, but still have to pass it as an argument
+    else
+        bcvalue = BCValues(interpolation, default_interpolation(celltype), eltype(dbc.faces))
+    end
     _add!(ch, dbc, dbc.faces, interpolation, field_dim, field_offset(ch.dh, dbc.field_name), bcvalue)
+
     return ch
 end
 
-function _add!(ch::ConstraintHandler, dbc::Dirichlet, bcfaces::Set{Tuple{Int,Int}}, interpolation::Interpolation, field_dim::Int, offset::Int, bcvalue::BCValues, cellset::Set{Int}=Set{Int}(1:getncells(ch.dh.grid)))
+function _add!(ch::ConstraintHandler, dbc::Dirichlet, bcfaces::Set{Index}, interpolation::Interpolation, field_dim::Int, offset::Int, bcvalue::BCValues, cellset::Set{Int}=Set{Int}(1:getncells(ch.dh.grid))) where Index<:BoundaryIndex
     # calculate which local dof index live on each face
     # face `i` have dofs `local_face_dofs[local_face_dofs_offset[i]:local_face_dofs_offset[i+1]-1]
     local_face_dofs = Int[]
     local_face_dofs_offset = Int[1]
-    for (i, face) in enumerate(faces(interpolation))
+    boundary = boundaryfunction(eltype(bcfaces))
+    for (i, face) in enumerate(boundary(interpolation))
         for fdof in face, d in 1:field_dim
             if d ∈ dbc.components # skip unless this component should be constrained
                 push!(local_face_dofs, (fdof-1)*field_dim + d + offset)
@@ -191,10 +202,10 @@ function _add!(ch::ConstraintHandler, dbc::Dirichlet, bcfaces::Set{Tuple{Int,Int
 
     # loop over all the faces in the set and add the global dofs to `constrained_dofs`
     constrained_dofs = Int[]
-    redundant_faces = NTuple{2, Int}[]
+    redundant_faces = Index[]
     for (cellidx, faceidx) in bcfaces
         if cellidx ∉ cellset
-            push!(redundant_faces, (cellidx, faceidx)) # will be removed from dbc
+            push!(redundant_faces, Index(cellidx, faceidx)) # will be removed from dbc
             continue # skip faces that are not part of the cellset
         end
         _celldofs = fill(0, ndofs_per_cell(ch.dh, cellidx))
@@ -242,7 +253,7 @@ function _add!(ch::ConstraintHandler, dbc::Dirichlet, bcnodes::Set{Int}, interpo
     sizehint!(dbc.local_face_dofs, length(bcnodes))
     for node in bcnodes
         if !visited[node]
-    # either the node belongs to another field handler or it does not have dofs in the constrained field
+            # either the node belongs to another field handler or it does not have dofs in the constrained field
             continue
         end
         for i in 1:ncomps
@@ -269,7 +280,7 @@ function update!(ch::ConstraintHandler, time::Real=0.0)
 end
 
 # for faces
-function _update!(values::Vector{Float64}, f::Function, faces::Set{Tuple{Int,Int}}, field::Symbol, local_face_dofs::Vector{Int}, local_face_dofs_offset::Vector{Int},
+function _update!(values::Vector{Float64}, f::Function, faces::Set{<:BoundaryIndex}, field::Symbol, local_face_dofs::Vector{Int}, local_face_dofs_offset::Vector{Int},
                   components::Vector{Int}, dh::AbstractDofHandler, facevalues::BCValues,
                   dofmapping::Dict{Int,Int}, time::T) where {T}
 
@@ -342,9 +353,10 @@ function WriteVTK.vtk_point_data(vtkfile, ch::ConstraintHandler)
         data = zeros(Float64, nd, getnnodes(ch.dh.grid))
         for dbc in ch.dbcs
             dbc.field_name != field && continue
-            if eltype(dbc.faces) <: Tuple
+            if eltype(dbc.faces) <: BoundaryIndex
+                functype = boundaryfunction(eltype(dbc.faces))
                 for (cellidx, faceidx) in dbc.faces
-                    for facenode in faces(ch.dh.grid.cells[cellidx])[faceidx]
+                    for facenode in functype(ch.dh.grid.cells[cellidx])[faceidx]
                         for component in dbc.components
                             data[component, facenode] = 1
                         end
@@ -449,17 +461,24 @@ end
 function add!(ch::ConstraintHandler, fh::FieldHandler, dbc::Dirichlet)
     _check_cellset_dirichlet(ch.dh.grid, fh.cellset, dbc.faces)
 
+    celltype = getcelltype(ch.dh.grid, first(fh.cellset)) #Assume same celltype of all cells in fh.cellset
+
     field_idx = find_field(fh, dbc.field_name)
     # Extract stuff for the field
     interpolation = getfieldinterpolations(fh)[field_idx]
     field_dim = getfielddims(fh)[field_idx]
-    bcvalue = fh.bc_values[field_idx]
+
+    if eltype(dbc.faces)==Int #Special case when dbc.faces is a nodeset
+        bcvalue = BCValues(interpolation, default_interpolation(celltype), FaceIndex) #Not used by node bcs, but still have to pass it as an argument
+    else
+        bcvalue = BCValues(interpolation, default_interpolation(celltype), eltype(dbc.faces))
+    end
 
     Ferrite._add!(ch, dbc, dbc.faces, interpolation, field_dim, field_offset(fh, dbc.field_name), bcvalue, fh.cellset)
     return ch
 end
 
-function _check_cellset_dirichlet(::AbstractGrid, cellset::Set{Int}, faceset::Set{Tuple{Int,Int}})
+function _check_cellset_dirichlet(::AbstractGrid, cellset::Set{Int}, faceset::Set{<:BoundaryIndex})
     for (cellid,faceid) in faceset
         if !(cellid in cellset)
             @warn("You are trying to add a constraint to a face that is not in the cellset of the fieldhandler. The face will be skipped.")
diff --git a/src/Dofs/DofHandler.jl b/src/Dofs/DofHandler.jl
index 4b09ed05d6..01d5f2b931 100644
--- a/src/Dofs/DofHandler.jl
+++ b/src/Dofs/DofHandler.jl
@@ -17,7 +17,7 @@ struct DofHandler{dim,C,T} <: AbstractDofHandler
     # TODO: field_interpolations can probably be better typed: We should at least require
     #       all the interpolations to have the same dimension and reference shape
     field_interpolations::Vector{Interpolation}
-    bc_values::Vector{BCValues{T}} # for boundary conditions
+    bc_values::Vector{BCValues{T}} # TODO: BcValues is created/handeld by the constrainthandler, so this can be removed
     cell_dofs::Vector{Int}
     cell_dofs_offset::Vector{Int}
     closed::ScalarWrapper{Bool}
@@ -26,6 +26,7 @@ struct DofHandler{dim,C,T} <: AbstractDofHandler
 end
 
 function DofHandler(grid::Grid)
+    isconcretetype(getcelltype(grid)) || error("Grid includes different celltypes. Use MixedDofHandler instead of DofHandler")
     DofHandler(Symbol[], Int[], Interpolation[], BCValues{Float64}[], Int[], Int[], ScalarWrapper(false), grid, Ferrite.ScalarWrapper(-1))
 end
 
diff --git a/src/Dofs/MixedDofHandler.jl b/src/Dofs/MixedDofHandler.jl
index 6a25e1f03c..a83568fb2b 100644
--- a/src/Dofs/MixedDofHandler.jl
+++ b/src/Dofs/MixedDofHandler.jl
@@ -2,7 +2,6 @@
 mutable struct FieldHandler
     fields::Vector{Field}
     cellset::Set{Int}
-    bc_values::Vector{BCValues} # for boundary conditions
 end
 """
     FieldHandler(fields::Vector{Field}, cellset)
@@ -10,9 +9,7 @@ end
 Construct a `FieldHandler` based on an array of `Field`s and assigns it a set of cells.
 """
 function FieldHandler(fields::Vector{Field}, cellset)
-    # TODO for now, only accept isoparamtric mapping
-    bc_values = [BCValues(field.interpolation, field.interpolation) for field in fields]
-    FieldHandler(fields, cellset, bc_values)
+    return FieldHandler(fields, cellset)
 end
 
 struct CellVector{T}
@@ -138,19 +135,20 @@ end
 function Base.push!(dh::MixedDofHandler, name::Symbol, dim::Int, ip::Interpolation)
     @assert !isclosed(dh)
 
+    celltype = getcelltype(dh.grid)
+    @assert isconcretetype(celltype)
+
     if length(dh.fieldhandlers) == 0
         cellset = Set(1:getncells(dh.grid))
-        push!(dh.fieldhandlers, FieldHandler(Field[],cellset,BCValues[]))
+        push!(dh.fieldhandlers, FieldHandler(Field[], cellset))
     elseif length(dh.fieldhandlers) > 1
         error("If you have more than one FieldHandler, you must specify field")
     end
     fh = first(dh.fieldhandlers)
 
     field = Field(name,ip,dim)
-    bc_value = BCValues(field.interpolation, field.interpolation)
 
     push!(fh.fields, field)
-    push!(fh.bc_values, bc_value)
 
     return dh
 end
@@ -414,4 +412,3 @@ find_field(dh::MixedDofHandler, field_name::Symbol) = find_field(first(dh.fieldh
 field_offset(dh::MixedDofHandler, field_name::Symbol) = field_offset(first(dh.fieldhandlers), field_name)
 getfieldinterpolation(dh::MixedDofHandler, field_idx::Int) = dh.fieldhandlers[1].fields[field_idx].interpolation
 getfielddim(dh::MixedDofHandler, field_idx::Int) = dh.fieldhandlers[1].fields[field_idx].dim
-getbcvalue(dh::MixedDofHandler, field_idx::Int) =  dh.fieldhandlers[1].bc_values[field_idx]
diff --git a/src/Export/VTK.jl b/src/Export/VTK.jl
index 54979a4ed5..098f613796 100644
--- a/src/Export/VTK.jl
+++ b/src/Export/VTK.jl
@@ -1,7 +1,10 @@
 cell_to_vtkcell(::Type{Line}) = VTKCellTypes.VTK_LINE
+cell_to_vtkcell(::Type{Line2D}) = VTKCellTypes.VTK_LINE
+cell_to_vtkcell(::Type{Line3D}) = VTKCellTypes.VTK_LINE
 cell_to_vtkcell(::Type{QuadraticLine}) = VTKCellTypes.VTK_QUADRATIC_EDGE
 
 cell_to_vtkcell(::Type{Quadrilateral}) = VTKCellTypes.VTK_QUAD
+cell_to_vtkcell(::Type{Quadrilateral3D}) = VTKCellTypes.VTK_QUAD
 cell_to_vtkcell(::Type{QuadraticQuadrilateral}) = VTKCellTypes.VTK_BIQUADRATIC_QUAD
 cell_to_vtkcell(::Type{Triangle}) = VTKCellTypes.VTK_TRIANGLE
 cell_to_vtkcell(::Type{QuadraticTriangle}) = VTKCellTypes.VTK_QUADRATIC_TRIANGLE
diff --git a/src/FEValues/face_values.jl b/src/FEValues/face_values.jl
index 428494c601..df0f88f97e 100644
--- a/src/FEValues/face_values.jl
+++ b/src/FEValues/face_values.jl
@@ -208,22 +208,27 @@ Return the normal at the quadrature point `qp` for the active face of the
 """
 getnormal(fv::FaceValues, qp::Int) = fv.normals[qp]
 
-# like FaceScalarValues but contains only the parts needed
-# to calculate the x-coordinate for the dof locations.
+"""
+    BCValues(func_interpol::Interpolation, geom_interpol::Interpolation, boundary_type::Union{Type{<:BoundaryIndex}})
+
+`BCValues` stores the shape values at all faces/edges/vertices (depending on `boundary_type`) for the geomatric interpolation (`geom_interpol`), 
+for each dof-position determined by the `func_interpol`. Used mainly by the `ConstrainHandler`.
+"""
 struct BCValues{T}
     M::Array{T,3}
     current_face::ScalarWrapper{Int}
 end
 
-BCValues(func_interpol::Interpolation, geom_interpol::Interpolation) =
-    BCValues(Float64, func_interpol, geom_interpol)
+BCValues(func_interpol::Interpolation, geom_interpol::Interpolation, boundary_type::Type{<:BoundaryIndex} = Ferrite.FaceIndex) =
+    BCValues(Float64, func_interpol, geom_interpol, boundary_type)
 
-function BCValues(::Type{T}, func_interpol::Interpolation{dim,refshape}, geom_interpol::Interpolation{dim,refshape}) where {T,dim,refshape}
+function BCValues(::Type{T}, func_interpol::Interpolation{dim,refshape}, geom_interpol::Interpolation{dim,refshape}, boundary_type::Type{<:BoundaryIndex} = Ferrite.FaceIndex) where {T,dim,refshape}
     # set up quadrature rules for each face with dof-positions
     # (determined by func_interpol) as the quadrature points
     interpolation_coords = reference_coordinates(func_interpol)
 
     qrs = QuadratureRule{dim,refshape,T}[]
+    faces = boundaryfunction(boundary_type) # faces, edges or vertices
     for face in faces(func_interpol)
         dofcoords = Vec{dim,T}[]
         for facedof in face
diff --git a/src/Ferrite.jl b/src/Ferrite.jl
index e4e67d2715..d8054f97f1 100644
--- a/src/Ferrite.jl
+++ b/src/Ferrite.jl
@@ -27,6 +27,11 @@ abstract type Values{dim,T,refshape} end
 abstract type CellValues{dim,T,refshape} <: Values{dim,T,refshape} end
 abstract type FaceValues{dim,T,refshape} <: Values{dim,T,refshape} end
 
+"""
+Abstract type which is used as identifier for faces, edges and verices
+"""
+abstract type BoundaryIndex end
+
 include("utils.jl")
 
 # Interpolations
diff --git a/src/Grid/grid.jl b/src/Grid/grid.jl
index 9a903fe122..849ea4ca04 100644
--- a/src/Grid/grid.jl
+++ b/src/Grid/grid.jl
@@ -51,7 +51,21 @@ end
 """
 A `FaceIndex` wraps an (Int, Int) and defines a face by pointing to a (cell, face).
 """
-struct FaceIndex
+struct FaceIndex <: BoundaryIndex
+    idx::Tuple{Int,Int} # cell and side
+end
+
+"""
+A `EdgeIndex` wraps an (Int, Int) and defines a face by pointing to a (cell, edge).
+"""
+struct EdgeIndex <: BoundaryIndex
+    idx::Tuple{Int,Int} # cell and side
+end
+
+"""
+A `VertexIndex` wraps an (Int, Int) and defines a face by pointing to a (cell, vert).
+"""
+struct VertexIndex <: BoundaryIndex
     idx::Tuple{Int,Int} # cell and side
 end
 
@@ -66,7 +80,9 @@ mutable struct Grid{dim,C<:AbstractCell,T<:Real} <: AbstractGrid{dim}
     # Sets
     cellsets::Dict{String,Set{Int}}
     nodesets::Dict{String,Set{Int}}
-    facesets::Dict{String,Set{Tuple{Int,Int}}} # TODO: This could be Set{FaceIndex} which could result in nicer use later
+    facesets::Dict{String,Set{FaceIndex}} 
+    edgesets::Dict{String,Set{EdgeIndex}} 
+    vertexsets::Dict{String,Set{VertexIndex}} 
     # Boundary matrix (faces per cell × cell)
     boundary_matrix::SparseMatrixCSC{Bool,Int}
 end
@@ -75,9 +91,11 @@ 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}}(),
-              facesets::Dict{String,Set{Tuple{Int,Int}}}=Dict{String,Set{Tuple{Int,Int}}}(),
+              facesets::Dict{String,Set{FaceIndex}}=Dict{String,Set{FaceIndex}}(),
+              edgesets::Dict{String,Set{EdgeIndex}}=Dict{String,Set{EdgeIndex}}(),
+              vertexsets::Dict{String,Set{VertexIndex}}=Dict{String,Set{VertexIndex}}(),
               boundary_matrix::SparseMatrixCSC{Bool,Int}=spzeros(Bool, 0, 0)) where {dim,C,T}
-    return Grid(cells, nodes, cellsets, nodesets, facesets, boundary_matrix)
+    return Grid(cells, nodes, cellsets, nodesets, facesets, edgesets, vertexsets, boundary_matrix)
 end
 
 ##########################
@@ -89,6 +107,7 @@ end
 @inline getcells(grid::AbstractGrid, set::String) = grid.cells[collect(grid.cellsets[set])]
 @inline getncells(grid::AbstractGrid) = length(grid.cells)
 @inline getcelltype(grid::AbstractGrid) = eltype(grid.cells)
+@inline getcelltype(grid::AbstractGrid, i::Int) = typeof(grid.cells[i])
 
 @inline getnodes(grid::AbstractGrid) = grid.nodes
 @inline getnodes(grid::AbstractGrid, v::Union{Int, Vector{Int}}) = grid.nodes[v]
@@ -105,6 +124,12 @@ end
 @inline getfaceset(grid::AbstractGrid, set::String) = grid.facesets[set]
 @inline getfacesets(grid::AbstractGrid) = grid.facesets
 
+@inline getedgeset(grid::AbstractGrid, set::String) = grid.edgesets[set]
+@inline getedgesets(grid::AbstractGrid) = grid.edgesets
+
+@inline getvertexset(grid::AbstractGrid, set::String) = grid.vertexsets[set]
+@inline getvertexsets(grid::AbstractGrid) = grid.vertexsets
+
 n_faces_per_cell(grid::Grid) = nfaces(eltype(grid.cells))
 
 
@@ -165,29 +190,41 @@ function addcellset!(grid::AbstractGrid, name::String, f::Function; all::Bool=tr
     grid
 end
 
-function addfaceset!(grid::AbstractGrid, name::String, faceid::Set{Tuple{Int,Int}})
-    _check_setname(grid.facesets, name)
-    faceset = Set(faceid)
-    _warn_emptyset(faceset)
-    grid.facesets[name] = faceset
+addfaceset!(grid::Grid, name::String, set::Union{Set{FaceIndex},Vector{FaceIndex}}) = 
+    _addset!(grid, name, set, grid.facesets)
+addedgeset!(grid::Grid, name::String, set::Union{Set{EdgeIndex},Vector{EdgeIndex}}) = 
+    _addset!(grid, name, set, grid.edgesets)
+addvertexset!(grid::Grid, name::String, set::Union{Set{VertexIndex},Vector{VertexIndex}}) = 
+    _addset!(grid, name, set, grid.vertexsets)
+function _addset!(grid::AbstractGrid, name::String, _set, dict::Dict)
+    _check_setname(dict, name)
+    set = Set(_set)
+    _warn_emptyset(set)
+    dict[name] = set
     grid
 end
 
-function addfaceset!(grid::AbstractGrid, name::String, f::Function; all::Bool=true)
-    _check_setname(grid.facesets, name)
-    faceset = Set{Tuple{Int,Int}}()
+addfaceset!(grid::AbstractGrid, name::String, f::Function; all::Bool=true) = 
+    _addset!(grid, name, f, Ferrite.faces, grid.facesets, FaceIndex; all=all)
+addedgeset!(grid::AbstractGrid, name::String, f::Function; all::Bool=true) = 
+    _addset!(grid, name, f, Ferrite.edges, grid.edgesets, EdgeIndex; all=all)
+addvertexset!(grid::AbstractGrid, name::String, f::Function; all::Bool=true) = 
+    _addset!(grid, name, f, Ferrite.vertices, grid.vertexsets, VertexIndex; all=all)
+function _addset!(grid::AbstractGrid, name::String, f::Function, _ftype::Function, dict::Dict, _indextype::Type; all::Bool=true)
+    _check_setname(dict, name)
+    _set = Set{_indextype}()
     for (cell_idx, cell) in enumerate(getcells(grid))
-        for (face_idx, face) in enumerate(faces(cell))
+        for (face_idx, face) in enumerate(_ftype(cell))
             pass = all
             for node_idx in face
                 v = f(grid.nodes[node_idx].x)
                 all ? (!v && (pass = false; break)) : (v && (pass = true; break))
             end
-            pass && push!(faceset, (cell_idx, face_idx))
+            pass && push!(_set, _indextype(cell_idx, face_idx))
         end
     end
-    _warn_emptyset(faceset)
-    grid.facesets[name] = faceset
+    _warn_emptyset(_set)
+    dict[name] = _set
     grid
 end
 
@@ -300,3 +337,22 @@ default_interpolation(::Type{Tetrahedron}) = Lagrange{3,RefTetrahedron,1}()
 default_interpolation(::Type{QuadraticTetrahedron}) = Lagrange{3,RefTetrahedron,2}()
 default_interpolation(::Type{Hexahedron}) = Lagrange{3,RefCube,1}()
 default_interpolation(::Type{QuadraticHexahedron}) = Lagrange{3,RefCube,2}()
+
+boundaryfunction(::Type{FaceIndex}) = Ferrite.faces
+boundaryfunction(::Type{EdgeIndex}) = Ferrite.edges
+boundaryfunction(::Type{VertexIndex}) = Ferrite.vertices
+
+for INDEX in (:VertexIndex, :EdgeIndex, :FaceIndex)
+    @eval begin  
+        #Constructor
+        ($INDEX)(a::Int, b::Int) = ($INDEX)((a,b))
+
+        Base.getindex(I::($INDEX), i::Int) = I.idx[i]
+        
+        #To be able to do a,b = faceidx
+        Base.iterate(I::($INDEX), state::Int=1) = (state==3) ?  nothing : (I[state], state+1)
+
+        #For (cellid, faceidx) in faceset
+        Base.in(v::Tuple{Int, Int}, s::Set{$INDEX}) = in($INDEX(v), s)
+    end
+end
\ No newline at end of file
diff --git a/src/Grid/grid_generators.jl b/src/Grid/grid_generators.jl
index 7161010058..9067c6668d 100644
--- a/src/Grid/grid_generators.jl
+++ b/src/Grid/grid_generators.jl
@@ -4,7 +4,7 @@ function boundaries_to_sparse(boundary)
     J = Vector{Int}(undef, n)
     V = Vector{Bool}(undef, n)
     for (idx, el) in enumerate(boundary)
-        cell, face = el
+        cell, face = el.idx
         I[idx] = face
         J[idx] = cell
         V[idx] = true
@@ -42,14 +42,14 @@ function generate_grid(::Type{Line}, nel::NTuple{1,Int}, left::Vec{1,T}=Vec{1}((
 
 
     # Cell faces
-    boundary = Vector([(1, 1),
-                       (nel_x, 2)])
+    boundary = Vector([FaceIndex(1, 1),
+                       FaceIndex(nel_x, 2)])
 
     boundary_matrix = boundaries_to_sparse(boundary)
 
     # Cell face sets
-    facesets = Dict("left"  => Set{Tuple{Int,Int}}([boundary[1]]),
-                    "right" => Set{Tuple{Int,Int}}([boundary[2]]))
+    facesets = Dict("left"  => Set{FaceIndex}([boundary[1]]),
+                    "right" => Set{FaceIndex}([boundary[2]]))
     return Grid(cells, nodes, facesets=facesets, boundary_matrix=boundary_matrix)
 end
 
@@ -72,14 +72,14 @@ function generate_grid(::Type{QuadraticLine}, nel::NTuple{1,Int}, left::Vec{1,T}
     end
 
     # Cell faces
-    boundary = Tuple{Int,Int}[(1, 1),
-                         (nel_x, 2)]
+    boundary = FaceIndex[FaceIndex(1, 1),
+                         FaceIndex(nel_x, 2)]
 
     boundary_matrix = boundaries_to_sparse(boundary)
 
     # Cell face sets
-    facesets = Dict("left"  => Set{Tuple{Int,Int}}([boundary[1]]),
-                    "right" => Set{Tuple{Int,Int}}([boundary[2]]))
+    facesets = Dict("left"  => Set{FaceIndex}([boundary[1]]),
+                    "right" => Set{FaceIndex}([boundary[2]]))
     return Grid(cells, nodes, facesets=facesets, boundary_matrix=boundary_matrix)
 end
 
@@ -135,20 +135,20 @@ function generate_grid(C::Type{Quadrilateral}, nel::NTuple{2,Int}, LL::Vec{2,T},
 
     # Cell faces
     cell_array = reshape(collect(1:nel_tot),(nel_x, nel_y))
-    boundary = Tuple{Int,Int}[[(cl, 1) for cl in cell_array[:,1]];
-                              [(cl, 2) for cl in cell_array[end,:]];
-                              [(cl, 3) for cl in cell_array[:,end]];
-                              [(cl, 4) for cl in cell_array[1,:]]]
+    boundary = FaceIndex[[FaceIndex(cl, 1) for cl in cell_array[:,1]];
+                              [FaceIndex(cl, 2) for cl in cell_array[end,:]];
+                              [FaceIndex(cl, 3) for cl in cell_array[:,end]];
+                              [FaceIndex(cl, 4) for cl in cell_array[1,:]]]
 
     boundary_matrix = boundaries_to_sparse(boundary)
 
     # Cell face sets
     offset = 0
-    facesets = Dict{String, Set{Tuple{Int,Int}}}()
-    facesets["bottom"] = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[:,1]))   .+ offset]); offset += length(cell_array[:,1])
-    facesets["right"]  = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[end,:])) .+ offset]); offset += length(cell_array[end,:])
-    facesets["top"]    = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[:,end])) .+ offset]); offset += length(cell_array[:,end])
-    facesets["left"]   = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[1,:]))   .+ offset]); offset += length(cell_array[1,:])
+    facesets = Dict{String, Set{FaceIndex}}()
+    facesets["bottom"] = Set{FaceIndex}(boundary[(1:length(cell_array[:,1]))   .+ offset]); offset += length(cell_array[:,1])
+    facesets["right"]  = Set{FaceIndex}(boundary[(1:length(cell_array[end,:])) .+ offset]); offset += length(cell_array[end,:])
+    facesets["top"]    = Set{FaceIndex}(boundary[(1:length(cell_array[:,end])) .+ offset]); offset += length(cell_array[:,end])
+    facesets["left"]   = Set{FaceIndex}(boundary[(1:length(cell_array[1,:]))   .+ offset]); offset += length(cell_array[1,:])
 
     return Grid(cells, nodes, facesets=facesets, boundary_matrix=boundary_matrix)
 end
@@ -174,20 +174,20 @@ function generate_grid(::Type{QuadraticQuadrilateral}, nel::NTuple{2,Int}, LL::V
 
     # Cell faces
     cell_array = reshape(collect(1:nel_tot),(nel_x, nel_y))
-    boundary = Tuple{Int,Int}[[(cl, 1) for cl in cell_array[:,1]];
-                              [(cl, 2) for cl in cell_array[end,:]];
-                              [(cl, 3) for cl in cell_array[:,end]];
-                              [(cl, 4) for cl in cell_array[1,:]]]
+    boundary = FaceIndex[[FaceIndex(cl, 1) for cl in cell_array[:,1]];
+                              [FaceIndex(cl, 2) for cl in cell_array[end,:]];
+                              [FaceIndex(cl, 3) for cl in cell_array[:,end]];
+                              [FaceIndex(cl, 4) for cl in cell_array[1,:]]]
 
     boundary_matrix = boundaries_to_sparse(boundary)
 
     # Cell face sets
     offset = 0
-    facesets = Dict{String, Set{Tuple{Int,Int}}}()
-    facesets["bottom"] = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[:,1]))   .+ offset]); offset += length(cell_array[:,1])
-    facesets["right"]  = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[end,:])) .+ offset]); offset += length(cell_array[end,:])
-    facesets["top"]    = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[:,end])) .+ offset]); offset += length(cell_array[:,end])
-    facesets["left"]   = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[1,:]))   .+ offset]); offset += length(cell_array[1,:])
+    facesets = Dict{String, Set{FaceIndex}}()
+    facesets["bottom"] = Set{FaceIndex}(boundary[(1:length(cell_array[:,1]))   .+ offset]); offset += length(cell_array[:,1])
+    facesets["right"]  = Set{FaceIndex}(boundary[(1:length(cell_array[end,:])) .+ offset]); offset += length(cell_array[end,:])
+    facesets["top"]    = Set{FaceIndex}(boundary[(1:length(cell_array[:,end])) .+ offset]); offset += length(cell_array[:,end])
+    facesets["left"]   = Set{FaceIndex}(boundary[(1:length(cell_array[1,:]))   .+ offset]); offset += length(cell_array[1,:])
 
     return Grid(cells, nodes, facesets=facesets, boundary_matrix=boundary_matrix)
 end
@@ -217,24 +217,24 @@ function generate_grid(::Type{Hexahedron}, nel::NTuple{3,Int}, left::Vec{3,T}=Ve
 
     # Cell faces
     cell_array = reshape(collect(1:nel_tot),(nel_x, nel_y, nel_z))
-    boundary = Tuple{Int,Int}[[(cl, 1) for cl in cell_array[:,:,1][:]];
-                              [(cl, 2) for cl in cell_array[:,1,:][:]];
-                              [(cl, 3) for cl in cell_array[end,:,:][:]];
-                              [(cl, 4) for cl in cell_array[:,end,:][:]];
-                              [(cl, 5) for cl in cell_array[1,:,:][:]];
-                              [(cl, 6) for cl in cell_array[:,:,end][:]]]
+    boundary = FaceIndex[[FaceIndex(cl, 1) for cl in cell_array[:,:,1][:]];
+                              [FaceIndex(cl, 2) for cl in cell_array[:,1,:][:]];
+                              [FaceIndex(cl, 3) for cl in cell_array[end,:,:][:]];
+                              [FaceIndex(cl, 4) for cl in cell_array[:,end,:][:]];
+                              [FaceIndex(cl, 5) for cl in cell_array[1,:,:][:]];
+                              [FaceIndex(cl, 6) for cl in cell_array[:,:,end][:]]]
 
     boundary_matrix = boundaries_to_sparse(boundary)
 
     # Cell face sets
     offset = 0
-    facesets = Dict{String,Set{Tuple{Int,Int}}}()
-    facesets["bottom"] = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[:,:,1][:]))   .+ offset]); offset += length(cell_array[:,:,1][:])
-    facesets["front"]  = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[:,1,:][:]))   .+ offset]); offset += length(cell_array[:,1,:][:])
-    facesets["right"]  = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[end,:,:][:])) .+ offset]); offset += length(cell_array[end,:,:][:])
-    facesets["back"]   = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[:,end,:][:])) .+ offset]); offset += length(cell_array[:,end,:][:])
-    facesets["left"]   = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[1,:,:][:]))   .+ offset]); offset += length(cell_array[1,:,:][:])
-    facesets["top"]    = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[:,:,end][:])) .+ offset]); offset += length(cell_array[:,:,end][:])
+    facesets = Dict{String,Set{FaceIndex}}()
+    facesets["bottom"] = Set{FaceIndex}(boundary[(1:length(cell_array[:,:,1][:]))   .+ offset]); offset += length(cell_array[:,:,1][:])
+    facesets["front"]  = Set{FaceIndex}(boundary[(1:length(cell_array[:,1,:][:]))   .+ offset]); offset += length(cell_array[:,1,:][:])
+    facesets["right"]  = Set{FaceIndex}(boundary[(1:length(cell_array[end,:,:][:])) .+ offset]); offset += length(cell_array[end,:,:][:])
+    facesets["back"]   = Set{FaceIndex}(boundary[(1:length(cell_array[:,end,:][:])) .+ offset]); offset += length(cell_array[:,end,:][:])
+    facesets["left"]   = Set{FaceIndex}(boundary[(1:length(cell_array[1,:,:][:]))   .+ offset]); offset += length(cell_array[1,:,:][:])
+    facesets["top"]    = Set{FaceIndex}(boundary[(1:length(cell_array[:,:,end][:])) .+ offset]); offset += length(cell_array[:,:,end][:])
 
     return Grid(cells, nodes, facesets=facesets, boundary_matrix=boundary_matrix)
 end
@@ -259,20 +259,20 @@ function generate_grid(::Type{Triangle}, nel::NTuple{2,Int}, LL::Vec{2,T}, LR::V
 
     # Cell faces
     cell_array = reshape(collect(1:nel_tot),(2, nel_x, nel_y))
-    boundary = Tuple{Int,Int}[[(cl, 1) for cl in cell_array[1,:,1]];
-                               [(cl, 1) for cl in cell_array[2,end,:]];
-                               [(cl, 2) for cl in cell_array[2,:,end]];
-                               [(cl, 3) for cl in cell_array[1,1,:]]]
+    boundary = FaceIndex[[FaceIndex(cl, 1) for cl in cell_array[1,:,1]];
+                               [FaceIndex(cl, 1) for cl in cell_array[2,end,:]];
+                               [FaceIndex(cl, 2) for cl in cell_array[2,:,end]];
+                               [FaceIndex(cl, 3) for cl in cell_array[1,1,:]]]
 
     boundary_matrix = boundaries_to_sparse(boundary)
 
     # Cell face sets
     offset = 0
-    facesets = Dict{String,Set{Tuple{Int,Int}}}()
-    facesets["bottom"] = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[1,:,1]))   .+ offset]); offset += length(cell_array[1,:,1])
-    facesets["right"]  = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[2,end,:])) .+ offset]); offset += length(cell_array[2,end,:])
-    facesets["top"]    = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[2,:,end])) .+ offset]); offset += length(cell_array[2,:,end])
-    facesets["left"]   = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[1,1,:]))   .+ offset]); offset += length(cell_array[1,1,:])
+    facesets = Dict{String,Set{FaceIndex}}()
+    facesets["bottom"] = Set{FaceIndex}(boundary[(1:length(cell_array[1,:,1]))   .+ offset]); offset += length(cell_array[1,:,1])
+    facesets["right"]  = Set{FaceIndex}(boundary[(1:length(cell_array[2,end,:])) .+ offset]); offset += length(cell_array[2,end,:])
+    facesets["top"]    = Set{FaceIndex}(boundary[(1:length(cell_array[2,:,end])) .+ offset]); offset += length(cell_array[2,:,end])
+    facesets["left"]   = Set{FaceIndex}(boundary[(1:length(cell_array[1,1,:]))   .+ offset]); offset += length(cell_array[1,1,:])
 
     return Grid(cells, nodes, facesets=facesets, boundary_matrix=boundary_matrix)
 end
@@ -299,20 +299,20 @@ function generate_grid(::Type{QuadraticTriangle}, nel::NTuple{2,Int}, LL::Vec{2,
 
     # Cell faces
     cell_array = reshape(collect(1:nel_tot),(2, nel_x, nel_y))
-    boundary = Tuple{Int,Int}[[(cl, 1) for cl in cell_array[1,:,1]];
-                              [(cl, 1) for cl in cell_array[2,end,:]];
-                              [(cl, 2) for cl in cell_array[2,:,end]];
-                              [(cl, 3) for cl in cell_array[1,1,:]]]
+    boundary = FaceIndex[[FaceIndex(cl, 1) for cl in cell_array[1,:,1]];
+                              [FaceIndex(cl, 1) for cl in cell_array[2,end,:]];
+                              [FaceIndex(cl, 2) for cl in cell_array[2,:,end]];
+                              [FaceIndex(cl, 3) for cl in cell_array[1,1,:]]]
 
     boundary_matrix = boundaries_to_sparse(boundary)
 
     # Cell face sets
     offset = 0
-    facesets = Dict{String,Set{Tuple{Int,Int}}}()
-    facesets["bottom"] = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[1,:,1]))   .+ offset]); offset += length(cell_array[1,:,1])
-    facesets["right"]  = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[2,end,:])) .+ offset]); offset += length(cell_array[2,end,:])
-    facesets["top"]    = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[2,:,end])) .+ offset]); offset += length(cell_array[2,:,end])
-    facesets["left"]   = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[1,1,:]))   .+ offset]); offset += length(cell_array[1,1,:])
+    facesets = Dict{String,Set{FaceIndex}}()
+    facesets["bottom"] = Set{FaceIndex}(boundary[(1:length(cell_array[1,:,1]))   .+ offset]); offset += length(cell_array[1,:,1])
+    facesets["right"]  = Set{FaceIndex}(boundary[(1:length(cell_array[2,end,:])) .+ offset]); offset += length(cell_array[2,end,:])
+    facesets["top"]    = Set{FaceIndex}(boundary[(1:length(cell_array[2,:,end])) .+ offset]); offset += length(cell_array[2,:,end])
+    facesets["left"]   = Set{FaceIndex}(boundary[(1:length(cell_array[1,1,:]))   .+ offset]); offset += length(cell_array[1,1,:])
 
     return Grid(cells, nodes, facesets=facesets, boundary_matrix=boundary_matrix)
 end
@@ -376,12 +376,12 @@ function generate_grid(::Type{Tetrahedron}, cells_per_dim::NTuple{3,Int}, left::
     # Order the cells as c_nxyz[n, x, y, z] such that we can look up boundary cells
     c_nxyz = reshape(1:total_elements, (cells_per_cube, cells_per_dim...))
 
-    @views le = [map(x -> (x,4), c_nxyz[1, 1, :, :][:])   ; map(x -> (x,2), c_nxyz[2, 1, :, :][:])]
-    @views ri = [map(x -> (x,1), c_nxyz[4, end, :, :][:]) ; map(x -> (x,1), c_nxyz[6, end, :, :][:])]
-    @views fr = [map(x -> (x,1), c_nxyz[2, :, 1, :][:])   ; map(x -> (x,1), c_nxyz[5, :, 1, :][:])]
-    @views ba = [map(x -> (x,3), c_nxyz[3, :, end, :][:]) ; map(x -> (x,3), c_nxyz[4, :, end, :][:])]
-    @views bo = [map(x -> (x,1), c_nxyz[1, :, :, 1][:])   ; map(x -> (x,1), c_nxyz[3, :, :, 1][:])]
-    @views to = [map(x -> (x,3), c_nxyz[5, :, :, end][:]) ; map(x -> (x,3), c_nxyz[6, :, :, end][:])]
+    @views le = [map(x -> FaceIndex(x,4), c_nxyz[1, 1, :, :][:])   ; map(x -> FaceIndex(x,2), c_nxyz[2, 1, :, :][:])]
+    @views ri = [map(x -> FaceIndex(x,1), c_nxyz[4, end, :, :][:]) ; map(x -> FaceIndex(x,1), c_nxyz[6, end, :, :][:])]
+    @views fr = [map(x -> FaceIndex(x,1), c_nxyz[2, :, 1, :][:])   ; map(x -> FaceIndex(x,1), c_nxyz[5, :, 1, :][:])]
+    @views ba = [map(x -> FaceIndex(x,3), c_nxyz[3, :, end, :][:]) ; map(x -> FaceIndex(x,3), c_nxyz[4, :, end, :][:])]
+    @views bo = [map(x -> FaceIndex(x,1), c_nxyz[1, :, :, 1][:])   ; map(x -> FaceIndex(x,1), c_nxyz[3, :, :, 1][:])]
+    @views to = [map(x -> FaceIndex(x,3), c_nxyz[5, :, :, end][:]) ; map(x -> FaceIndex(x,3), c_nxyz[6, :, :, end][:])]
 
     boundary_matrix = boundaries_to_sparse([le; ri; bo; to; fr; ba])
 
diff --git a/src/exports.jl b/src/exports.jl
index 8f4dc23f4c..b70cd7437b 100644
--- a/src/exports.jl
+++ b/src/exports.jl
@@ -56,6 +56,8 @@ export
     QuadraticHexahedron,
     CellIndex,
     FaceIndex,
+    EdgeIndex,
+    VertexIndex,
     getcells,
     getncells,
     getnodes,
@@ -64,15 +66,21 @@ export
     getcellset,
     getnodeset,
     getfaceset,
+    getedgeset,
+    getvertexset,
     getcoordinates,
     getcoordinates!,
     getcellsets,
     getnodesets,
     getfacesets,
+    getedgesets,
+    getvertexsets,
     onboundary,
     nfaces,
     addnodeset!,
     addfaceset!,
+    addedgeset!,
+    addvertexset!,
     addcellset!,
     transform!,
     generate_grid,
diff --git a/src/interpolations.jl b/src/interpolations.jl
index 9fdab87371..07aad589ff 100644
--- a/src/interpolations.jl
+++ b/src/interpolations.jl
@@ -282,6 +282,7 @@ getnbasefunctions(::Lagrange{3,RefTetrahedron,1}) = 4
 nvertexdofs(::Lagrange{3,RefTetrahedron,1}) = 1
 
 faces(::Lagrange{3,RefTetrahedron,1}) = ((1,2,3), (1,2,4), (2,3,4), (1,4,3))
+edges(::Lagrange{3,RefTetrahedron,1}) = ((1,2), (2,3), (3,1), (1,4), (2,4), (3,4))
 
 function reference_coordinates(::Lagrange{3,RefTetrahedron,1})
     return [Vec{3, Float64}((0.0, 0.0, 0.0)),
@@ -309,6 +310,7 @@ nvertexdofs(::Lagrange{3,RefTetrahedron,2}) = 1
 nedgedofs(::Lagrange{3,RefTetrahedron,2}) = 1
 
 faces(::Lagrange{3,RefTetrahedron,2}) = ((1,2,3,5,6,7), (1,2,4,5,9,8), (2,3,4,6,10,9), (1,4,3,8,10,7))
+edges(::Lagrange{3,RefTetrahedron,2}) = ((1,5,2), (2,6,3), (3,7,1), (1,8,4), (2,9,4), (3,10,4))
 
 function reference_coordinates(::Lagrange{3,RefTetrahedron,2})
     return [Vec{3, Float64}((0.0, 0.0, 0.0)),
@@ -349,6 +351,7 @@ getnbasefunctions(::Lagrange{3,RefCube,1}) = 8
 nvertexdofs(::Lagrange{3,RefCube,1}) = 1
 
 faces(::Lagrange{3,RefCube,1}) = ((1,4,3,2), (1,2,6,5), (2,3,7,6), (3,4,8,7), (1,5,8,4), (5,6,7,8))
+edges(::Lagrange{3,RefCube,1}) = ((1,2), (2,3), (3,4), (4,1), (1,5), (2,6), (3,7), (4,8), (5,6), (6,7), (7,8), (8,5))
 
 function reference_coordinates(::Lagrange{3,RefCube,1})
     return [Vec{3, Float64}((-1.0, -1.0, -1.0)),
diff --git a/test/test_constraints.jl b/test/test_constraints.jl
index ff5cee07e7..8738f9255d 100644
--- a/test/test_constraints.jl
+++ b/test/test_constraints.jl
@@ -70,3 +70,46 @@
    @test ch.prescribed_dofs == [1,3,9,13,14]
    @test ch.values == [1.0, 1.0, 1.0, 2.0, 2.0]
 end
+
+@testset "edge bc" begin
+    grid = generate_grid(Hexahedron, (1, 1, 1))
+    addedgeset!(grid, "edge", x-> x[1] ≈ -1.0 && x[3] ≈ -1.0)
+
+    dh = DofHandler(grid)
+    push!(dh, :u, 3)
+    push!(dh, :p, 1)
+    close!(dh)
+
+    ch = ConstraintHandler(dh)
+    dbc1 = Dirichlet(:u, getedgeset(grid, "edge"), (x,t) -> x, [1, 2, 3])
+    add!(ch, dbc1)
+    close!(ch)
+    update!(ch)
+
+    @test ch.prescribed_dofs == [1,2,3,10,11,12]
+    @test ch.values == [-1.0, -1.0, -1.0, -1.0, 1.0, -1.0]
+
+
+    #Shell mesh edge bcs
+    nodes = [Node{3,Float64}(Vec(0.0,0.0,0.0)), Node{3,Float64}(Vec(1.0,0.0,0.0)), 
+             Node{3,Float64}(Vec(1.0,1.0,0.0)), Node{3,Float64}(Vec(0.0,1.0,0.0)),
+             Node{3,Float64}(Vec(2.0,0.0,0.0)), Node{3,Float64}(Vec(2.0,2.0,0.0))]
+
+    cells = [Quadrilateral3D((1,2,3,4)), Quadrilateral3D((2,5,6,3))]
+    grid = Grid(cells,nodes)
+
+    #3d quad with 1st order 2d interpolation
+    dh = DofHandler(grid)
+    push!(dh, :u, 1, Lagrange{2,RefCube,2}())
+    push!(dh, :θ, 1, Lagrange{2,RefCube,2}())
+    close!(dh)
+
+    addedgeset!(grid, "edge", x -> x[2] ≈ 0.0) #bottom edge
+    ch = ConstraintHandler(dh)
+    dbc1 = Dirichlet(:θ, getedgeset(grid, "edge"), (x,t) -> (0.0,), [1])
+    add!(ch, dbc1)
+    close!(ch)
+    update!(ch)
+
+    @test ch.prescribed_dofs == [10, 11, 14, 25, 27]
+end
\ No newline at end of file
diff --git a/test/test_dofs.jl b/test/test_dofs.jl
index 7cee0fba65..9e832b299b 100644
--- a/test/test_dofs.jl
+++ b/test/test_dofs.jl
@@ -65,7 +65,7 @@ end
 @testset "Dofs for quad in 3d (shell)" begin
 
 nodes = [Node{3,Float64}(Vec(0.0,0.0,0.0)), Node{3,Float64}(Vec(1.0,0.0,0.0)), 
-            Node{3,Float64}(Vec(1.0,1.0,0.0)), Node{3,Float64}(Vec(1.0,0.0,0.0)),
+            Node{3,Float64}(Vec(1.0,1.0,0.0)), Node{3,Float64}(Vec(0.0,1.0,0.0)),
             Node{3,Float64}(Vec(2.0,0.0,0.0)), Node{3,Float64}(Vec(2.0,2.0,0.0))]
 
 cells = [Quadrilateral3D((1,2,3,4)), Quadrilateral3D((2,5,6,3))]
@@ -81,12 +81,14 @@ close!(dh)
 @test celldofs(dh,2) == [4,5,6,25,26,27,28,29,30,7,8,9, # u
                          16,17,18,31,32,33,34,35,36,19,20,21]# θ
 
-#3d quad with 2nd order 1d interpolation
+#3d quads with two quadratic interpolations fields
+#Only 1 dim per field for simplicity...
 dh = DofHandler(grid)
-push!(dh, :x, 3, Lagrange{2,RefCube,2}()) 
+push!(dh, :u, 1, Lagrange{2,RefCube,2}())
+push!(dh, :θ, 1, Lagrange{2,RefCube,2}())
 close!(dh)
 
-@test celldofs(dh,1) == [1,2,3,4,5,6,7,8,9,10,11,12, 13,14,15,16,17,18,19,20,21,22,23,24, 25,26,27]
-@test celldofs(dh,2) == [4,5,6,28,29,30,31,32,33,7,8,9, 34,35,36,37,38,39,40,41,42,16,17,18, 43,44,45]
+@test celldofs(dh,1) == collect(1:18)
+@test celldofs(dh,2) == [2, 19, 20, 3, 21, 22, 23, 6, 24, 11, 25, 26, 12, 27, 28, 29, 15, 30]
 
 end
\ No newline at end of file
diff --git a/test/test_grid_dofhandler_vtk.jl b/test/test_grid_dofhandler_vtk.jl
index a8292cca60..acf3a20e9d 100644
--- a/test/test_grid_dofhandler_vtk.jl
+++ b/test/test_grid_dofhandler_vtk.jl
@@ -115,7 +115,6 @@ end
 
     @test collect(getcoordinates(getnodes(grid, 5)).data) ≈ [0.5, 0.5]
 
-
     @test getcells(grid, "cell_set") == [getcells(grid, 1)]
 
     f(x) = Tensor{1,1,Float64}((1 + x[1]^2 + 2x[2]^2, ))
@@ -142,3 +141,26 @@ end
     end
     @test n == div(getncells(grid)*(getncells(grid) + 1), 2)
 end
+
+@testset "Grid sets" begin
+
+    grid = Ferrite.generate_grid(Hexahedron, (1, 1, 1), Vec((0.,0., 0.)), Vec((1.,1.,1.)))
+
+    #Test manual add
+    addcellset!(grid, "cell_set", [1]);
+    addnodeset!(grid, "node_set", [1])
+    addfaceset!(grid, "face_set", [FaceIndex(1,1)])
+    addedgeset!(grid, "edge_set", [EdgeIndex(1,1)])
+    addvertexset!(grid, "vert_set", [VertexIndex(1,1)])
+
+    #Test function add
+    addfaceset!(grid, "left_face", (x)-> x[1] ≈ 0.0)
+    addedgeset!(grid, "left_lower_edge", (x)-> x[1] ≈ 0.0 && x[3] ≈ 0.0)
+    addvertexset!(grid, "left_corner", (x)-> x[1] ≈ 0.0 && x[2] ≈ 0.0 && x[3] ≈ 0.0)
+
+    @test 1 in Ferrite.getnodeset(grid, "node_set")
+    @test FaceIndex(1,5) in getfaceset(grid, "left_face")
+    @test EdgeIndex(1,4) in getedgeset(grid, "left_lower_edge")
+    @test VertexIndex(1,1) in getvertexset(grid, "left_corner")
+
+end