diff --git a/CHANGELOG.md b/CHANGELOG.md
index fa565f8e60..480cbc9b5e 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -6,6 +6,17 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
 ## [Unreleased]
+### Added
+ - Support for classical trilinear and triquadratic wedge elements. ([#581][github-581])
+ - Symmetric quadrature rules up to order 10 for prismatic elements. ([#581][github-581])
+ - Finer granulation of dof distribution, allowing to distribute different amounts of dofs
+   per entity. ([#581][github-581])
+### Fixed
+ - Dof distribution for embedded elements. ([#581][github-581])
+### Other improvements
+ - To clarify the dof management `vertices(ip)`, `edges(ip)` and `faces(ip)` has been
+   deprecated in favor of `vertexdof_indices(ip)`, `edgedof_indices(ip)` and
+   `facedof_indices(ip)`. ([#578][github-578])
 
 ## [0.3.12] - 2023-02-28
 ### Added
@@ -324,6 +335,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 [github-574]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/574
 [github-575]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/575
 [github-578]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/578
+[github-581]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/581
 [github-583]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/583
 [github-588]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/588
 [github-591]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/591
diff --git a/docs/src/devdocs/interpolations.md b/docs/src/devdocs/interpolations.md
index 5dfe0bfe5f..08fbd570d0 100644
--- a/docs/src/devdocs/interpolations.md
+++ b/docs/src/devdocs/interpolations.md
@@ -13,19 +13,27 @@ Ferrite.getrefshape(::Interpolation)
 Ferrite.getorder(::Interpolation)
 Ferrite.value(::Interpolation{dim}, ::Vec{dim,T}) where {dim,T}
 Ferrite.derivative(::Interpolation{dim}, ::Vec{dim}) where {dim}
+Ferrite.boundarydof_indices
 ```
 
 ### Required methods to implement for all subtypes of `Interpolation` to define a new finite element
 
+Depending on the dimension of the reference element the following functions have to be implemented
+
 ```@docs
 Ferrite.value(::Interpolation, ::Int, ::Vec)
-Ferrite.vertices(::Interpolation)
-Ferrite.nvertexdofs(::Interpolation)
-Ferrite.faces(::Interpolation)
-Ferrite.nfacedofs(::Interpolation)
-Ferrite.edges(::Interpolation)
-Ferrite.nedgedofs(::Interpolation)
-Ferrite.ncelldofs(::Interpolation)
+Ferrite.vertexdof_indices(::Interpolation)
+Ferrite.facedof_indices(::Interpolation)
+Ferrite.facedof_interior_indices(::Interpolation)
+Ferrite.edgedof_indices(::Interpolation)
+Ferrite.edgedof_interior_indices(::Interpolation)
+Ferrite.celldof_interior_indices(::Interpolation)
 Ferrite.getnbasefunctions(::Interpolation)
 Ferrite.reference_coordinates(::Interpolation)
 ```
+
+for all entities which exist on that reference element. The dof functions default to having no
+dofs defined on a specific entity. Hence, not overloading of the dof functions will result in an 
+element with zero dofs. Also, it should always be double checked that everything is consistent as 
+specified in the docstring of the corresponding function, as inconsistent implementations can
+lead to bugs which are really difficult to track down.
diff --git a/src/Dofs/ConstraintHandler.jl b/src/Dofs/ConstraintHandler.jl
index c95452b115..635e24e7ca 100644
--- a/src/Dofs/ConstraintHandler.jl
+++ b/src/Dofs/ConstraintHandler.jl
@@ -308,7 +308,7 @@ end
 
 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}
     local_face_dofs, local_face_dofs_offset =
-        _local_face_dofs_for_bc(interpolation, field_dim, dbc.components, offset, boundaryfunction(eltype(bcfaces)))
+        _local_face_dofs_for_bc(interpolation, field_dim, dbc.components, offset, boundarydof_indices(eltype(bcfaces)))
     copy!(dbc.local_face_dofs, local_face_dofs)
     copy!(dbc.local_face_dofs_offset, local_face_dofs_offset)
 
@@ -337,7 +337,7 @@ end
 
 # 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]
-function _local_face_dofs_for_bc(interpolation, field_dim, components, offset, boundaryfunc::F=faces) where F
+function _local_face_dofs_for_bc(interpolation, field_dim, components, offset, boundaryfunc::F=facedof_indices) where F
     @assert issorted(components)
     local_face_dofs = Int[]
     local_face_dofs_offset = Int[1]
@@ -442,24 +442,23 @@ function update!(ch::ConstraintHandler, time::Real=0.0)
     return nothing
 end
 
-# for faces
-function _update!(inhomogeneities::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,
+# for vertices, faces and edges
+function _update!(inhomogeneities::Vector{Float64}, f::Function, boundary_entities::Set{<:BoundaryIndex}, field::Symbol, local_face_dofs::Vector{Int}, local_face_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}
 
     cc = CellCache(dh, UpdateFlags(; nodes=false, coords=true, dofs=true))
-    for (cellidx, faceidx) in faces
+    for (cellidx, entityidx) in boundary_entities
         reinit!(cc, cellidx)
 
-        # no need to reinit!, enough to update current_face since we only need geometric shape functions M
-        facevalues.current_face[] = faceidx
+        # no need to reinit!, enough to update current_entity since we only need geometric shape functions M
+        boundaryvalues.current_entity[] = entityidx
 
         # local dof-range for this face
-        r = local_face_dofs_offset[faceidx]:(local_face_dofs_offset[faceidx+1]-1)
+        r = local_face_dofs_offset[entityidx]:(local_face_dofs_offset[entityidx+1]-1)
         counter = 1
-
-        for location in 1:getnquadpoints(facevalues)
-            x = spatial_coordinate(facevalues, location, cc.coords)
+        for location in 1:getnquadpoints(boundaryvalues)
+            x = spatial_coordinate(boundaryvalues, location, cc.coords)
             bc_value = f(x, time)
             @assert length(bc_value) == length(components)
 
@@ -1302,7 +1301,7 @@ end
 function mirror_local_dofs(local_face_dofs, local_face_dofs_offset, ip::Lagrange{2,<:Union{RefCube,RefTetrahedron}}, n::Int)
     # For 2D we always permute since Ferrite defines dofs counter-clockwise
     ret = collect(1:length(local_face_dofs))
-    for (i, f) in enumerate(faces(ip))
+    for (i, f) in enumerate(facedof_indices(ip))
         this_offset = local_face_dofs_offset[i]
         other_offset = this_offset + n
         for d in 1:n
@@ -1323,7 +1322,7 @@ function mirror_local_dofs(local_face_dofs, local_face_dofs_offset, ip::Lagrange
     ret = collect(1:length(local_face_dofs))
 
     # Mirror by changing from counter-clockwise to clockwise
-    for (i, f) in enumerate(faces(ip))
+    for (i, f) in enumerate(facedof_indices(ip))
         r = local_face_dofs_offset[i]:(local_face_dofs_offset[i+1] - 1)
         # 1. Rotate the corners
         vertex_range = r[1:(N*n)]
diff --git a/src/Dofs/DofHandler.jl b/src/Dofs/DofHandler.jl
index 2cc46d4de7..aafd6f255a 100644
--- a/src/Dofs/DofHandler.jl
+++ b/src/Dofs/DofHandler.jl
@@ -9,13 +9,12 @@ Operates slightly faster than [`MixedDofHandler`](@ref). Supports:
 - `Grid`s with a single concrete cell type.
 - One or several fields on the whole domaine.
 """
-struct DofHandler{dim,T,G<:AbstractGrid{dim}} <: AbstractDofHandler
+struct DofHandler{dim,G<:AbstractGrid{dim}} <: AbstractDofHandler
     field_names::Vector{Symbol}
     field_dims::Vector{Int}
     # 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}} # 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}
@@ -25,7 +24,7 @@ end
 
 function DofHandler(grid::AbstractGrid)
     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))
+    DofHandler(Symbol[], Int[], Interpolation[], Int[], Int[], ScalarWrapper(false), grid, Ferrite.ScalarWrapper(-1))
 end
 
 function Base.show(io::IO, ::MIME"text/plain", dh::DofHandler)
@@ -70,7 +69,6 @@ end
 
 getfieldinterpolation(dh::DofHandler, field_idx::Int) = dh.field_interpolations[field_idx]
 getfielddim(dh::DofHandler, field_idx::Int) = dh.field_dims[field_idx]
-getbcvalue(dh::DofHandler, field_idx::Int) = dh.bc_values[field_idx]
 
 function getfielddim(dh::DofHandler, name::Symbol)
     field_pos = findfirst(i->i == name, getfieldnames(dh))
@@ -108,9 +106,9 @@ end
 
 Add a `dim`-dimensional `Field` called `name` which is approximated by `ip` to `dh`.
 
-The field is added to all cells of the underlying grid. In case no interpolation `ip` is given,
-the default interpolation of the grid's celltype is used. 
-If the grid uses several celltypes, [`add!(dh::MixedDofHandler, fh::FieldHandler)`](@ref) must be used instead.
+The field is added to all cells of the underlying grid. In case no interpolation `ip` is
+given, the default interpolation of the grid's celltype is used. If the grid uses several
+celltypes, [`add!(dh::MixedDofHandler, fh::FieldHandler)`](@ref) must be used instead.
 """
 function add!(dh::DofHandler, name::Symbol, dim::Int, ip::Interpolation=default_interpolation(getcelltype(dh.grid)))
     @assert !isclosed(dh)
@@ -118,10 +116,14 @@ function add!(dh::DofHandler, name::Symbol, dim::Int, ip::Interpolation=default_
     push!(dh.field_names, name)
     push!(dh.field_dims, dim)
     push!(dh.field_interpolations, ip)
-    push!(dh.bc_values, BCValues(ip, default_interpolation(getcelltype(dh.grid))))
     return dh
 end
 
+# Method for supporting dim=1 default
+function add!(dh::DofHandler, name::Symbol, ip::Interpolation=default_interpolation(getcelltype(dh.grid)))
+    return add!(dh, name, 1, ip)
+end
+
 """
     sortedge(edge::Tuple{Int,Int})
 
@@ -203,7 +205,7 @@ function __close!(dh::DofHandler{dim}) where {dim}
     end
 
     # not implemented yet: more than one facedof per face in 3D
-    dim == 3 && @assert(!any(x->x.nfacedofs > 1, interpolation_infos))
+    dim == 3 && @assert(!any(x->any(y->y > 1, x.nfacedofs), interpolation_infos))
 
     nextdof = 1 # next free dof to distribute
     push!(dh.cell_dofs_offset, 1) # dofs for the first cell start at 1
@@ -211,89 +213,112 @@ function __close!(dh::DofHandler{dim}) where {dim}
     # loop over all the cells, and distribute dofs for all the fields
     for (ci, cell) in enumerate(getcells(dh.grid))
         @debug println("cell #$ci")
-        for fi in 1:nfields(dh)
-            interpolation_info = interpolation_infos[fi]
-            @debug println("  field: $(dh.field_names[fi])")
-            if interpolation_info.nvertexdofs > 0
-                for vertex in vertices(cell)
+        for field_idx in 1:nfields(dh)
+            interpolation_info = interpolation_infos[field_idx]
+            @debug println("  field: $(dh.field_names[field_idx])")
+            for (vi, vertex) in enumerate(vertices(cell))
+                if interpolation_info.nvertexdofs[vi] > 0
                     @debug println("    vertex#$vertex")
-                    token = Base.ht_keyindex2!(vertexdicts[fi], vertex)
+                    token = Base.ht_keyindex2!(vertexdicts[field_idx], vertex)
                     if token > 0 # haskey(vertexdicts[fi], vertex) # reuse dofs
-                        reuse_dof = vertexdicts[fi].vals[token] # vertexdicts[fi][vertex]
-                        for d in 1:dh.field_dims[fi]
+                        reuse_dof = vertexdicts[field_idx].vals[token] # vertexdicts[fi][vertex]
+                        for d in 1:dh.field_dims[field_idx]
                             @debug println("      reusing dof #$(reuse_dof + (d-1))")
                             push!(dh.cell_dofs, reuse_dof + (d-1))
                         end
                     else # token <= 0, distribute new dofs
-                        for vertexdof in 1:interpolation_info.nvertexdofs
-                            Base._setindex!(vertexdicts[fi], nextdof, vertex, -token) # vertexdicts[fi][vertex] = nextdof
-                            for d in 1:dh.field_dims[fi]
+                        for vertexdof in 1:interpolation_info.nvertexdofs[vi]
+                            Base._setindex!(vertexdicts[field_idx], nextdof, vertex, -token) # vertexdicts[field_idx][vertex] = nextdof
+                            for d in 1:dh.field_dims[field_idx]
                                 @debug println("      adding dof#$nextdof")
                                 push!(dh.cell_dofs, nextdof)
                                 nextdof += 1
                             end
                         end
                     end
-                end # vertex loop
-            end
+                end
+            end # vertex loop
             if dim == 3 # edges only in 3D
-                if interpolation_info.nedgedofs > 0
-                    for edge in edges(cell)
+                for (ei, edge) in enumerate(edges(cell))
+                    if interpolation_info.dim == 3 && interpolation_info.nedgedofs[ei] > 0
                         sedge, dir = sortedge(edge)
                         @debug println("    edge#$sedge dir: $(dir)")
-                        token = Base.ht_keyindex2!(edgedicts[fi], sedge)
-                        if token > 0 # haskey(edgedicts[fi], sedge), reuse dofs
-                            startdof, olddir = edgedicts[fi].vals[token] # edgedicts[fi][sedge] # first dof for this edge (if dir == true)
-                            for edgedof in (dir == olddir ? (1:interpolation_info.nedgedofs) : (interpolation_info.nedgedofs:-1:1))
-                                for d in 1:dh.field_dims[fi]
-                                    reuse_dof = startdof + (d-1) + (edgedof-1)*dh.field_dims[fi]
+                        token = Base.ht_keyindex2!(edgedicts[field_idx], sedge)
+                        if token > 0 # haskey(edgedicts[field_idx], sedge), reuse dofs
+                            startdof, olddir = edgedicts[field_idx].vals[token] # edgedicts[field_idx][sedge] # first dof for this edge (if dir == true)
+                            for edgedof in (dir == olddir ? (1:interpolation_info.nedgedofs[ei]) : (interpolation_info.nedgedofs[ei]:-1:1))
+                                for d in 1:dh.field_dims[field_idx]
+                                    reuse_dof = startdof + (d-1) + (edgedof-1)*dh.field_dims[field_idx]
                                     @debug println("      reusing dof#$(reuse_dof)")
                                     push!(dh.cell_dofs, reuse_dof)
                                 end
                             end
                         else # token <= 0, distribute new dofs
-                            Base._setindex!(edgedicts[fi], (nextdof, dir), sedge, -token) # edgedicts[fi][sedge] = (nextdof, dir),  store only the first dof for the edge
-                            for edgedof in 1:interpolation_info.nedgedofs
-                                for d in 1:dh.field_dims[fi]
+                            Base._setindex!(edgedicts[field_idx], (nextdof, dir), sedge, -token) # edgedicts[field_idx][sedge] = (nextdof, dir),  store only the first dof for the edge
+                            for edgedof in 1:interpolation_info.nedgedofs[ei]
+                                for d in 1:dh.field_dims[field_idx]
                                     @debug println("      adding dof#$nextdof")
                                     push!(dh.cell_dofs, nextdof)
                                     nextdof += 1
                                 end
                             end
                         end
-                    end # edge loop
-                end
+                    elseif interpolation_info.dim == 2 && interpolation_info.nfacedofs[ei] > 0 # hotfix for current embedded elements
+                        sedge, dir = sortedge(edge)
+                        @debug println("    edge#$sedge dir: $(dir)")
+                        token = Base.ht_keyindex2!(edgedicts[field_idx], sedge)
+                        if token > 0 # haskey(edgedicts[field_idx], sedge), reuse dofs
+                            startdof, olddir = edgedicts[field_idx].vals[token] # edgedicts[field_idx][sedge] # first dof for this edge (if dir == true)
+                            for edgedof in (dir == olddir ? (1:interpolation_info.nfacedofs[ei]) : (interpolation_info.nfacedofs[ei]:-1:1))
+                                for d in 1:dh.field_dims[field_idx]
+                                    reuse_dof = startdof + (d-1) + (edgedof-1)*dh.field_dims[field_idx]
+                                    @debug println("      reusing dof#$(reuse_dof)")
+                                    push!(dh.cell_dofs, reuse_dof)
+                                end
+                            end
+                        else # token <= 0, distribute new dofs
+                            Base._setindex!(edgedicts[field_idx], (nextdof, dir), sedge, -token) # edgedicts[field_idx][sedge] = (nextdof, dir),  store only the first dof for the edge
+                            for edgedof in 1:interpolation_info.nfacedofs[ei]
+                                for d in 1:dh.field_dims[field_idx]
+                                    @debug println("      adding dof#$nextdof")
+                                    push!(dh.cell_dofs, nextdof)
+                                    nextdof += 1
+                                end
+                            end
+                        end
+                    end 
+                end # edge loop
             end
-            if interpolation_info.nfacedofs > 0 && (interpolation_info.dim == dim)
-                for face in faces(cell)
-                    sface = sortface(face) # TODO: faces(cell) may as well just return the sorted list
+            for (fi, face) in enumerate(faces(cell))
+                if interpolation_info.nfacedofs[fi] > 0 && (interpolation_info.dim == dim)
+                    sface = sortface(face)
                     @debug println("    face#$sface")
-                    token = Base.ht_keyindex2!(facedicts[fi], sface)
-                    if token > 0 # haskey(facedicts[fi], sface), reuse dofs
-                        startdof = facedicts[fi].vals[token] # facedicts[fi][sface]
-                        for facedof in interpolation_info.nfacedofs:-1:1 # always reverse (YOLO)
-                            for d in 1:dh.field_dims[fi]
-                                reuse_dof = startdof + (d-1) + (facedof-1)*dh.field_dims[fi]
+                    token = Base.ht_keyindex2!(facedicts[field_idx], sface)
+                    if token > 0 # haskey(facedicts[field_idx], sface), reuse dofs
+                        startdof = facedicts[field_idx].vals[token] # facedicts[field_idx][sface]
+                        for facedof in interpolation_info.nfacedofs[fi]:-1:1 # always reverse (YOLO)
+                            for d in 1:dh.field_dims[field_idx]
+                                reuse_dof = startdof + (d-1) + (facedof-1)*dh.field_dims[field_idx]
                                 @debug println("      reusing dof#$(reuse_dof)")
                                 push!(dh.cell_dofs, reuse_dof)
                             end
                         end
                     else # distribute new dofs
-                        Base._setindex!(facedicts[fi], nextdof, sface, -token)# facedicts[fi][sface] = nextdof,  store the first dof for this face
-                        for facedof in 1:interpolation_info.nfacedofs
-                            for d in 1:dh.field_dims[fi]
+                        Base._setindex!(facedicts[field_idx], nextdof, sface, -token)# facedicts[field_idx][sface] = nextdof,  store the first dof for this face
+                        for facedof in 1:interpolation_info.nfacedofs[fi]
+                            for d in 1:dh.field_dims[field_idx]
                                 @debug println("      adding dof#$nextdof")
                                 push!(dh.cell_dofs, nextdof)
                                 nextdof += 1
                             end
                         end
                     end
-                end # face loop
-            end
+                end
+            end # face loop
             if interpolation_info.ncelldofs > 0 # always distribute new dofs for cell
                 @debug println("    cell#$ci")
                 for celldof in 1:interpolation_info.ncelldofs
-                    for d in 1:dh.field_dims[fi]
+                    for d in 1:dh.field_dims[field_idx]
                         @debug println("      adding dof#$nextdof")
                         push!(dh.cell_dofs, nextdof)
                         nextdof += 1
diff --git a/src/Dofs/MixedDofHandler.jl b/src/Dofs/MixedDofHandler.jl
index 1bc5973286..d6468a0d39 100644
--- a/src/Dofs/MixedDofHandler.jl
+++ b/src/Dofs/MixedDofHandler.jl
@@ -58,7 +58,7 @@ end
 
 function MixedDofHandler(grid::Grid{dim,C,T}) where {dim,C,T}
     ncells = getncells(grid)
-    MixedDofHandler{dim,T,typeof(grid)}(FieldHandler[], CellVector(Int[],zeros(Int,ncells),zeros(Int,ncells)), CellVector(Int[],Int[],Int[]), CellVector(Vec{dim,T}[],Int[],Int[]), Ferrite.ScalarWrapper(false), grid, Ferrite.ScalarWrapper(-1))
+    MixedDofHandler{dim,T,typeof(grid)}(FieldHandler[], CellVector(Int[],zeros(Int,ncells),zeros(Int,ncells)), CellVector(Int[],Int[],Int[]), CellVector(Vec{dim,T}[],Int[],Int[]), ScalarWrapper(false), grid, ScalarWrapper(-1))
 end
 
 getfieldnames(fh::FieldHandler) = [field.name for field in fh.fields]
@@ -205,9 +205,8 @@ function close!(dh::MixedDofHandler)
 end
 
 function __close!(dh::MixedDofHandler{dim}) where {dim}
-
-    @assert !Ferrite.isclosed(dh)
-    field_names = Ferrite.getfieldnames(dh)  # all the fields in the problem
+    @assert !isclosed(dh)
+    field_names = getfieldnames(dh)  # all the fields in the problem
     numfields =  length(field_names)
 
     # Create dicts that store created dofs
@@ -228,9 +227,9 @@ function __close!(dh::MixedDofHandler{dim}) where {dim}
             dh,
             cellnumbers,
             field_names,
-            Ferrite.getfieldnames(fh),
-            Ferrite.getfielddims(fh),
-            Ferrite.getfieldinterpolations(fh),
+            getfieldnames(fh),
+            getfielddims(fh),
+            getfieldinterpolations(fh),
             nextdof,
             vertexdicts,
             edgedicts,
@@ -259,15 +258,12 @@ function __close!(dh::MixedDofHandler{dim}) where {dim}
 end
 
 function _close!(dh::MixedDofHandler{dim}, cellnumbers, global_field_names, field_names, field_dims, field_interpolations, nextdof, vertexdicts, edgedicts, facedicts, celldicts) where {dim}
-
-    ip_infos = Ferrite.InterpolationInfo[]
+    ip_infos = InterpolationInfo[]
     for interpolation in field_interpolations
-        ip_info = Ferrite.InterpolationInfo(interpolation)
+        ip_info = InterpolationInfo(interpolation)
         # these are not implemented yet (or have not been tested)
-        @assert(ip_info.nvertexdofs <= 1)
-        @assert(ip_info.nedgedofs <= 1)
-        @assert(ip_info.nfacedofs <= 1)
-        @assert(ip_info.ncelldofs <= 1)  # not tested but probably works
+        @assert(all(ip_info.nedgedofs .<= 1))
+        @assert(all(ip_info.nfacedofs .<= 1))
         push!(ip_infos, ip_info)
     end
 
@@ -285,22 +281,25 @@ function _close!(dh::MixedDofHandler{dim}, cellnumbers, global_field_names, fiel
             @debug "\tfield: $(field_name)"
             ip_info = ip_infos[local_num]
 
-            if ip_info.nvertexdofs > 0
-                nextdof = add_vertex_dofs(cell_dofs, cell, vertexdicts[fi], field_dims[local_num], ip_info.nvertexdofs, nextdof)
-            end
+            # We first distribute the vertex dofs
+            nextdof = add_vertex_dofs(cell_dofs, cell, vertexdicts[fi], field_dims[local_num], ip_info.nvertexdofs, nextdof)
 
-            if ip_info.nedgedofs > 0 && dim == 3 #Edges only in 3d
-                nextdof = add_edge_dofs(cell_dofs, cell, edgedicts[fi], field_dims[local_num], ip_info.nedgedofs, nextdof)
+            # Then the edge dofs
+            if dim == 3
+                if ip_info.dim == 3 # Regular 3D element
+                    nextdof = add_edge_dofs(cell_dofs, cell, edgedicts[fi], field_dims[local_num], ip_info.nedgedofs, nextdof)
+                elseif ip_info.dim == 2 # 2D embedded element in 3D
+                    nextdof = add_edge_dofs(cell_dofs, cell, edgedicts[fi], field_dims[local_num], ip_info.nfacedofs, nextdof)
+                end
             end
 
-            if ip_info.nfacedofs > 0 && (ip_info.dim == dim)
+            # Then the face dofs
+            if ip_info.dim == dim # Regular 3D element
                 nextdof = add_face_dofs(cell_dofs, cell, facedicts[fi], field_dims[local_num], ip_info.nfacedofs, nextdof)
             end
 
-            if ip_info.ncelldofs > 0
-                nextdof = add_cell_dofs(cell_dofs, ci, celldicts[fi], field_dims[local_num], ip_info.ncelldofs, nextdof)
-            end
-
+            # And finally the celldofs
+            nextdof = add_cell_dofs(cell_dofs, ci, celldicts[fi], field_dims[local_num], ip_info.ncelldofs, nextdof)
         end
         # after done creating dofs for the cell, push them to the global list
         append!(dh.cell_dofs.values, cell_dofs)
@@ -316,13 +315,11 @@ Returns the next global dof number and an array of dofs.
 If dofs have already been created for the object (vertex, face) then simply return those, otherwise create new dofs.
 """
 function get_or_create_dofs!(nextdof, field_dim; dict, key)
-
     token = Base.ht_keyindex2!(dict, key)
     if token > 0  # vertex, face etc. visited before
         # reuse stored dofs (TODO unless field is discontinuous)
         @debug "\t\tkey: $key dofs: $(dict[key])  (reused dofs)"
         return nextdof, dict[key]
-
     else  # create new dofs
         dofs = nextdof : (nextdof + field_dim-1)
         @debug "\t\tkey: $key dofs: $dofs"
@@ -333,33 +330,39 @@ function get_or_create_dofs!(nextdof, field_dim; dict, key)
 end
 
 function add_vertex_dofs(cell_dofs, cell, vertexdict, field_dim, nvertexdofs, nextdof)
-    for vertex in Ferrite.vertices(cell)
-        @debug "\tvertex #$vertex"
-        nextdof, dofs = get_or_create_dofs!(nextdof, field_dim, dict=vertexdict, key=vertex)
-        append!(cell_dofs, dofs)
+    for (vi, vertex) in enumerate(vertices(cell))
+        if nvertexdofs[vi] > 0
+            @debug "\tvertex #$vertex"
+            nextdof, dofs = get_or_create_dofs!(nextdof, field_dim, dict=vertexdict, key=vertex)
+            append!(cell_dofs, dofs)
+        end
     end
     return nextdof
 end
 
 function add_face_dofs(cell_dofs, cell, facedict, field_dim, nfacedofs, nextdof)
-    @assert nfacedofs == 1 "Currently only supports interpolations with nfacedofs = 1"
-
-    for face in Ferrite.faces(cell)
-        sface = Ferrite.sortface(face)
-        @debug "\tface #$sface"
-        nextdof, dofs = get_or_create_dofs!(nextdof, field_dim, dict=facedict, key=sface)
-        append!(cell_dofs, dofs)
+    @debug @assert all(nfacedofs .<= 1) "Currently only supports interpolations with less that 2 dofs per face"
+
+    for (fi,face) in enumerate(faces(cell))
+        if nfacedofs[fi] > 0
+            sface = sortface(face)
+            @debug "\tface #$sface"
+            nextdof, dofs = get_or_create_dofs!(nextdof, field_dim, dict=facedict, key=sface)
+            # TODO permutate dofs according to face orientation
+            append!(cell_dofs, dofs)
+        end
     end
     return nextdof
 end
 
 function add_edge_dofs(cell_dofs, cell, edgedict, field_dim, nedgedofs, nextdof)
-    @assert nedgedofs == 1 "Currently only supports interpolations with nedgedofs = 1"
-    for edge in Ferrite.edges(cell)
-        sedge, dir = Ferrite.sortedge(edge)
-        @debug "\tedge #$sedge"
-        nextdof, dofs = get_or_create_dofs!(nextdof, field_dim, dict=edgedict, key=sedge)
-        append!(cell_dofs, dofs)
+    for (ei,edge) in enumerate(edges(cell))
+        if nedgedofs[ei] > 0
+            sedge, dir = sortedge(edge)
+            @debug "\tedge #$sedge"
+            nextdof, dofs = get_or_create_dofs!(nextdof, field_dim, dict=edgedict, key=sedge)
+            append!(cell_dofs, dofs)
+        end
     end
     return nextdof
 end
@@ -389,9 +392,9 @@ function field_offset(fh::FieldHandler, field_name::Symbol)
     return offset
 end
 
-function Ferrite.dof_range(fh::FieldHandler, field_name::Symbol)
-    f = Ferrite.find_field(fh, field_name)
-    offset = Ferrite.field_offset(fh, field_name)
+function dof_range(fh::FieldHandler, field_name::Symbol)
+    f = find_field(fh, field_name)
+    offset = field_offset(fh, field_name)
     field_interpolation = fh.fields[f].interpolation
     field_dim = fh.fields[f].dim
     n_field_dofs = getnbasefunctions(field_interpolation)::Int * field_dim
@@ -406,9 +409,8 @@ getfielddim(fh::FieldHandler, field_idx::Int) = fh.fields[field_idx].dim
 getfielddim(dh::MixedDofHandler, field_idx::Int) = dh.fieldhandlers[1].fields[field_idx].dim
 
 function reshape_to_nodes(dh::MixedDofHandler, u::Vector{T}, fieldname::Symbol) where T
-
     # make sure the field exists
-    fieldname ∈ Ferrite.getfieldnames(dh) || error("Field $fieldname not found.")
+    fieldname ∈ getfieldnames(dh) || error("Field $fieldname not found.")
 
     field_dim = getfielddim(dh, fieldname)
     space_dim = field_dim == 2 ? 3 : field_dim
diff --git a/src/Export/VTK.jl b/src/Export/VTK.jl
index 1977d8a692..2d54ad99c4 100644
--- a/src/Export/VTK.jl
+++ b/src/Export/VTK.jl
@@ -14,6 +14,7 @@ cell_to_vtkcell(::Type{Hexahedron}) = VTKCellTypes.VTK_HEXAHEDRON
 cell_to_vtkcell(::Type{Cell{3,20,6}}) = VTKCellTypes.VTK_QUADRATIC_HEXAHEDRON
 cell_to_vtkcell(::Type{Tetrahedron}) = VTKCellTypes.VTK_TETRA
 cell_to_vtkcell(::Type{QuadraticTetrahedron}) = VTKCellTypes.VTK_QUADRATIC_TETRA
+cell_to_vtkcell(::Type{Wedge}) = VTKCellTypes.VTK_WEDGE
 
 nodes_to_vtkorder(cell::AbstractCell) = collect(cell.nodes)
 
diff --git a/src/FEValues/face_values.jl b/src/FEValues/face_values.jl
index 25e941b18d..6b65180627 100644
--- a/src/FEValues/face_values.jl
+++ b/src/FEValues/face_values.jl
@@ -232,48 +232,50 @@ for each dof-position determined by the `func_interpol`. Used mainly by the `Con
 """
 struct BCValues{T}
     M::Array{T,3}
-    current_face::ScalarWrapper{Int}
+    nqp::Array{Int}
+    current_entity::ScalarWrapper{Int}
 end
 
 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}, boundary_type::Type{<:BoundaryIndex} = Ferrite.FaceIndex) where {T,dim,refshape}
-    # set up quadrature rules for each face with dof-positions
+    # set up quadrature rules for each boundary entity 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)
+    for boundarydofs in boundarydof_indices(boundary_type)(func_interpol)
         dofcoords = Vec{dim,T}[]
-        for facedof in face
-            push!(dofcoords, interpolation_coords[facedof])
+        for boundarydof in boundarydofs
+            push!(dofcoords, interpolation_coords[boundarydof])
         end
         qrf = QuadratureRule{dim,refshape,T}(fill(T(NaN), length(dofcoords)), dofcoords) # weights will not be used
         push!(qrs, qrf)
     end
 
-    n_faces = length(qrs)
-    n_qpoints = n_faces == 0 ? 0 : length(getweights(qrs[1])) # assume same in all
+    n_boundary_entities = length(qrs)
+    n_qpoints = n_boundary_entities == 0 ? 0 : maximum(qr->length(getweights(qr)), qrs) # Bound number of qps correctly.
     n_geom_basefuncs = getnbasefunctions(geom_interpol)
-    M =    fill(zero(T)           * T(NaN), n_geom_basefuncs, n_qpoints, n_faces)
+    M   = fill(zero(T) * T(NaN), n_geom_basefuncs, n_qpoints, n_boundary_entities)
+    nqp = zeros(Int,n_boundary_entities)
 
-    for face in 1:n_faces, (qp, ξ) in enumerate(qrs[face].points)
-        for i in 1:n_geom_basefuncs
-            M[i, qp, face] = value(geom_interpol, i, ξ)
+    for n_boundary_entity in 1:n_boundary_entities
+        for (qp, ξ) in enumerate(qrs[n_boundary_entity].points), i in 1:n_geom_basefuncs
+            M[i, qp, n_boundary_entity] = value(geom_interpol, i, ξ)
         end
+        nqp[n_boundary_entity] = length(qrs[n_boundary_entity].points)
     end
 
-    BCValues{T}(M, ScalarWrapper(0))
+    BCValues{T}(M, nqp, ScalarWrapper(0))
 end
 
-getnquadpoints(bcv::BCValues) = size(bcv.M, 2)
+getnquadpoints(bcv::BCValues) = bcv.nqp[bcv.current_entity.x]
 function spatial_coordinate(bcv::BCValues, q_point::Int, xh::AbstractVector{Vec{dim,T}}) where {dim,T}
     n_base_funcs = size(bcv.M, 1)
     length(xh) == n_base_funcs || throw_incompatible_coord_length(length(xh), n_base_funcs)
     x = zero(Vec{dim,T})
-    face = bcv.current_face[]
+    face = bcv.current_entity[]
     @inbounds for i in 1:n_base_funcs
         x += bcv.M[i,q_point,face] * xh[i] # geometric_value(fe_v, q_point, i) * xh[i]
     end
diff --git a/src/Ferrite.jl b/src/Ferrite.jl
index 63f26d8f46..3ee75b6e95 100644
--- a/src/Ferrite.jl
+++ b/src/Ferrite.jl
@@ -20,6 +20,7 @@ abstract type AbstractRefShape end
 
 struct RefTetrahedron <: AbstractRefShape end
 struct RefCube <: AbstractRefShape end
+struct RefPrism <: AbstractRefShape end
 
 """
 Abstract type which has `CellValues` and `FaceValues` as subtypes
@@ -33,6 +34,34 @@ Abstract type which is used as identifier for faces, edges and verices
 """
 abstract type BoundaryIndex end
 
+"""
+A `CellIndex` wraps an Int and corresponds to a cell with that number in the mesh
+"""
+struct CellIndex
+    idx::Int
+end
+
+"""
+A `FaceIndex` wraps an (Int, Int) and defines a local face by pointing to a (cell, face).
+"""
+struct FaceIndex <: BoundaryIndex
+    idx::Tuple{Int,Int} # cell and side
+end
+
+"""
+A `EdgeIndex` wraps an (Int, Int) and defines a local edge 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 local vertex by pointing to a (cell, vert).
+"""
+struct VertexIndex <: BoundaryIndex
+    idx::Tuple{Int,Int} # cell and side
+end
+
 include("utils.jl")
 
 # Matrix/Vector utilities
diff --git a/src/Grid/grid.jl b/src/Grid/grid.jl
index b1d24a395f..ee862abab1 100644
--- a/src/Grid/grid.jl
+++ b/src/Grid/grid.jl
@@ -100,36 +100,10 @@ const implemented_celltypes = (
     (const QuadraticTetrahedron = Cell{3,10,4}),
     
     (const Hexahedron = Cell{3,8,6}),
-    (Cell{2,20,6})
-)
-
-"""
-A `CellIndex` wraps an Int and corresponds to a cell with that number in the mesh
-"""
-struct CellIndex
-    idx::Int
-end
-
-"""
-A `FaceIndex` wraps an (Int, Int) and defines a local face by pointing to a (cell, face).
-"""
-struct FaceIndex <: BoundaryIndex
-    idx::Tuple{Int,Int} # cell and side
-end
+    (Cell{2,20,6}),
 
-"""
-A `EdgeIndex` wraps an (Int, Int) and defines a local edge 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 local vertex by pointing to a (cell, vert).
-"""
-struct VertexIndex <: BoundaryIndex
-    idx::Tuple{Int,Int} # cell and side
-end
+    (const Wedge = Cell{3,6,5})
+)
 
 struct EntityNeighborhood{T<:Union{BoundaryIndex,CellIndex}}
     neighbor_info::Vector{T}
@@ -863,6 +837,11 @@ faces(c::Union{Hexahedron,Cell{3,20,6}}) = ((c.nodes[1],c.nodes[4],c.nodes[3],c.
 edges(c::Union{Quadrilateral3D}) = ((c.nodes[1],c.nodes[2]), (c.nodes[2],c.nodes[3]), (c.nodes[3],c.nodes[4]), (c.nodes[4],c.nodes[1]))
 faces(c::Union{Quadrilateral3D}) = ((c.nodes[1],c.nodes[2],c.nodes[3],c.nodes[4]),)
 
+vertices(c::Wedge) = (c.nodes[1], c.nodes[2], c.nodes[3], c.nodes[4], c.nodes[5], c.nodes[6])
+edges(c::Wedge) = ((c.nodes[2],c.nodes[1]), (c.nodes[1],c.nodes[3]), (c.nodes[1],c.nodes[4]), (c.nodes[3],c.nodes[2]), (c.nodes[2],c.nodes[5]), (c.nodes[3],c.nodes[6]), (c.nodes[4],c.nodes[5]), (c.nodes[4],c.nodes[6]), (c.nodes[6],c.nodes[5]))
+faces(c::Wedge) = ((c.nodes[1],c.nodes[3],c.nodes[2]), (c.nodes[1],c.nodes[2],c.nodes[5],c.nodes[4]), (c.nodes[3],c.nodes[1],c.nodes[4],c.nodes[6]), (c.nodes[2],c.nodes[3],c.nodes[6],c.nodes[5]), (c.nodes[4],c.nodes[5],c.nodes[6]))
+
+# random stuff
 default_interpolation(::Union{Type{Line},Type{Line2D},Type{Line3D}}) = Lagrange{1,RefCube,1}()
 default_interpolation(::Type{QuadraticLine}) = Lagrange{1,RefCube,2}()
 default_interpolation(::Type{Triangle}) = Lagrange{2,RefTetrahedron,1}()
@@ -873,6 +852,7 @@ 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{Cell{3,20,6}}) = Serendipity{3,RefCube,2}()
+default_interpolation(::Type{Wedge}) = Lagrange{3,RefPrism,1}()
 
 """
     boundaryfunction(::Type{<:BoundaryIndex})
diff --git a/src/Grid/grid_generators.jl b/src/Grid/grid_generators.jl
index 7b2b0e424a..d5ca1c8f36 100644
--- a/src/Grid/grid_generators.jl
+++ b/src/Grid/grid_generators.jl
@@ -239,6 +239,55 @@ function generate_grid(::Type{Hexahedron}, nel::NTuple{3,Int}, left::Vec{3,T}=Ve
     return Grid(cells, nodes, facesets=facesets, boundary_matrix=boundary_matrix)
 end 
 
+# Wedge
+function generate_grid(::Type{Wedge}, nel::NTuple{3,Int}, left::Vec{3,T}=Vec{3}((-1.0,-1.0,-1.0)), right::Vec{3,T}=Vec{3}((1.0,1.0,1.0))) where {T}
+    nel_x = nel[1]; nel_y = nel[2]; nel_z = nel[3]; nel_tot = nel_x*nel_y*nel_z
+    n_nodes_x = nel_x + 1; n_nodes_y = nel_y + 1; n_nodes_z = nel_z + 1
+    n_nodes = n_nodes_x * n_nodes_y * n_nodes_z
+
+    # Generate nodes
+    coords_x = range(left[1], stop=right[1], length=n_nodes_x)
+    coords_y = range(left[2], stop=right[2], length=n_nodes_y)
+    coords_z = range(left[3], stop=right[3], length=n_nodes_z)
+    nodes = Node{3,T}[]
+    for k in 1:n_nodes_z, j in 1:n_nodes_y, i in 1:n_nodes_x
+        push!(nodes, Node((coords_x[i], coords_y[j], coords_z[k])))
+    end
+
+    # Generate cells
+    node_array = reshape(collect(1:n_nodes), (n_nodes_x, n_nodes_y, n_nodes_z))
+    cells = Wedge[]
+    for k in 1:nel_z, j in 1:nel_y, i in 1:nel_x
+        push!(cells, Wedge((node_array[i,j,k], node_array[i+1,j,k], node_array[i,j+1,k], 
+                            node_array[i,j,k+1], node_array[i+1,j,k+1], node_array[i,j+1,k+1]))) # ◺
+        push!(cells, Wedge((node_array[i+1,j,k], node_array[i+1,j+1,k], node_array[i,j+1,k], 
+                            node_array[i+1,j,k+1], node_array[i+1,j+1,k+1], node_array[i,j+1,k+1]))) # ◹
+    end
+
+    # Order the cells as c_nxyz[2, x, y, z] such that we can look up boundary cells
+    c_nxyz = reshape(1:length(cells), (2, nel...))
+
+    @views le = map(x -> FaceIndex(x,3), c_nxyz[1,   1, :, :][:])
+    @views ri = map(x -> FaceIndex(x,2), c_nxyz[2, end, :, :][:])
+    @views fr = map(x -> FaceIndex(x,2), c_nxyz[1, :, 1, :][:])  
+    @views ba = map(x -> FaceIndex(x,4), c_nxyz[2, :, end, :][:])
+    @views bo = [map(x -> FaceIndex(x,1), c_nxyz[1, :, :, 1][:])   ; map(x -> FaceIndex(x,1), c_nxyz[2, :, :, 1][:])]
+    @views to = [map(x -> FaceIndex(x,5), c_nxyz[1, :, :, end][:]) ; map(x -> FaceIndex(x,5), c_nxyz[2, :, :, end][:])]
+
+    boundary_matrix = boundaries_to_sparse([le; ri; bo; to; fr; ba])
+
+    facesets = Dict(
+        "left" => Set(le),
+        "right" => Set(ri),
+        "front" => Set(fr),
+        "back" => Set(ba),
+        "bottom" => Set(bo),
+        "top" => Set(to),
+    )
+
+    return Grid(cells, nodes, facesets=facesets, boundary_matrix=boundary_matrix)
+end 
+
 function Ferrite.generate_grid(::Type{Cell{3,20,6}}, nel::NTuple{3,Int}, left::Vec{3,T}=Vec{3}((-1.0,-1.0,-1.0)), right::Vec{3,T}=Vec{3}((1.0,1.0,1.0))) where {T}
     nel_x = nel[1]; nel_y = nel[2]; nel_z = nel[3]; nel_tot = nel_x*nel_y*nel_z
     nnode_x = 2nel_x + 1; nnode_y = 2nel_y + 1; nnode_z = 2nel_z + 1 #Note: not the actually number of nodes in x/y/z, just a temporary variables
diff --git a/src/Quadrature/gaussquad_prism_table.jl b/src/Quadrature/gaussquad_prism_table.jl
new file mode 100644
index 0000000000..8fe0c88724
--- /dev/null
+++ b/src/Quadrature/gaussquad_prism_table.jl
@@ -0,0 +1,350 @@
+# Symmetric quadrature rules takes from
+#   Witherden, Freddie D., and Peter E. Vincent. "On the identification of 
+#   symmetric quadrature rules for finite element methods." Computers & 
+#   Mathematics with Applications 69.10 (2015): 1232-1241.
+# Note that the original rule is defined on [-1,1]^3 while our reference prism is defined on [0,1]^3, hence we transform in the end.
+function _get_gauss_prismdata_polyquad(n::Int)
+    if n == 1
+        return reshape([1/3 1/3 1/2 1/2], (1,4))
+    elseif n == 2
+        xw = [
+            -0.33333333333333333333333333333333333333  -0.33333333333333333333333333333333333333  -0.99999999999999985119303811076522004183   0.66666666666666686507594918564641756458
+            -0.33333333333333333333333333333333333333  -0.33333333333333333333333333333333333333   0.99999999999999985119303811076522004183   0.66666666666666686507594918564641756458
+            -0.74158162379719638007464124598498266439   0.48316324759439276014928249196996532878                                          0   0.88888888888888875661603387623572162361
+             0.48316324759439276014928249196996532878  -0.74158162379719638007464124598498266439                                          0   0.88888888888888875661603387623572162361
+            -0.74158162379719638007464124598498266439  -0.74158162379719638007464124598498266439                                          0   0.88888888888888875661603387623572162361
+        ]
+    elseif n == 3
+        xw = [
+            -0.33333333333333333333333333333333333333  -0.33333333333333333333333333333333333333  -0.86066920446477793249332885107223920359   0.89998695257960196814387903509173005101
+            -0.33333333333333333333333333333333333333  -0.33333333333333333333333333333333333333   0.86066920446477793249332885107223920359   0.89998695257960196814387903509173005101
+            -0.52223312808684327863699399822724711768   0.52222829574978345631488370728601281383                                          0     0.366671015806799343952040321636089983
+             0.52222829574978345631488370728601281383  -0.52223312808684327863699399822724711768                                          0     0.366671015806799343952040321636089983
+            -0.99999516766294017767788970905876569615   0.52222829574978345631488370728601281383                                          0     0.366671015806799343952040321636089983
+             0.52222829574978345631488370728601281383  -0.99999516766294017767788970905876569615                                          0     0.366671015806799343952040321636089983
+            -0.99999516766294017767788970905876569615  -0.52223312808684327863699399822724711768                                          0     0.366671015806799343952040321636089983
+            -0.52223312808684327863699399822724711768  -0.99999516766294017767788970905876569615                                          0     0.366671015806799343952040321636089983
+        ]
+    elseif n == 4
+        xw = [
+            -0.33333333333333333333333333333333333333   -0.33333333333333333333333333333333333333   -0.86686197400903479841583223937151297546    0.43164789926214201915560571539043719684
+            -0.33333333333333333333333333333333333333   -0.33333333333333333333333333333333333333    0.86686197400903479841583223937151297546    0.43164789926214201915560571539043719684
+           -0.062688380276009575455822705856210492672   -0.87462323944798084908835458828757901466                                           0    0.54565845042191057226153838655080393041
+            -0.87462323944798084908835458828757901466  -0.062688380276009575455822705856210492672                                           0    0.54565845042191057226153838655080393041
+           -0.062688380276009575455822705856210492672  -0.062688380276009575455822705856210492672                                           0    0.54565845042191057226153838655080393041
+            -0.79851918840217872745314255427215462609    0.59703837680435745490628510854430925218   -0.67563982368225971739661759355428029031    0.24995480836833070748402890159445230251
+             0.59703837680435745490628510854430925218   -0.79851918840217872745314255427215462609   -0.67563982368225971739661759355428029031    0.24995480836833070748402890159445230251
+            -0.79851918840217872745314255427215462609   -0.79851918840217872745314255427215462609   -0.67563982368225971739661759355428029031    0.24995480836833070748402890159445230251
+            -0.79851918840217872745314255427215462609    0.59703837680435745490628510854430925218    0.67563982368225971739661759355428029031    0.24995480836833070748402890159445230251
+             0.59703837680435745490628510854430925218   -0.79851918840217872745314255427215462609    0.67563982368225971739661759355428029031    0.24995480836833070748402890159445230251
+            -0.79851918840217872745314255427215462609   -0.79851918840217872745314255427215462609    0.67563982368225971739661759355428029031    0.24995480836833070748402890159445230251
+        ]
+    elseif n == 5
+        xw = [
+            -0.33333333333333333333333333333333333333   -0.33333333333333333333333333333333333333                                           0    0.71145555593148684996197343022099556338
+            -0.025400070899508667767082370124706051859   -0.94919985820098266446583525975058789628                                           0    0.22471006722826685202639587426081419504
+             -0.94919985820098266446583525975058789628  -0.025400070899508667767082370124706051859                                           0    0.22471006722826685202639587426081419504
+            -0.025400070899508667767082370124706051859  -0.025400070899508667767082370124706051859                                           0    0.22471006722826685202639587426081419504
+             -0.10880379065925599247444934872713650457   -0.78239241868148801505110130254572699087    -0.8710029348654440527292590830738210619    0.18566142131615813987257487985503458339
+             -0.78239241868148801505110130254572699087   -0.10880379065925599247444934872713650457    -0.8710029348654440527292590830738210619    0.18566142131615813987257487985503458339
+             -0.10880379065925599247444934872713650457   -0.10880379065925599247444934872713650457    -0.8710029348654440527292590830738210619    0.18566142131615813987257487985503458339
+             -0.10880379065925599247444934872713650457   -0.78239241868148801505110130254572699087     0.8710029348654440527292590830738210619    0.18566142131615813987257487985503458339
+             -0.78239241868148801505110130254572699087   -0.10880379065925599247444934872713650457     0.8710029348654440527292590830738210619    0.18566142131615813987257487985503458339
+             -0.10880379065925599247444934872713650457   -0.10880379065925599247444934872713650457     0.8710029348654440527292590830738210619    0.18566142131615813987257487985503458339
+             -0.79828210803458293816870337538129563058    0.59656421606916587633740675076259126116   -0.57042698070515927206196786842334325646    0.25007428574779395912056494464439239186
+              0.59656421606916587633740675076259126116   -0.79828210803458293816870337538129563058   -0.57042698070515927206196786842334325646    0.25007428574779395912056494464439239186
+             -0.79828210803458293816870337538129563058   -0.79828210803458293816870337538129563058   -0.57042698070515927206196786842334325646    0.25007428574779395912056494464439239186
+             -0.79828210803458293816870337538129563058    0.59656421606916587633740675076259126116    0.57042698070515927206196786842334325646    0.25007428574779395912056494464439239186
+              0.59656421606916587633740675076259126116   -0.79828210803458293816870337538129563058    0.57042698070515927206196786842334325646    0.25007428574779395912056494464439239186
+             -0.79828210803458293816870337538129563058   -0.79828210803458293816870337538129563058    0.57042698070515927206196786842334325646    0.25007428574779395912056494464439239186
+        ]
+    elseif n == 6
+        wx = [
+            -0.33333333333333333333333333333333333333   -0.33333333333333333333333333333333333333   -0.52662270480497475869798007748379244322    0.22269313409222502250038629788629976598
+            -0.33333333333333333333333333333333333333   -0.33333333333333333333333333333333333333    0.52662270480497475869798007748379244322    0.22269313409222502250038629788629976598
+            -0.33333333333333333333333333333333333333   -0.33333333333333333333333333333333333333   -0.99083630081924474869286718300205167763    0.13839610937412451172575590081432159364
+            -0.33333333333333333333333333333333333333   -0.33333333333333333333333333333333333333    0.99083630081924474869286718300205167763    0.13839610937412451172575590081432159364
+            -0.66959346994381506531292392577889361708    0.33918693988763013062584785155778723416                                           0    0.31182801481333532900066338577135735737
+             0.33918693988763013062584785155778723416   -0.66959346994381506531292392577889361708                                           0    0.31182801481333532900066338577135735737
+            -0.66959346994381506531292392577889361708   -0.66959346994381506531292392577889361708                                           0    0.31182801481333532900066338577135735737
+            -0.95905346077548616195214584157965381704    0.91810692155097232390429168315930763408                                           0   0.050146535238696554002545078611329875622
+             0.91810692155097232390429168315930763408   -0.95905346077548616195214584157965381704                                           0   0.050146535238696554002545078611329875622
+            -0.95905346077548616195214584157965381704   -0.95905346077548616195214584157965381704                                           0   0.050146535238696554002545078611329875622
+           -0.067776668819680981202280000599855057023   -0.86444666236063803759543999880028988595   -0.48167413488438751267079566597126764678    0.19605056020608216687021916179149498019
+            -0.86444666236063803759543999880028988595  -0.067776668819680981202280000599855057023   -0.48167413488438751267079566597126764678    0.19605056020608216687021916179149498019
+           -0.067776668819680981202280000599855057023  -0.067776668819680981202280000599855057023   -0.48167413488438751267079566597126764678    0.19605056020608216687021916179149498019
+           -0.067776668819680981202280000599855057023   -0.86444666236063803759543999880028988595    0.48167413488438751267079566597126764678    0.19605056020608216687021916179149498019
+            -0.86444666236063803759543999880028988595  -0.067776668819680981202280000599855057023    0.48167413488438751267079566597126764678    0.19605056020608216687021916179149498019
+           -0.067776668819680981202280000599855057023  -0.067776668819680981202280000599855057023    0.48167413488438751267079566597126764678    0.19605056020608216687021916179149498019
+            -0.59485220654953844609181984939099585377    0.52549348166744454603903998678937679713   -0.80928579325583275231645432806703865787   0.084632875139559356776397936558477141722
+             0.52549348166744454603903998678937679713   -0.59485220654953844609181984939099585377   -0.80928579325583275231645432806703865787   0.084632875139559356776397936558477141722
+            -0.93064127511790609994722013739838094335    0.52549348166744454603903998678937679713   -0.80928579325583275231645432806703865787   0.084632875139559356776397936558477141722
+             0.52549348166744454603903998678937679713   -0.93064127511790609994722013739838094335   -0.80928579325583275231645432806703865787   0.084632875139559356776397936558477141722
+            -0.93064127511790609994722013739838094335   -0.59485220654953844609181984939099585377   -0.80928579325583275231645432806703865787   0.084632875139559356776397936558477141722
+            -0.59485220654953844609181984939099585377   -0.93064127511790609994722013739838094335   -0.80928579325583275231645432806703865787   0.084632875139559356776397936558477141722
+            -0.59485220654953844609181984939099585377    0.52549348166744454603903998678937679713    0.80928579325583275231645432806703865787   0.084632875139559356776397936558477141722
+             0.52549348166744454603903998678937679713   -0.59485220654953844609181984939099585377    0.80928579325583275231645432806703865787   0.084632875139559356776397936558477141722
+            -0.93064127511790609994722013739838094335    0.52549348166744454603903998678937679713    0.80928579325583275231645432806703865787   0.084632875139559356776397936558477141722
+             0.52549348166744454603903998678937679713   -0.93064127511790609994722013739838094335    0.80928579325583275231645432806703865787   0.084632875139559356776397936558477141722
+            -0.93064127511790609994722013739838094335   -0.59485220654953844609181984939099585377    0.80928579325583275231645432806703865787   0.084632875139559356776397936558477141722
+            -0.59485220654953844609181984939099585377   -0.93064127511790609994722013739838094335    0.80928579325583275231645432806703865787   0.084632875139559356776397936558477141722
+        ]
+    elseif n == 7
+        wx = [
+            -0.33333333333333333333333333333333333333   -0.33333333333333333333333333333333333333   -0.98022806959089160171504882914128008603    0.11378272445075715563392222208600913026
+            -0.33333333333333333333333333333333333333   -0.33333333333333333333333333333333333333    0.98022806959089160171504882914128008603    0.11378272445075715563392222208600913026
+            -0.98332480906795705560590831565479581293    0.96664961813591411121181663130959162585                                           0   0.024472968692380158977782543188506640362
+             0.96664961813591411121181663130959162585   -0.98332480906795705560590831565479581293                                           0   0.024472968692380158977782543188506640362
+            -0.98332480906795705560590831565479581293   -0.98332480906795705560590831565479581293                                           0   0.024472968692380158977782543188506640362
+           -0.036956049341726798726440754485546402356   -0.92608790131654640254711849102890719529   -0.84138735420652599295096060100359489758   0.086203456344905850875965416526967283444
+            -0.92608790131654640254711849102890719529  -0.036956049341726798726440754485546402356   -0.84138735420652599295096060100359489758   0.086203456344905850875965416526967283444
+           -0.036956049341726798726440754485546402356  -0.036956049341726798726440754485546402356   -0.84138735420652599295096060100359489758   0.086203456344905850875965416526967283444
+           -0.036956049341726798726440754485546402356   -0.92608790131654640254711849102890719529    0.84138735420652599295096060100359489758   0.086203456344905850875965416526967283444
+            -0.92608790131654640254711849102890719529  -0.036956049341726798726440754485546402356    0.84138735420652599295096060100359489758   0.086203456344905850875965416526967283444
+           -0.036956049341726798726440754485546402356  -0.036956049341726798726440754485546402356    0.84138735420652599295096060100359489758   0.086203456344905850875965416526967283444
+            -0.80903350325702126768802840491559333013    0.61806700651404253537605680983118666027   -0.79584909058698306626909009833843112408    0.11671409960839383433248991439345238969
+             0.61806700651404253537605680983118666027   -0.80903350325702126768802840491559333013   -0.79584909058698306626909009833843112408    0.11671409960839383433248991439345238969
+            -0.80903350325702126768802840491559333013   -0.80903350325702126768802840491559333013   -0.79584909058698306626909009833843112408    0.11671409960839383433248991439345238969
+            -0.80903350325702126768802840491559333013    0.61806700651404253537605680983118666027    0.79584909058698306626909009833843112408    0.11671409960839383433248991439345238969
+             0.61806700651404253537605680983118666027   -0.80903350325702126768802840491559333013    0.79584909058698306626909009833843112408    0.11671409960839383433248991439345238969
+            -0.80903350325702126768802840491559333013   -0.80903350325702126768802840491559333013    0.79584909058698306626909009833843112408    0.11671409960839383433248991439345238969
+            -0.97570173680324342035479366133901540544    0.48599336414579127066992364784746385989                                           0    0.10205942534059700976786329599159184945
+             0.48599336414579127066992364784746385989   -0.97570173680324342035479366133901540544                                           0    0.10205942534059700976786329599159184945
+            -0.51029162734254785031512998650844845445    0.48599336414579127066992364784746385989                                           0    0.10205942534059700976786329599159184945
+             0.48599336414579127066992364784746385989   -0.51029162734254785031512998650844845445                                           0    0.10205942534059700976786329599159184945
+            -0.51029162734254785031512998650844845445   -0.97570173680324342035479366133901540544                                           0    0.10205942534059700976786329599159184945
+            -0.97570173680324342035479366133901540544   -0.51029162734254785031512998650844845445                                           0    0.10205942534059700976786329599159184945
+            -0.38968756713554786569859973469754143293   -0.69403080315040476305982696739689230156    -0.4039345605321099116998238341362959742    0.15576281310483042016174134706586605691
+            -0.69403080315040476305982696739689230156   -0.38968756713554786569859973469754143293    -0.4039345605321099116998238341362959742    0.15576281310483042016174134706586605691
+            0.083718370285952628758426702094433734489   -0.69403080315040476305982696739689230156    -0.4039345605321099116998238341362959742    0.15576281310483042016174134706586605691
+            -0.69403080315040476305982696739689230156   0.083718370285952628758426702094433734489    -0.4039345605321099116998238341362959742    0.15576281310483042016174134706586605691
+            0.083718370285952628758426702094433734489   -0.38968756713554786569859973469754143293    -0.4039345605321099116998238341362959742    0.15576281310483042016174134706586605691
+            -0.38968756713554786569859973469754143293   0.083718370285952628758426702094433734489    -0.4039345605321099116998238341362959742    0.15576281310483042016174134706586605691
+            -0.38968756713554786569859973469754143293   -0.69403080315040476305982696739689230156     0.4039345605321099116998238341362959742    0.15576281310483042016174134706586605691
+            -0.69403080315040476305982696739689230156   -0.38968756713554786569859973469754143293     0.4039345605321099116998238341362959742    0.15576281310483042016174134706586605691
+            0.083718370285952628758426702094433734489   -0.69403080315040476305982696739689230156     0.4039345605321099116998238341362959742    0.15576281310483042016174134706586605691
+            -0.69403080315040476305982696739689230156   0.083718370285952628758426702094433734489     0.4039345605321099116998238341362959742    0.15576281310483042016174134706586605691
+            0.083718370285952628758426702094433734489   -0.38968756713554786569859973469754143293     0.4039345605321099116998238341362959742    0.15576281310483042016174134706586605691
+            -0.38968756713554786569859973469754143293   0.083718370285952628758426702094433734489     0.4039345605321099116998238341362959742    0.15576281310483042016174134706586605691
+        ]
+    elseif n == 8
+        xw = [
+            -0.33333333333333333333333333333333333333   -0.33333333333333333333333333333333333333   -0.22443386610599057286592461466529935585   0.082939746571421834033180571412028667484
+            -0.33333333333333333333333333333333333333   -0.33333333333333333333333333333333333333    0.22443386610599057286592461466529935585   0.082939746571421834033180571412028667484
+            -0.33333333333333333333333333333333333333   -0.33333333333333333333333333333333333333   -0.68182544157086581170215168002374039522    0.20435770844872619128158195584935122282
+            -0.33333333333333333333333333333333333333   -0.33333333333333333333333333333333333333    0.68182544157086581170215168002374039522    0.20435770844872619128158195584935122282
+            -0.89317529812618572284389086719897089711    0.78635059625237144568778173439794179421                                           0   0.072608995136598867700988758862181257017
+             0.78635059625237144568778173439794179421   -0.89317529812618572284389086719897089711                                           0   0.072608995136598867700988758862181257017
+            -0.89317529812618572284389086719897089711   -0.89317529812618572284389086719897089711                                           0   0.072608995136598867700988758862181257017
+           -0.079822074372578761440178266060788079463   -0.84035585125484247711964346787842384107                                           0     0.2107989796260286741645220068699562506
+            -0.84035585125484247711964346787842384107  -0.079822074372578761440178266060788079463                                           0     0.2107989796260286741645220068699562506
+           -0.079822074372578761440178266060788079463  -0.079822074372578761440178266060788079463                                           0     0.2107989796260286741645220068699562506
+           -0.082861862418097326120000870246452008956   -0.83427627516380534775999825950709598209   -0.87619591000025424474025327234232345955   0.084346037174345969123924723988899397111
+            -0.83427627516380534775999825950709598209  -0.082861862418097326120000870246452008956   -0.87619591000025424474025327234232345955   0.084346037174345969123924723988899397111
+           -0.082861862418097326120000870246452008956  -0.082861862418097326120000870246452008956   -0.87619591000025424474025327234232345955   0.084346037174345969123924723988899397111
+           -0.082861862418097326120000870246452008956   -0.83427627516380534775999825950709598209    0.87619591000025424474025327234232345955   0.084346037174345969123924723988899397111
+            -0.83427627516380534775999825950709598209  -0.082861862418097326120000870246452008956    0.87619591000025424474025327234232345955   0.084346037174345969123924723988899397111
+           -0.082861862418097326120000870246452008956  -0.082861862418097326120000870246452008956    0.87619591000025424474025327234232345955   0.084346037174345969123924723988899397111
+            -0.65187678415125922476703025933216051413    0.30375356830251844953406051866432102826   -0.40608422935459024220472507696440370377    0.16251705063595730345650668876189645828
+             0.30375356830251844953406051866432102826   -0.65187678415125922476703025933216051413   -0.40608422935459024220472507696440370377    0.16251705063595730345650668876189645828
+            -0.65187678415125922476703025933216051413   -0.65187678415125922476703025933216051413   -0.40608422935459024220472507696440370377    0.16251705063595730345650668876189645828
+            -0.65187678415125922476703025933216051413    0.30375356830251844953406051866432102826    0.40608422935459024220472507696440370377    0.16251705063595730345650668876189645828
+             0.30375356830251844953406051866432102826   -0.65187678415125922476703025933216051413    0.40608422935459024220472507696440370377    0.16251705063595730345650668876189645828
+            -0.65187678415125922476703025933216051413   -0.65187678415125922476703025933216051413    0.40608422935459024220472507696440370377    0.16251705063595730345650668876189645828
+            -0.90552242832046124932977560924577476395     0.8110448566409224986595512184915495279   -0.82881544305868020493486275746014093078   0.028581431014154463281085418946775972209
+              0.8110448566409224986595512184915495279   -0.90552242832046124932977560924577476395   -0.82881544305868020493486275746014093078   0.028581431014154463281085418946775972209
+            -0.90552242832046124932977560924577476395   -0.90552242832046124932977560924577476395   -0.82881544305868020493486275746014093078   0.028581431014154463281085418946775972209
+            -0.90552242832046124932977560924577476395     0.8110448566409224986595512184915495279    0.82881544305868020493486275746014093078   0.028581431014154463281085418946775972209
+              0.8110448566409224986595512184915495279   -0.90552242832046124932977560924577476395    0.82881544305868020493486275746014093078   0.028581431014154463281085418946775972209
+            -0.90552242832046124932977560924577476395   -0.90552242832046124932977560924577476395    0.82881544305868020493486275746014093078   0.028581431014154463281085418946775972209
+            -0.68050147211482206696479393170526142918    0.36100294422964413392958786341052285835   -0.96586516654065435448682371732369762245   0.044572109393861233332081245114974812183
+             0.36100294422964413392958786341052285835   -0.68050147211482206696479393170526142918   -0.96586516654065435448682371732369762245   0.044572109393861233332081245114974812183
+            -0.68050147211482206696479393170526142918   -0.68050147211482206696479393170526142918   -0.96586516654065435448682371732369762245   0.044572109393861233332081245114974812183
+            -0.68050147211482206696479393170526142918    0.36100294422964413392958786341052285835    0.96586516654065435448682371732369762245   0.044572109393861233332081245114974812183
+             0.36100294422964413392958786341052285835   -0.68050147211482206696479393170526142918    0.96586516654065435448682371732369762245   0.044572109393861233332081245114974812183
+            -0.68050147211482206696479393170526142918   -0.68050147211482206696479393170526142918    0.96586516654065435448682371732369762245   0.044572109393861233332081245114974812183
+             -0.4743723986175179830847346899909858658    0.45719614360199995718682043629988661246   -0.57735026918962576450914878050195745565   0.054590116363492292384362848950462321489
+             0.45719614360199995718682043629988661246    -0.4743723986175179830847346899909858658   -0.57735026918962576450914878050195745565   0.054590116363492292384362848950462321489
+            -0.98282374498448197410208574630890074666    0.45719614360199995718682043629988661246   -0.57735026918962576450914878050195745565   0.054590116363492292384362848950462321489
+             0.45719614360199995718682043629988661246   -0.98282374498448197410208574630890074666   -0.57735026918962576450914878050195745565   0.054590116363492292384362848950462321489
+            -0.98282374498448197410208574630890074666    -0.4743723986175179830847346899909858658   -0.57735026918962576450914878050195745565   0.054590116363492292384362848950462321489
+             -0.4743723986175179830847346899909858658   -0.98282374498448197410208574630890074666   -0.57735026918962576450914878050195745565   0.054590116363492292384362848950462321489
+             -0.4743723986175179830847346899909858658    0.45719614360199995718682043629988661246    0.57735026918962576450914878050195745565   0.054590116363492292384362848950462321489
+             0.45719614360199995718682043629988661246    -0.4743723986175179830847346899909858658    0.57735026918962576450914878050195745565   0.054590116363492292384362848950462321489
+            -0.98282374498448197410208574630890074666    0.45719614360199995718682043629988661246    0.57735026918962576450914878050195745565   0.054590116363492292384362848950462321489
+             0.45719614360199995718682043629988661246   -0.98282374498448197410208574630890074666    0.57735026918962576450914878050195745565   0.054590116363492292384362848950462321489
+            -0.98282374498448197410208574630890074666    -0.4743723986175179830847346899909858658    0.57735026918962576450914878050195745565   0.054590116363492292384362848950462321489
+             -0.4743723986175179830847346899909858658   -0.98282374498448197410208574630890074666    0.57735026918962576450914878050195745565   0.054590116363492292384362848950462321489
+        ]
+    elseif n == 9
+        xw = [
+            -0.33333333333333333333333333333333333333    -0.33333333333333333333333333333333333333                                            0     0.19114798844620108508537046584425003736
+            -0.33333333333333333333333333333333333333    -0.33333333333333333333333333333333333333    -0.87583915879987761660772171906042408394      0.1060969195893273614653564550346218398
+            -0.33333333333333333333333333333333333333    -0.33333333333333333333333333333333333333     0.87583915879987761660772171906042408394      0.1060969195893273614653564550346218398
+          -0.0010699296586281549894456411507705393926    -0.99786014068274369002110871769845892121                                            0    0.049462405441895223670279761487537762472
+            -0.99786014068274369002110871769845892121  -0.0010699296586281549894456411507705393926                                            0    0.049462405441895223670279761487537762472
+          -0.0010699296586281549894456411507705393926  -0.0010699296586281549894456411507705393926                                            0    0.049462405441895223670279761487537762472
+            -0.10986499991499267966311447997302837166    -0.78027000017001464067377104005394325668    -0.52963355520607675696300739352733925666     0.15080793261275860374457825128943921743
+            -0.78027000017001464067377104005394325668    -0.10986499991499267966311447997302837166    -0.52963355520607675696300739352733925666     0.15080793261275860374457825128943921743
+            -0.10986499991499267966311447997302837166    -0.10986499991499267966311447997302837166    -0.52963355520607675696300739352733925666     0.15080793261275860374457825128943921743
+            -0.10986499991499267966311447997302837166    -0.78027000017001464067377104005394325668     0.52963355520607675696300739352733925666     0.15080793261275860374457825128943921743
+            -0.78027000017001464067377104005394325668    -0.10986499991499267966311447997302837166     0.52963355520607675696300739352733925666     0.15080793261275860374457825128943921743
+            -0.10986499991499267966311447997302837166    -0.10986499991499267966311447997302837166     0.52963355520607675696300739352733925666     0.15080793261275860374457825128943921743
+            -0.57769113257897895683607913964943464903     0.15538226515795791367215827929886929806    -0.32397688216877005165250824998597314548     0.09338246974553840864616210022748459774
+             0.15538226515795791367215827929886929806    -0.57769113257897895683607913964943464903    -0.32397688216877005165250824998597314548     0.09338246974553840864616210022748459774
+            -0.57769113257897895683607913964943464903    -0.57769113257897895683607913964943464903    -0.32397688216877005165250824998597314548     0.09338246974553840864616210022748459774
+            -0.57769113257897895683607913964943464903     0.15538226515795791367215827929886929806     0.32397688216877005165250824998597314548     0.09338246974553840864616210022748459774
+             0.15538226515795791367215827929886929806    -0.57769113257897895683607913964943464903     0.32397688216877005165250824998597314548     0.09338246974553840864616210022748459774
+            -0.57769113257897895683607913964943464903    -0.57769113257897895683607913964943464903     0.32397688216877005165250824998597314548     0.09338246974553840864616210022748459774
+            -0.91231503711701435232450671034134809014     0.82463007423402870464901342068269618029    -0.40165166857078872104499428459378049449     0.03838828638520642757038156568015552951
+             0.82463007423402870464901342068269618029    -0.91231503711701435232450671034134809014    -0.40165166857078872104499428459378049449     0.03838828638520642757038156568015552951
+            -0.91231503711701435232450671034134809014    -0.91231503711701435232450671034134809014    -0.40165166857078872104499428459378049449     0.03838828638520642757038156568015552951
+            -0.91231503711701435232450671034134809014     0.82463007423402870464901342068269618029     0.40165166857078872104499428459378049449     0.03838828638520642757038156568015552951
+             0.82463007423402870464901342068269618029    -0.91231503711701435232450671034134809014     0.40165166857078872104499428459378049449     0.03838828638520642757038156568015552951
+            -0.91231503711701435232450671034134809014    -0.91231503711701435232450671034134809014     0.40165166857078872104499428459378049449     0.03838828638520642757038156568015552951
+            -0.03450013665826621486161695118676603904    -0.93099972668346757027676609762646792192     -0.9677921254642105035410819146400981344    0.026528431202320433845464036710085263608
+            -0.93099972668346757027676609762646792192    -0.03450013665826621486161695118676603904     -0.9677921254642105035410819146400981344    0.026528431202320433845464036710085263608
+            -0.03450013665826621486161695118676603904    -0.03450013665826621486161695118676603904     -0.9677921254642105035410819146400981344    0.026528431202320433845464036710085263608
+            -0.03450013665826621486161695118676603904    -0.93099972668346757027676609762646792192      0.9677921254642105035410819146400981344    0.026528431202320433845464036710085263608
+            -0.93099972668346757027676609762646792192    -0.03450013665826621486161695118676603904      0.9677921254642105035410819146400981344    0.026528431202320433845464036710085263608
+            -0.03450013665826621486161695118676603904    -0.03450013665826621486161695118676603904      0.9677921254642105035410819146400981344    0.026528431202320433845464036710085263608
+             -0.6456435117766454460710959796587760834     0.29128702355329089214219195931755216681    -0.89078174045697066578383528302064143989    0.066719942573448407007300296577585702464
+             0.29128702355329089214219195931755216681     -0.6456435117766454460710959796587760834    -0.89078174045697066578383528302064143989    0.066719942573448407007300296577585702464
+             -0.6456435117766454460710959796587760834     -0.6456435117766454460710959796587760834    -0.89078174045697066578383528302064143989    0.066719942573448407007300296577585702464
+             -0.6456435117766454460710959796587760834     0.29128702355329089214219195931755216681     0.89078174045697066578383528302064143989    0.066719942573448407007300296577585702464
+             0.29128702355329089214219195931755216681     -0.6456435117766454460710959796587760834     0.89078174045697066578383528302064143989    0.066719942573448407007300296577585702464
+             -0.6456435117766454460710959796587760834     -0.6456435117766454460710959796587760834     0.89078174045697066578383528302064143989    0.066719942573448407007300296577585702464
+            -0.90753763917811234690103494985844043039     0.81507527835622469380206989971688086078     -0.9577062038929269606887790227035291939    0.012131173836407148944270419030858660558
+             0.81507527835622469380206989971688086078    -0.90753763917811234690103494985844043039     -0.9577062038929269606887790227035291939    0.012131173836407148944270419030858660558
+            -0.90753763917811234690103494985844043039    -0.90753763917811234690103494985844043039     -0.9577062038929269606887790227035291939    0.012131173836407148944270419030858660558
+            -0.90753763917811234690103494985844043039     0.81507527835622469380206989971688086078      0.9577062038929269606887790227035291939    0.012131173836407148944270419030858660558
+             0.81507527835622469380206989971688086078    -0.90753763917811234690103494985844043039      0.9577062038929269606887790227035291939    0.012131173836407148944270419030858660558
+            -0.90753763917811234690103494985844043039    -0.90753763917811234690103494985844043039      0.9577062038929269606887790227035291939    0.012131173836407148944270419030858660558
+            -0.88126598051326876444819587362986475675     0.44939053910599895722040631866985064464                                            0    0.087691161312839084087580790508946261895
+             0.44939053910599895722040631866985064464    -0.88126598051326876444819587362986475675                                            0    0.087691161312839084087580790508946261895
+            -0.56812455859273019277221044503998588789     0.44939053910599895722040631866985064464                                            0    0.087691161312839084087580790508946261895
+             0.44939053910599895722040631866985064464    -0.56812455859273019277221044503998588789                                            0    0.087691161312839084087580790508946261895
+            -0.56812455859273019277221044503998588789    -0.88126598051326876444819587362986475675                                            0    0.087691161312839084087580790508946261895
+            -0.88126598051326876444819587362986475675    -0.56812455859273019277221044503998588789                                            0    0.087691161312839084087580790508946261895
+             -0.5517449084932878730050536869300298641     0.49681091913207739922639584227322548803    -0.69522055387163478698981196744187399375      0.0495312141698622864915543816230467997
+             0.49681091913207739922639584227322548803     -0.5517449084932878730050536869300298641    -0.69522055387163478698981196744187399375      0.0495312141698622864915543816230467997
+            -0.94506601063878952622134215534319562393     0.49681091913207739922639584227322548803    -0.69522055387163478698981196744187399375      0.0495312141698622864915543816230467997
+             0.49681091913207739922639584227322548803    -0.94506601063878952622134215534319562393    -0.69522055387163478698981196744187399375      0.0495312141698622864915543816230467997
+            -0.94506601063878952622134215534319562393     -0.5517449084932878730050536869300298641    -0.69522055387163478698981196744187399375      0.0495312141698622864915543816230467997
+             -0.5517449084932878730050536869300298641    -0.94506601063878952622134215534319562393    -0.69522055387163478698981196744187399375      0.0495312141698622864915543816230467997
+             -0.5517449084932878730050536869300298641     0.49681091913207739922639584227322548803     0.69522055387163478698981196744187399375      0.0495312141698622864915543816230467997
+             0.49681091913207739922639584227322548803     -0.5517449084932878730050536869300298641     0.69522055387163478698981196744187399375      0.0495312141698622864915543816230467997
+            -0.94506601063878952622134215534319562393     0.49681091913207739922639584227322548803     0.69522055387163478698981196744187399375      0.0495312141698622864915543816230467997
+             0.49681091913207739922639584227322548803    -0.94506601063878952622134215534319562393     0.69522055387163478698981196744187399375      0.0495312141698622864915543816230467997
+            -0.94506601063878952622134215534319562393     -0.5517449084932878730050536869300298641     0.69522055387163478698981196744187399375      0.0495312141698622864915543816230467997
+             -0.5517449084932878730050536869300298641    -0.94506601063878952622134215534319562393     0.69522055387163478698981196744187399375      0.0495312141698622864915543816230467997
+        ]
+    elseif n == 10
+        xw = [
+            -0.33333333333333333333333333333333333333    -0.33333333333333333333333333333333333333                                           -1    0.038043108084240763009572271119163370293
+            -0.33333333333333333333333333333333333333    -0.33333333333333333333333333333333333333                                            1    0.038043108084240763009572271119163370293
+            -0.33333333333333333333333333333333333333    -0.33333333333333333333333333333333333333     -0.2838364013985715684657060577249329797     0.12804920643255787761445727849451714328
+            -0.33333333333333333333333333333333333333    -0.33333333333333333333333333333333333333      0.2838364013985715684657060577249329797     0.12804920643255787761445727849451714328
+            -0.58509925002156576804484663776694671178     0.17019850004313153608969327553389342357                                            0    0.045423158480136973856219994221831113886
+             0.17019850004313153608969327553389342357    -0.58509925002156576804484663776694671178                                            0    0.045423158480136973856219994221831113886
+            -0.58509925002156576804484663776694671178    -0.58509925002156576804484663776694671178                                            0    0.045423158480136973856219994221831113886
+           -0.087745274263836468012344092631834458624    -0.82450945147232706397531181473633108275                                            0     0.12078661810034629566292168130821741978
+            -0.82450945147232706397531181473633108275   -0.087745274263836468012344092631834458624                                            0     0.12078661810034629566292168130821741978
+           -0.087745274263836468012344092631834458624   -0.087745274263836468012344092631834458624                                            0     0.12078661810034629566292168130821741978
+            -0.83336149569989450879622615827585222698     0.66672299139978901759245231655170445397                                            0    0.070523638001611290833549221721909796384
+             0.66672299139978901759245231655170445397    -0.83336149569989450879622615827585222698                                            0    0.070523638001611290833549221721909796384
+            -0.83336149569989450879622615827585222698    -0.83336149569989450879622615827585222698                                            0    0.070523638001611290833549221721909796384
+            -0.58566729572033783230931886738991607412     0.17133459144067566461863773477983214824    -0.47829026090036813970247985116366265217     0.10435318130434033235265655019864541098
+             0.17133459144067566461863773477983214824    -0.58566729572033783230931886738991607412    -0.47829026090036813970247985116366265217     0.10435318130434033235265655019864541098
+            -0.58566729572033783230931886738991607412    -0.58566729572033783230931886738991607412    -0.47829026090036813970247985116366265217     0.10435318130434033235265655019864541098
+            -0.58566729572033783230931886738991607412     0.17133459144067566461863773477983214824     0.47829026090036813970247985116366265217     0.10435318130434033235265655019864541098
+             0.17133459144067566461863773477983214824    -0.58566729572033783230931886738991607412     0.47829026090036813970247985116366265217     0.10435318130434033235265655019864541098
+            -0.58566729572033783230931886738991607412    -0.58566729572033783230931886738991607412     0.47829026090036813970247985116366265217     0.10435318130434033235265655019864541098
+            -0.71401760691431835996629386931253079667     0.42803521382863671993258773862506159333     -0.8520666639991396363392553699935906628    0.045353049770314160210925207718221960406
+             0.42803521382863671993258773862506159333    -0.71401760691431835996629386931253079667     -0.8520666639991396363392553699935906628    0.045353049770314160210925207718221960406
+            -0.71401760691431835996629386931253079667    -0.71401760691431835996629386931253079667     -0.8520666639991396363392553699935906628    0.045353049770314160210925207718221960406
+            -0.71401760691431835996629386931253079667     0.42803521382863671993258773862506159333      0.8520666639991396363392553699935906628    0.045353049770314160210925207718221960406
+             0.42803521382863671993258773862506159333    -0.71401760691431835996629386931253079667      0.8520666639991396363392553699935906628    0.045353049770314160210925207718221960406
+            -0.71401760691431835996629386931253079667    -0.71401760691431835996629386931253079667      0.8520666639991396363392553699935906628    0.045353049770314160210925207718221960406
+          -0.0059182255205082649355625415633530872236    -0.98816354895898347012887491687329382555    -0.62749011164795345095028595959683589121    0.028126573654953943549801873601641305154
+            -0.98816354895898347012887491687329382555  -0.0059182255205082649355625415633530872236    -0.62749011164795345095028595959683589121    0.028126573654953943549801873601641305154
+          -0.0059182255205082649355625415633530872236  -0.0059182255205082649355625415633530872236    -0.62749011164795345095028595959683589121    0.028126573654953943549801873601641305154
+          -0.0059182255205082649355625415633530872236    -0.98816354895898347012887491687329382555     0.62749011164795345095028595959683589121    0.028126573654953943549801873601641305154
+            -0.98816354895898347012887491687329382555  -0.0059182255205082649355625415633530872236     0.62749011164795345095028595959683589121    0.028126573654953943549801873601641305154
+          -0.0059182255205082649355625415633530872236  -0.0059182255205082649355625415633530872236     0.62749011164795345095028595959683589121    0.028126573654953943549801873601641305154
+            -0.14443508029301579714242731927948333445    -0.71112983941396840571514536144103333111    -0.76897188255294104672722599221030791768     0.09730006751474251906683045094709204535
+            -0.71112983941396840571514536144103333111    -0.14443508029301579714242731927948333445    -0.76897188255294104672722599221030791768     0.09730006751474251906683045094709204535
+            -0.14443508029301579714242731927948333445    -0.14443508029301579714242731927948333445    -0.76897188255294104672722599221030791768     0.09730006751474251906683045094709204535
+            -0.14443508029301579714242731927948333445    -0.71112983941396840571514536144103333111     0.76897188255294104672722599221030791768     0.09730006751474251906683045094709204535
+            -0.71112983941396840571514536144103333111    -0.14443508029301579714242731927948333445     0.76897188255294104672722599221030791768     0.09730006751474251906683045094709204535
+            -0.14443508029301579714242731927948333445    -0.14443508029301579714242731927948333445     0.76897188255294104672722599221030791768     0.09730006751474251906683045094709204535
+            -0.92655062106089238633040808286062465799     0.85310124212178477266081616572124931598    -0.94536185592175700771401046143327478453   0.0075697428403356175236518069263715176348
+             0.85310124212178477266081616572124931598    -0.92655062106089238633040808286062465799    -0.94536185592175700771401046143327478453   0.0075697428403356175236518069263715176348
+            -0.92655062106089238633040808286062465799    -0.92655062106089238633040808286062465799    -0.94536185592175700771401046143327478453   0.0075697428403356175236518069263715176348
+            -0.92655062106089238633040808286062465799     0.85310124212178477266081616572124931598     0.94536185592175700771401046143327478453   0.0075697428403356175236518069263715176348
+             0.85310124212178477266081616572124931598    -0.92655062106089238633040808286062465799     0.94536185592175700771401046143327478453   0.0075697428403356175236518069263715176348
+            -0.92655062106089238633040808286062465799    -0.92655062106089238633040808286062465799     0.94536185592175700771401046143327478453   0.0075697428403356175236518069263715176348
+            -0.99621476080190120325860648054065438352     0.99582008190465265332555390082623777504                                            0   0.0031208636401827010315915455956052660089
+             0.99582008190465265332555390082623777504    -0.99621476080190120325860648054065438352                                            0   0.0031208636401827010315915455956052660089
+            -0.99960532110275145006694742028558339152     0.99582008190465265332555390082623777504                                            0   0.0031208636401827010315915455956052660089
+             0.99582008190465265332555390082623777504    -0.99960532110275145006694742028558339152                                            0   0.0031208636401827010315915455956052660089
+            -0.99960532110275145006694742028558339152    -0.99621476080190120325860648054065438352                                            0   0.0031208636401827010315915455956052660089
+            -0.99621476080190120325860648054065438352    -0.99960532110275145006694742028558339152                                            0   0.0031208636401827010315915455956052660089
+            -0.91342063244101618001407082339238449491      0.3255064864791910018042153093720481789    -0.94236346745083749482627300953387466221    0.024183540649806648348002441154357646858
+              0.3255064864791910018042153093720481789    -0.91342063244101618001407082339238449491    -0.94236346745083749482627300953387466221    0.024183540649806648348002441154357646858
+            -0.41208585403817482179014448597966368399      0.3255064864791910018042153093720481789    -0.94236346745083749482627300953387466221    0.024183540649806648348002441154357646858
+              0.3255064864791910018042153093720481789    -0.41208585403817482179014448597966368399    -0.94236346745083749482627300953387466221    0.024183540649806648348002441154357646858
+            -0.41208585403817482179014448597966368399    -0.91342063244101618001407082339238449491    -0.94236346745083749482627300953387466221    0.024183540649806648348002441154357646858
+            -0.91342063244101618001407082339238449491    -0.41208585403817482179014448597966368399    -0.94236346745083749482627300953387466221    0.024183540649806648348002441154357646858
+            -0.91342063244101618001407082339238449491      0.3255064864791910018042153093720481789     0.94236346745083749482627300953387466221    0.024183540649806648348002441154357646858
+              0.3255064864791910018042153093720481789    -0.91342063244101618001407082339238449491     0.94236346745083749482627300953387466221    0.024183540649806648348002441154357646858
+            -0.41208585403817482179014448597966368399      0.3255064864791910018042153093720481789     0.94236346745083749482627300953387466221    0.024183540649806648348002441154357646858
+              0.3255064864791910018042153093720481789    -0.41208585403817482179014448597966368399     0.94236346745083749482627300953387466221    0.024183540649806648348002441154357646858
+            -0.41208585403817482179014448597966368399    -0.91342063244101618001407082339238449491     0.94236346745083749482627300953387466221    0.024183540649806648348002441154357646858
+            -0.91342063244101618001407082339238449491    -0.41208585403817482179014448597966368399     0.94236346745083749482627300953387466221    0.024183540649806648348002441154357646858
+             -0.9234169616555928413056660388640430162     0.37513146835671456389900781390965761438    -0.33533610866739435446116498131265141599     0.05460310873669147926531122821390556603
+             0.37513146835671456389900781390965761438     -0.9234169616555928413056660388640430162    -0.33533610866739435446116498131265141599     0.05460310873669147926531122821390556603
+            -0.45171450670112172259334177504561459817     0.37513146835671456389900781390965761438    -0.33533610866739435446116498131265141599     0.05460310873669147926531122821390556603
+             0.37513146835671456389900781390965761438    -0.45171450670112172259334177504561459817    -0.33533610866739435446116498131265141599     0.05460310873669147926531122821390556603
+            -0.45171450670112172259334177504561459817     -0.9234169616555928413056660388640430162    -0.33533610866739435446116498131265141599     0.05460310873669147926531122821390556603
+             -0.9234169616555928413056660388640430162    -0.45171450670112172259334177504561459817    -0.33533610866739435446116498131265141599     0.05460310873669147926531122821390556603
+             -0.9234169616555928413056660388640430162     0.37513146835671456389900781390965761438     0.33533610866739435446116498131265141599     0.05460310873669147926531122821390556603
+             0.37513146835671456389900781390965761438     -0.9234169616555928413056660388640430162     0.33533610866739435446116498131265141599     0.05460310873669147926531122821390556603
+            -0.45171450670112172259334177504561459817     0.37513146835671456389900781390965761438     0.33533610866739435446116498131265141599     0.05460310873669147926531122821390556603
+             0.37513146835671456389900781390965761438    -0.45171450670112172259334177504561459817     0.33533610866739435446116498131265141599     0.05460310873669147926531122821390556603
+            -0.45171450670112172259334177504561459817     -0.9234169616555928413056660388640430162     0.33533610866739435446116498131265141599     0.05460310873669147926531122821390556603
+             -0.9234169616555928413056660388640430162    -0.45171450670112172259334177504561459817     0.33533610866739435446116498131265141599     0.05460310873669147926531122821390556603
+            -0.96792121274553310722108842071500064645     0.74052721727873797554079834546410281777    -0.61432350187707741015732547519719343903    0.024769538519410488660113297222678366236
+             0.74052721727873797554079834546410281777    -0.96792121274553310722108842071500064645    -0.61432350187707741015732547519719343903    0.024769538519410488660113297222678366236
+            -0.77260600453320486831970992474910217133     0.74052721727873797554079834546410281777    -0.61432350187707741015732547519719343903    0.024769538519410488660113297222678366236
+             0.74052721727873797554079834546410281777    -0.77260600453320486831970992474910217133    -0.61432350187707741015732547519719343903    0.024769538519410488660113297222678366236
+            -0.77260600453320486831970992474910217133    -0.96792121274553310722108842071500064645    -0.61432350187707741015732547519719343903    0.024769538519410488660113297222678366236
+            -0.96792121274553310722108842071500064645    -0.77260600453320486831970992474910217133    -0.61432350187707741015732547519719343903    0.024769538519410488660113297222678366236
+            -0.96792121274553310722108842071500064645     0.74052721727873797554079834546410281777     0.61432350187707741015732547519719343903    0.024769538519410488660113297222678366236
+             0.74052721727873797554079834546410281777    -0.96792121274553310722108842071500064645     0.61432350187707741015732547519719343903    0.024769538519410488660113297222678366236
+            -0.77260600453320486831970992474910217133     0.74052721727873797554079834546410281777     0.61432350187707741015732547519719343903    0.024769538519410488660113297222678366236
+             0.74052721727873797554079834546410281777    -0.77260600453320486831970992474910217133     0.61432350187707741015732547519719343903    0.024769538519410488660113297222678366236
+            -0.77260600453320486831970992474910217133    -0.96792121274553310722108842071500064645     0.61432350187707741015732547519719343903    0.024769538519410488660113297222678366236
+            -0.96792121274553310722108842071500064645    -0.77260600453320486831970992474910217133     0.61432350187707741015732547519719343903    0.024769538519410488660113297222678366236
+        ]
+    else
+        throw(ArgumentError("unsupported order for prism polyquad integration"))
+    end
+
+    # Transform from [-1,1] × [-1,1] × [-1,1] to [0,1] × [0,1] × [0,1]
+    #x
+    xw[:,1] .+= 1.0
+    xw[:,1] ./= 2.0
+    #y
+    xw[:,2] .+= 1.0
+    xw[:,2] ./= 2.0
+    #z
+    xw[:,3] .+= 1.0
+    xw[:,3] ./= 2.0
+    #w
+    xw[:,4] ./= 2^3
+
+    @show xw[:,3]
+
+    return xw
+end
diff --git a/src/Quadrature/quadrature.jl b/src/Quadrature/quadrature.jl
index a2c0de7375..69d9dbbbf3 100644
--- a/src/Quadrature/quadrature.jl
+++ b/src/Quadrature/quadrature.jl
@@ -1,5 +1,6 @@
 include("gaussquad_tri_table.jl")
 include("gaussquad_tet_table.jl")
+include("gaussquad_prism_table.jl")
 include("generate_quadrature.jl")
 
 import Base.Cartesian: @nloops, @nref, @ntuple, @nexprs
@@ -85,6 +86,7 @@ julia> getpoints(qr)
 getpoints(qr::QuadratureRule) = qr.points
 
 QuadratureRule{dim,shape}(order::Int) where {dim,shape} = QuadratureRule{dim,shape}(:legendre, order)
+QuadratureRule{3,RefPrism}(order::Int) = QuadratureRule{3,RefPrism}(:polyquad, order) # No legendre rule available yet for prism...
 
 # Special case for face integration of 1D problems
 function (::Type{QuadratureRule{0, RefCube}})(quad_type::Symbol, order::Int)
@@ -162,3 +164,20 @@ function (::Type{QuadratureRule{1,RefTetrahedron}})(quad_type::Symbol, order::In
     end
     return QuadratureRule{1,RefTetrahedron,Float64}(weights, points)
 end
+
+# Grab prism quadrature rule from table
+function (::Type{QuadratureRule{3, RefPrism}})(quad_type::Symbol, order::Int)
+    if quad_type == :polyquad
+        data = _get_gauss_prismdata_polyquad(order)
+    else
+        throw(ArgumentError("unsupported quadrature rule"))
+    end
+    n_points = size(data,1)
+    points = Vector{Vec{3,Float64}}(undef, n_points)
+
+    for p in 1:size(data, 1)
+        points[p] = Vec{3,Float64}(@ntuple 3 i -> data[p, i])
+    end
+    weights = data[:, 4]
+    QuadratureRule{3,RefPrism,Float64}(weights, points)
+end
diff --git a/src/deprecations.jl b/src/deprecations.jl
index da4358d54b..11b2a1436d 100644
--- a/src/deprecations.jl
+++ b/src/deprecations.jl
@@ -1,4 +1,9 @@
 Base.@deprecate_binding DirichletBoundaryConditions ConstraintHandler
 Base.@deprecate_binding DirichletBoundaryCondition Dirichlet
+
 import Base: push!
 @deprecate push!(dh::AbstractDofHandler, args...) add!(dh, args...)
+
+@deprecate vertices(ip::Interpolation) vertexdof_indices(ip) false
+@deprecate faces(ip::Interpolation) facedof_indices(ip) false
+@deprecate edges(ip::Interpolation) edgedof_indices(ip) false
diff --git a/src/exports.jl b/src/exports.jl
index 83e3d9d6ad..b735672237 100644
--- a/src/exports.jl
+++ b/src/exports.jl
@@ -3,6 +3,7 @@ export
     Interpolation,
     RefCube,
     RefTetrahedron,
+    RefPrism,
     BubbleEnrichedLagrange,
     CrouzeixRaviart,
     Lagrange,
@@ -55,6 +56,7 @@ export
     QuadraticTetrahedron,
     Hexahedron,
     #QuadraticHexahedron,
+    Wedge,
     CellIndex,
     FaceIndex,
     EdgeIndex,
diff --git a/src/interpolations.jl b/src/interpolations.jl
index 8a10b7c301..6082a37fbf 100644
--- a/src/interpolations.jl
+++ b/src/interpolations.jl
@@ -24,6 +24,8 @@ The following interpolations are implemented:
 * `Lagrange{3,RefCube,2}`
 * `Lagrange{3,RefTetrahedron,1}`
 * `Lagrange{3,RefTetrahedron,2}`
+* `Lagrange{3,RefPrism,1}`
+* `Lagrange{3,RefPrism,2}`
 * `Serendipity{2,RefCube,2}`
 * `Serendipity{3,RefCube,2}`
 
@@ -38,6 +40,72 @@ julia> getnbasefunctions(ip)
 """
 abstract type Interpolation{dim,shape,order} end
 
+# TODO: Remove: this is a hotfix to apply constraints to embedded elements.
+edges(ip::Interpolation{2}) = faces(ip)
+edgedof_indices(ip::Interpolation{2}) = facedof_indices(ip)
+edgedof_interior_indices(ip::Interpolation{2}) = facedof_interior_indices(ip)
+facedof_indices(ip::Interpolation{1}) = vertexdof_indices(ip)
+
+"""
+    InterpolationInfo
+
+Gathers all the information needed to distribute dofs for a given interpolation. Note that
+this cache is of the same type no matter the interpolation: the purpose is to make
+dof-distribution type-stable.
+"""
+struct InterpolationInfo
+    nvertexdofs::Vector{Int}
+    nedgedofs::Vector{Int}
+    nfacedofs::Vector{Int}
+    ncelldofs::Int
+    dim::Int
+    function InterpolationInfo(interpolation::Interpolation{3})
+        new(
+            [length(i) for i ∈ vertexdof_indices(interpolation)],
+            [length(i) for i ∈ edgedof_interior_indices(interpolation)],
+            [length(i) for i ∈ facedof_interior_indices(interpolation)],
+            length(celldof_interior_indices(interpolation)),
+            3,
+        )
+    end
+    function InterpolationInfo(interpolation::Interpolation{2})
+        new(
+            [length(i) for i ∈ vertexdof_indices(interpolation)],
+            Int[],
+            [length(i) for i ∈ facedof_interior_indices(interpolation)],
+            length(celldof_interior_indices(interpolation)),
+            2,
+        )
+    end
+    function InterpolationInfo(interpolation::Interpolation{1})
+        new(
+            [length(i) for i ∈ vertexdof_indices(interpolation)],
+            Int[],
+            [0, 0], # Well. Apparently vertices are also faces. :)
+            length(celldof_interior_indices(interpolation)),
+            1,
+        )
+    end
+end
+
+# Some redundant information about the geometry of the reference cells.
+nfaces(::Interpolation{dim, RefCube}) where {dim}= 2*dim
+nfaces(::Interpolation{2, RefTetrahedron})       = 3
+nfaces(::Interpolation{3, RefTetrahedron})       = 4
+nfaces(::Interpolation{3, RefPrism})             = 5
+
+nedges(::Interpolation{1, RefCube})          =  0
+nedges(::Interpolation{2, RefCube})          =  0
+nedges(::Interpolation{3, RefCube})          = 12
+nedges(::Interpolation{2, RefTetrahedron})   =  0
+nedges(::Interpolation{3, RefTetrahedron})   =  6
+nedges(::Interpolation{3, RefPrism})         =  9
+
+nvertices(::Interpolation{dim, RefCube}) where {dim} = 2^dim
+nvertices(::Interpolation{2, RefTetrahedron})        = 3
+nvertices(::Interpolation{3, RefTetrahedron})        = 4
+nvertices(::Interpolation{3, RefPrism})              = 6
+
 Base.copy(ip::Interpolation) = ip
 
 """
@@ -94,52 +162,12 @@ Return the number of base functions for the interpolation `ip`.
 """
 getnbasefunctions(::Interpolation)
 
-# struct that gathers all the information needed to distribute
-# dofs for a given interpolation.
-struct InterpolationInfo
-    # TODO: Can be smaller than `Int` if that matters...
-    nvertexdofs::Int
-    nedgedofs::Int
-    nfacedofs::Int
-    ncelldofs::Int
-    dim::Int
-    InterpolationInfo(interpolation::Interpolation{dim}) where {dim} =
-        new(nvertexdofs(interpolation), nedgedofs(interpolation),
-            nfacedofs(interpolation),   ncelldofs(interpolation), dim)
-end
-
 # The following functions are used to distribute the dofs. Definitions:
 #   vertexdof: dof on a "corner" of the reference shape
 #   facedof: dof in the dim-1 dimension (line in 2D, surface in 3D)
 #   edgedof: dof on a line between 2 vertices (i.e. "corners") (3D only)
 #   celldof: dof that is local to the element
 
-# Necessary for correct distribution of dofs.
-"""
-    Ferrite.nvertexdofs(::Interpolation)
-
-Number of dofs per vertex.
-"""
-nvertexdofs(::Interpolation) = 0
-"""
-    Ferrite.nedgedofs(::Interpolation)
-
-Number of dofs on the interior of a single edge.
-"""
-nedgedofs(::Interpolation)   = 0
-"""
-    Ferrite.nfacedofs(::Interpolation)
-
-Number of dofs on the interior of a single face.
-"""
-nfacedofs(::Interpolation)   = 0
-"""
-    Ferrite.ncelldofs(::Interpolation)
-
-Total number of dofs in the interior of the element.
-"""
-ncelldofs(::Interpolation)   = 0
-
 """
     value(ip::Interpolation, i::Int, ξ::Vec)
 
@@ -154,9 +182,9 @@ indices of [`reference_coordinates(::Interpolation)`](@ref).
 value(ip::Interpolation, i::Int, ξ::Vec)
 
 """
-    reference_coordinates(::Interpolation)
+    reference_coordinates(ip::Interpolation)
 
-Returns a vector of coordinates with length [`getnbasefunctions(::Interpolation)`](@ref) 
+Returns a vector of coordinates with length [`getnbasefunctions(::Interpolation)`](@ref)
 and indices corresponding to the indices of a dof in [`vertices`](@ref), [`faces`](@ref) and
 [`edges`](@ref).
 
@@ -164,41 +192,78 @@ and indices corresponding to the indices of a dof in [`vertices`](@ref), [`faces
     
     TODO: Separate nodal and non-nodal interpolations.
 """
-reference_coordinates(ip::Interpolation)
+reference_coordinates(::Interpolation)
 
 """
-    vertices(::Interpolation)
+    vertexdof_indices(ip::Interpolation)
 
-A tuple containing tuples of local dof indices for the respective 
-vertex in local enumeration on a cell defined by [`vertices(::Cell)`](@ref). The vertex enumeration must 
+A tuple containing tuples of local dof indices for the respective vertex in local
+enumeration on a cell defined by [`vertices(::Cell)`](@ref). The vertex enumeration must
 match the vertex enumeration of the corresponding geometrical cell.
-The length of each tuple must be exactly [`nvertexdofs(::Interpolation)`](@ref).
 """
-vertices(::Interpolation)
+vertexdof_indices(ip::Interpolation) = ntuple(_ -> (), nvertices(ip))
+
+"""
+    edgedof_indices(ip::Interpolation)
+
+A tuple containing tuples of local dof indices for the respective edge in local enumeration
+on a cell defined by [`edges(::Cell)`](@ref). The edge enumeration must match the edge
+enumeration of the corresponding geometrical cell.
 """
-    edges(::Interpolation)
+edgedof_indices(::Interpolation)
 
-A tuple containing tuples of local dof indices for the respective 
-edge in local enumeration on a cell defined by [`edges(::Cell)`](@ref). The edge enumeration must 
-match the edge enumeration of the corresponding geometrical cell.
-Note that the vertex dofs are included here. The length of 
-each tuple must be exactly [`nvertexdofs(::Interpolation)`](@ref)+[`nedgedofs(::Interpolation)`](@ref).
 """
-edges(::Interpolation)
+    edgedof_interior_indices(ip::Interpolation)
+
+A tuple containing tuples of the local dof indices on the interior of the respective edge in
+local enumeration on a cell defined by [`edges(::Cell)`](@ref). The edge enumeration must
+match the edge enumeration of the corresponding geometrical cell. Note that the vertex dofs
+are included here.
 """
-    faces(::Interpolation)
+edgedof_interior_indices(::Interpolation)
 
-A tuple containing tuples of all local dof indices for the respective 
-face in local enumeration on a cell defined by [`faces(::Cell)`](@ref). The face enumeration must 
-match the face enumeration of the corresponding geometrical cell.
-Note that the vertex and edge dofs are included here. The
-length of each tuple must be exactly [`nvertexdofs(::Interpolation)`](@ref)+[`nedgedofs(::Interpolation)`](@ref)+[`nfacedofs(::Interpolation)`](@ref).
 """
-faces(::Interpolation)
+    facedof_indices(ip::Interpolation)
 
-# Needed for distributing dofs on shells correctly (face in 2d is edge in 3d)
-edges(ip::Interpolation{2}) = faces(ip)
-nedgedofs(ip::Interpolation{2}) = nfacedofs(ip)
+A tuple containing tuples of all local dof indices for the respective face in local
+enumeration on a cell defined by [`faces(::Cell)`](@ref). The face enumeration must match
+the face enumeration of the corresponding geometrical cell.
+"""
+facedof_indices(::Interpolation)
+
+"""
+    facedof_interior_indices(ip::Interpolation)
+
+A tuple containing tuples of the local dof indices on the interior of the respective face in
+local enumeration on a cell defined by [`faces(::Cell)`](@ref). The face enumeration must
+match the face enumeration of the corresponding geometrical cell. Note that the vertex and
+edge dofs are included here.
+"""
+facedof_interior_indices(::Interpolation) 
+
+"""
+    celldof_interior_indices(ip::Interpolation)
+
+Tuple containing the dof indices associated with the interior of the cell.
+"""
+celldof_interior_indices(::Interpolation) = ()
+
+# Some helpers to skip boilerplate
+edgedof_indices(ip::Interpolation{3}) = ntuple(_ -> (), nedges(ip))
+edgedof_interior_indices(ip::Interpolation{3}) = ntuple(_ -> (), nedges(ip))
+facedof_indices(ip::Union{Interpolation{2}, Interpolation{3}}) =  ntuple(_ -> (), nfaces(ip))
+facedof_interior_indices(ip::Union{Interpolation{2}, Interpolation{3}}) =  ntuple(_ -> (), nfaces(ip))
+
+"""
+    boundarydof_indices(::Type{<:BoundaryIndex})
+
+Helper function to generically dispatch on the correct dof sets of a boundary entity.
+"""
+boundarydof_indices(::Type{<:BoundaryIndex})
+
+boundarydof_indices(::Type{FaceIndex}) = Ferrite.facedof_indices
+boundarydof_indices(::Type{EdgeIndex}) = Ferrite.edgedof_indices
+boundarydof_indices(::Type{VertexIndex}) = Ferrite.vertexdof_indices
 
 #########################
 # DiscontinuousLagrange #
@@ -212,17 +277,11 @@ struct DiscontinuousLagrange{dim,shape,order} <: Interpolation{dim,shape,order}
 getlowerdim(::DiscontinuousLagrange{dim,shape,order}) where {dim,shape,order} = DiscontinuousLagrange{dim-1,shape,order}()
 getlowerorder(::DiscontinuousLagrange{dim,shape,order}) where {dim,shape,order} = DiscontinuousLagrange{dim,shape,order-1}()
 
-# TODO generalize properly
-# ncelldofs(::DiscontinuousLagrange{dim,RefCube,order}) where {dim,order} = (order+1)^dim
-# ncelldofs(::DiscontinuousLagrange{2,RefTetrahedron,order}) where {order} = (order+1)*(order+2)/2
-# ncelldofs(::DiscontinuousLagrange{3,RefTetrahedron,order}) where {order} = (order+1)*(order+2)/2
 getnbasefunctions(::DiscontinuousLagrange{dim,shape,order}) where {dim,shape,order} = getnbasefunctions(Lagrange{dim,shape,order}())
-ncelldofs(::DiscontinuousLagrange{dim,shape,order}) where {dim,shape,order} = getnbasefunctions(DiscontinuousLagrange{dim,shape,order}())
-
 getnbasefunctions(::DiscontinuousLagrange{dim,shape,0}) where {dim,shape} = 1
-ncelldofs(::DiscontinuousLagrange{dim,shape,0}) where {dim,shape} = 1
 
-faces(::DiscontinuousLagrange{dim,shape,order}) where {dim,shape,order} = ()
+# This just moves all dofs into the interior of the element.
+celldof_interior_indices(ip::DiscontinuousLagrange{dim,shape,order}) where {dim,shape,order} = ntuple(i->i, getnbasefunctions(ip))
 
 # Mirror the Lagrange element for now.
 function reference_coordinates(ip::DiscontinuousLagrange{dim,shape,order}) where {dim,shape,order}
@@ -246,6 +305,7 @@ function reference_coordinates(ip::DiscontinuousLagrange{3,RefTetrahedron,0})
 end
 
 function value(ip::DiscontinuousLagrange{dim,shape,0}, i::Int, ξ::Vec{dim}) where {dim,shape}
+    i > 1 && throw(ArgumentError("no shape function $i for interpolation $ip"))
     return 1.0
 end
 
@@ -255,11 +315,12 @@ end
 struct Lagrange{dim,shape,order} <: Interpolation{dim,shape,order} end
 
 # Vertices for all Lagrange interpolations are the same
-vertices(::Lagrange{2,RefCube,order}) where {order} = ((1,),(2,),(3,),(4,))
-vertices(::Lagrange{3,RefCube,order}) where {order} = ((1,),(2,),(3,),(4,),(5,),(6,),(7,),(8,))
-vertices(::Lagrange{2,RefTetrahedron,order}) where {order} = ((1,),(2,),(3,))
-vertices(::Lagrange{3,RefTetrahedron,order}) where {order} = ((1,),(2,),(3,),(4,))
-nvertexdofs(::Lagrange{dim,shape,order}) where {dim,shape,order} = 1
+vertexdof_indices(::Lagrange{1,RefCube,order}) where {order} = ((1,),(2,))
+vertexdof_indices(::Lagrange{2,RefCube,order}) where {order} = ((1,),(2,),(3,),(4,))
+vertexdof_indices(::Lagrange{3,RefCube,order}) where {order} = ((1,),(2,),(3,),(4,),(5,),(6,),(7,),(8,))
+vertexdof_indices(::Lagrange{2,RefTetrahedron,order}) where {order} = ((1,),(2,),(3,))
+vertexdof_indices(::Lagrange{3,RefTetrahedron,order}) where {order} = ((1,),(2,),(3,),(4,))
+vertexdof_indices(::Lagrange{3,RefPrism,order}) where {order} = ((1,), (2,), (3,), (4,), (5,), (6,))
 
 getlowerdim(::Lagrange{dim,shape,order}) where {dim,shape,order} = Lagrange{dim-1,shape,order}()
 getlowerorder(::Lagrange{dim,shape,order}) where {dim,shape,order} = Lagrange{dim,shape,order-1}()
@@ -270,8 +331,6 @@ getlowerorder(::Lagrange{dim,shape,1}) where {dim,shape} = DiscontinuousLagrange
 ##################################
 getnbasefunctions(::Lagrange{1,RefCube,1}) = 2
 
-faces(::Lagrange{1,RefCube,1}) = ((1,), (2,))
-
 function reference_coordinates(::Lagrange{1,RefCube,1})
     return [Vec{1, Float64}((-1.0,)),
             Vec{1, Float64}(( 1.0,))]
@@ -288,9 +347,9 @@ end
 # Lagrange dim 1 RefCube order 2 #
 ##################################
 getnbasefunctions(::Lagrange{1,RefCube,2}) = 3
-ncelldofs(::Lagrange{1,RefCube,2}) = 1
 
-faces(::Lagrange{1,RefCube,2}) = ((1,), (2,))
+facedof_indices(::Lagrange{1,RefCube,2}) = ((1,), (2,))
+celldof_interior_indices(::Lagrange{1,RefCube,2}) = (3,)
 
 function reference_coordinates(::Lagrange{1,RefCube,2})
     return [Vec{1, Float64}((-1.0,)),
@@ -311,7 +370,7 @@ end
 ##################################
 getnbasefunctions(::Lagrange{2,RefCube,1}) = 4
 
-faces(::Lagrange{2,RefCube,1}) = ((1,2), (2,3), (3,4), (4,1))
+facedof_indices(::Lagrange{2,RefCube,1}) = ((1,2), (2,3), (3,4), (4,1))
 
 function reference_coordinates(::Lagrange{2,RefCube,1})
     return [Vec{2, Float64}((-1.0, -1.0)),
@@ -334,10 +393,10 @@ end
 # Lagrange dim 2 RefCube order 2 #
 ##################################
 getnbasefunctions(::Lagrange{2,RefCube,2}) = 9
-nfacedofs(::Lagrange{2,RefCube,2}) = 1
-ncelldofs(::Lagrange{2,RefCube,2}) = 1
 
-faces(::Lagrange{2,RefCube,2}) = ((1,2,5), (2,3,6), (3,4,7), (4,1,8))
+facedof_indices(::Lagrange{2,RefCube,2}) = ((1,2, 5), (2,3, 6), (3,4, 7), (4,1, 8))
+facedof_interior_indices(::Lagrange{2,RefCube,2}) = ((5,), (6,), (7,), (8,))
+celldof_interior_indices(::Lagrange{2,RefCube,2}) = (9,)
 
 function reference_coordinates(::Lagrange{2,RefCube,2})
     return [Vec{2, Float64}((-1.0, -1.0)),
@@ -372,8 +431,7 @@ end
 getnbasefunctions(::Lagrange{2,RefTetrahedron,1}) = 3
 getlowerdim(::Lagrange{2, RefTetrahedron, order}) where {order} = Lagrange{1, RefCube, order}()
 
-vertices(::Lagrange{2,RefTetrahedron,1}) = (1,2,3)
-faces(::Lagrange{2,RefTetrahedron,1}) = ((1,2), (2,3), (3,1))
+facedof_indices(::Lagrange{2,RefTetrahedron,1}) = ((1,2), (2,3), (3,1))
 
 function reference_coordinates(::Lagrange{2,RefTetrahedron,1})
     return [Vec{2, Float64}((1.0, 0.0)),
@@ -394,10 +452,9 @@ end
 # Lagrange dim 2 RefTetrahedron order 2 #
 #########################################
 getnbasefunctions(::Lagrange{2,RefTetrahedron,2}) = 6
-nfacedofs(::Lagrange{2,RefTetrahedron,2}) = 1
 
-vertices(::Lagrange{2,RefTetrahedron,2}) = (1,2,3)
-faces(::Lagrange{2,RefTetrahedron,2}) = ((1,2,4), (2,3,5), (3,1,6))
+facedof_indices(::Lagrange{2,RefTetrahedron,2}) = ((1,2,4), (2,3,5), (3,1,6))
+facedof_interior_indices(::Lagrange{2,RefTetrahedron,2}) = ((4,), (5,), (6,))
 
 function reference_coordinates(::Lagrange{2,RefTetrahedron,2})
     return [Vec{2, Float64}((1.0, 0.0)),
@@ -436,16 +493,8 @@ function getnbasefunctions(ip::Lagrange2Tri345)
     return (order + 1) * (order + 2) ÷ 2
 end
 
-nvertexdofs(::Lagrange2Tri345) = 1
-nedgedofs(ip::Lagrange2Tri345) = getorder(ip) - 1
-nfacedofs(ip::Lagrange2Tri345) = getorder(ip) - 1
-function ncelldofs(ip::Lagrange2Tri345)
-    order = getorder(ip)
-    return (order + 1) * (order + 2) ÷ 2 - 3 * order
-end
-
 # Permutation to switch numbering to Ferrite ordering
-const permdof2D = Dict{Int,Vector{Int}}(
+const permdof2DLagrange2Tri345 = Dict{Int,Vector{Int}}(
     1 => [1, 2, 3],
     2 => [3, 6, 1, 5, 4, 2],
     3 => [4, 10, 1, 7, 9, 8, 5, 2, 3, 6],
@@ -453,23 +502,35 @@ const permdof2D = Dict{Int,Vector{Int}}(
     5 => [6, 21, 1, 11, 15, 18, 20, 19, 16, 12, 7, 2, 3, 4, 5, 8, 9, 10, 13, 14, 17],
 )
 
-vertices(::Lagrange2Tri345) = (1, 2, 3)
+vertexdof_indices(::Lagrange2Tri345) = (1, 2, 3)
 
-function faces(ip::Lagrange2Tri345)
+function facedof_indices(ip::Lagrange2Tri345)
     order = getorder(ip)
-    if order == 1
-        return ((1,2), (2,3), (3,1))
-    elseif order == 2
-        return ((1,2,4), (2,3,5), (3,1,6))
-    elseif order == 3
-        return ((1,2,4,5), (2,3,6,7), (3,1,8,9))
-    elseif order == 4
-        return ((1,2,4,5,6), (2,3,7,8,9), (3,1,10,11,12))
-    elseif order == 5
-        return ((1,2,4,5,6,7), (2,3,8,9,10,11), (3,1,12,13,14,15))
-    else
-        error("unreachable")
-    end
+    order == 1 && return ((1,2), (2,3), (3,1))
+    order == 2 && return ((1,2,4), (2,3,5), (3,1,6))
+    order == 3 && return ((1,2,4,5), (2,3,6,7), (3,1,8,9))
+    order == 4 && return ((1,2,4,5,6), (2,3,7,8,9), (3,1,10,11,12))
+    order == 5 && return ((1,2,4,5,6,7), (2,3,8,9,10,11), (3,1,12,13,14,15))
+
+    throw(ArgumentError("Unsupported order $order for Lagrange on triangles."))
+end
+
+function facedof_interior_indices(ip::Lagrange2Tri345)
+    order = getorder(ip)
+    order == 1 && return ((), (), ())
+    order == 2 && return ((4,), (5,), (6,))
+    order == 3 && return ((4,5), (6,7), (8,9))
+    order == 4 && return ((4,5,6), (7,8,9), (10,11,12))
+    order == 5 && return ((4,5,6,7), (8,9,10,11), (12,13,14,15))
+    
+    throw(ArgumentError("Unsupported order $order for Lagrange on triangles."))
+end
+
+function celldof_interior_indices(ip::Lagrange2Tri345)
+    order = getorder(ip)
+    ncellintdofs = (order + 1) * (order + 2) ÷ 2 - 3 * order
+    totaldofs = getnbasefunctions(ip)
+    return ntuple(i->totaldofs-ncellintdofs+i, ncellintdofs)
 end
 
 function reference_coordinates(ip::Lagrange2Tri345)
@@ -480,12 +541,15 @@ function reference_coordinates(ip::Lagrange2Tri345)
             push!(coordpts, Vec{2, Float64}((l / order, k / order)))
         end
     end
-    return permute!(coordpts, permdof2D[order])
+    return permute!(coordpts, permdof2DLagrange2Tri345[order])
 end
 
 function value(ip::Lagrange2Tri345, i::Int, ξ::Vec{2})
+    if !(0 < i <= getnbasefunctions(ip))
+        throw(ArgumentError("no shape function $i for interpolation $ip"))
+    end
     order = getorder(ip)
-    i = permdof2D[order][i]
+    i = permdof2DLagrange2Tri345[order][i]
     ξ_x = ξ[1]
     ξ_y = ξ[2]
     γ = 1. - ξ_x - ξ_y
@@ -517,8 +581,8 @@ end
 #########################################
 getnbasefunctions(::Lagrange{3,RefTetrahedron,1}) = 4
 
-faces(::Lagrange{3,RefTetrahedron,1}) = ((1,3,2), (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))
+facedof_indices(::Lagrange{3,RefTetrahedron,1}) = ((1,3,2), (1,2,4), (2,3,4), (1,4,3))
+edgedof_indices(::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)),
@@ -542,10 +606,10 @@ end
 # Lagrange dim 3 RefTetrahedron order 2 #
 #########################################
 getnbasefunctions(::Lagrange{3,RefTetrahedron,2}) = 10
-nedgedofs(::Lagrange{3,RefTetrahedron,2}) = 1
 
-faces(::Lagrange{3,RefTetrahedron,2}) = ((1,3,2,7,6,5), (1,2,4,5,9,8), (2,3,4,6,10,9), (1,4,3,8,10,7))
-edges(::Lagrange{3,RefTetrahedron,2}) = ((1,2,5), (2,3,6), (3,1,7), (1,4,8), (2,4,9), (3,4,10))
+facedof_indices(::Lagrange{3,RefTetrahedron,2}) = ((1,3,2,7,6,5), (1,2,4,5,9,8), (2,3,4,6,10,9), (1,4,3,8,10,7))
+edgedof_indices(::Lagrange{3,RefTetrahedron,2}) = ((1,2,5), (2,3,6), (3,1,7), (1,4,8), (2,4,9), (3,4,10))
+edgedof_interior_indices(::Lagrange{3,RefTetrahedron,2}) = ((5,), (6,), (7,), (8,), (9,), (10,))
 
 function reference_coordinates(::Lagrange{3,RefTetrahedron,2})
     return [Vec{3, Float64}((0.0, 0.0, 0.0)),
@@ -584,8 +648,8 @@ end
 ##################################
 getnbasefunctions(::Lagrange{3,RefCube,1}) = 8
 
-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))
+facedof_indices(::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))
+edgedof_indices(::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)),
@@ -619,11 +683,8 @@ end
 ##################################
 # Based on vtkTriQuadraticHexahedron (see https://kitware.github.io/vtk-examples/site/Cxx/GeometricObjects/IsoparametricCellsDemo/)
 getnbasefunctions(::Lagrange{3,RefCube,2}) = 27
-nedgedofs(::Lagrange{3,RefCube,2}) = 1
-nfacedofs(::Lagrange{3,RefCube,2}) = 1
-ncelldofs(::Lagrange{3,RefCube,2}) = 1
 
-faces(::Lagrange{3,RefCube,2}) = (
+facedof_indices(::Lagrange{3,RefCube,2}) = (
     (1,4,3,2, 12,11,10,9, 21),
     (1,2,6,5, 9,18,13,17, 22),
     (2,3,7,6, 10,19,14,18, 23),
@@ -631,7 +692,29 @@ faces(::Lagrange{3,RefCube,2}) = (
     (1,5,8,4, 17,16,20,12, 25),
     (5,6,7,8, 13,14,15,16, 26),
 )
-edges(::Lagrange{3,RefCube,2}) = ((1,2, 9), (2,3, 10), (3,4, 11), (4,1, 12), (1,5, 17), (2,6, 18), (3,7, 19), (4,8, 20), (5,6, 13), (6,7, 14), (7,8, 15), (8,5, 16))
+facedof_interior_indices(::Lagrange{3,RefCube,2}) = (
+    (21,), (22,), (23,), (24,), (25,), (26,),
+)
+
+edgedof_indices(::Lagrange{3,RefCube,2}) = (
+    (1,2, 9),
+    (2,3, 10),
+    (3,4, 11),
+    (4,1, 12),
+    (1,5, 17),
+    (2,6, 18),
+    (3,7, 19),
+    (4,8, 20),
+    (5,6, 13),
+    (6,7, 14),
+    (7,8, 15),
+    (8,5, 16)
+)
+edgedof_interior_indices(::Lagrange{3,RefCube,2}) = (
+    (9,), (10,), (11,), (12,), (17), (18,), (19,), (20,), (13,), (14,), (15,), (16,)
+)
+
+celldof_interior_indices(::Lagrange{3,RefCube,2}) = (27,)
 
 function reference_coordinates(::Lagrange{3,RefCube,2})
            # vertex
@@ -707,6 +790,131 @@ function value(ip::Lagrange{3,RefCube,2}, i::Int, ξ::Vec{3, T}) where {T}
     throw(ArgumentError("no shape function $i for interpolation $ip"))
 end
 
+
+###################################
+# Lagrange dim 3 RefPrism order 1 #
+###################################
+# Build on https://defelement.com/elements/examples/prism-Lagrange-1.html
+getnbasefunctions(::Lagrange{3,RefPrism,1}) = 6
+
+facedof_indices(::Lagrange{3,RefPrism,1}) = ((1,3,2), (1,2,5,4), (3,1,4,6), (2,3,6,5), (4,5,6))
+edgedof_indices(::Lagrange{3,RefPrism,1}) = ((2,1), (1,3), (1,4), (3,2), (2,5), (3,6), (4,5), (4,6), (6,5))
+
+function reference_coordinates(::Lagrange{3,RefPrism,1})
+    return [Vec{3, Float64}((0.0, 0.0, 0.0)),
+            Vec{3, Float64}((1.0, 0.0, 0.0)),
+            Vec{3, Float64}((0.0, 1.0, 0.0)),
+            Vec{3, Float64}((0.0, 0.0, 1.0)),
+            Vec{3, Float64}((1.0, 0.0, 1.0)),
+            Vec{3, Float64}((0.0, 1.0, 1.0))]
+end
+
+function value(ip::Lagrange{3,RefPrism,1}, i::Int, ξ::Vec{3})
+    (x,y,z) = ξ
+    i == 1 && return 1-x-y -z*(1-x-y)
+    i == 2 && return x*(1-z)
+    i == 3 && return y*(1-z)
+    i == 4 && return z*(1-x-y)
+    i == 5 && return x*z
+    i == 6 && return y*z
+    throw(ArgumentError("no shape function $i for interpolation $ip"))
+end
+
+###################################
+# Lagrange dim 3 RefPrism order 2 #
+###################################
+# Build on https://defelement.com/elements/examples/prism-Lagrange-2.html .
+# This is simply the tensor-product of a quadratic triangle with a quadratic line.
+getnbasefunctions(::Lagrange{3,RefPrism,2}) = 18
+
+facedof_indices(::Lagrange{3,RefPrism,2}) = (
+    #Vertices| Edges  | Face 
+    (1,3,2  , 8,10,7         ),
+    (1,2,5,4, 7,11,13,9,   16), 
+    (3,1,4,6, 8,9,14,12,   17),
+    (2,3,6,5, 10,12,15,11, 18),
+    (4,5,6  , 13,15,14       ),
+)
+facedof_interior_indices(::Lagrange{3,RefPrism,2}) = (
+    #Vertices| Edges  | Face 
+    (), 
+    (16,), 
+    (17,), 
+    (18,), 
+    (),
+)
+edgedof_indices(::Lagrange{3,RefPrism,2}) = (
+    #Vert|Edge
+    (2,1, 7),
+    (1,3, 8),
+    (1,4, 9),
+    (3,2, 10),
+    (2,5, 11),
+    (3,6, 12),
+    (4,5, 13),
+    (4,6, 14),
+    (6,5, 15),
+)
+edgedof_interior_indices(::Lagrange{3,RefPrism,2}) = (
+    #Vert|Edge
+    (7,),
+    (8,),
+    (9,),
+    (10,),
+    (11,),
+    (12,),
+    (13,),
+    (14,),
+    (15,),
+)
+
+function reference_coordinates(::Lagrange{3,RefPrism,2})
+    return [Vec{3, Float64}((0.0, 0.0, 0.0)),
+            Vec{3, Float64}((1.0, 0.0, 0.0)),
+            Vec{3, Float64}((0.0, 1.0, 0.0)),
+            Vec{3, Float64}((0.0, 0.0, 1.0)),
+            Vec{3, Float64}((1.0, 0.0, 1.0)),
+            Vec{3, Float64}((0.0, 1.0, 1.0)),
+            Vec{3, Float64}((1/2, 0.0, 0.0)),
+            Vec{3, Float64}((0.0, 1/2, 0.0)),
+            Vec{3, Float64}((0.0, 0.0, 1/2)),
+            Vec{3, Float64}((1/2, 1/2, 0.0)),
+            Vec{3, Float64}((1.0, 0.0, 1/2)),
+            Vec{3, Float64}((0.0, 1.0, 1/2)),
+            Vec{3, Float64}((1/2, 0.0, 1.0)),
+            Vec{3, Float64}((0.0, 1/2, 1.0)),
+            Vec{3, Float64}((1/2, 1/2, 1.0)),
+            Vec{3, Float64}((1/2, 0.0, 1/2)),
+            Vec{3, Float64}((0.0, 1/2, 1/2)),
+            Vec{3, Float64}((1/2, 1/2, 1/2)),]
+end
+
+function value(ip::Lagrange{3,RefPrism,2}, i::Int, ξ::Vec{3})
+    (x,y,z) = ξ
+    x² = x*x
+    y² = y*y
+    z² = z*z
+    i == 1  && return 4*x²*z² - 6x²*z +2x² +8x*y*z² -12x*y*z +4x*y -6x*z² +9x*z -3x +4y²*z² -6y²*z + 2y² -6y*z² +9y*z -3*y +2z² -3z +1
+    i == 2  && return x*(4x*z² -6x*z +2x -2z² +3z -1)
+    i == 3  && return y*(4y*z² -6y*z +2y -2z² +3z -1)
+    i == 4  && return z*(4x²*z -2x² + 8x*y*z -4x*y -6x*z +3x +4y²*z -2y² -6y*z +3y +2z -1)
+    i == 5  && return x*z*(4x*z -2x -2z +1)
+    i == 6  && return y*z*(4y*z -2y -2z +1)
+    i == 7  && return 4x*(-2x*z² +3x*z -x -2*y*z² +3y*z -y +2z² -3z +1)
+    i == 8  && return 4y*(-2x*z² +3x*z -x -2*y*z² +3y*z -y +2z² -3z +1)
+    i == 9  && return 4z*(-2x²*z +2x² -4x*y*z +4x*y +3x*z -3x -2y²*z +2y² +3y*z -3y -z +1)
+    i == 10 && return 4x*y*(2z² -3z +1)
+    i == 11 && return 4x*z*(-2x*z +2x +z -1)
+    i == 12 && return 4y*z*(-2y*z +2y +z -1)
+    i == 13 && return 4x*z*(-2x*z +x -2y*z +y +2z -1)
+    i == 14 && return 4y*z*(-2x*z +x -2y*z +y +2z -1)
+    i == 15 && return 4x*y*z*(2z -1)
+    i == 16 && return 16x*z*(x*z -x +y*z -y -z +1)
+    i == 17 && return 16y*z*(x*z -x +y*z -y -z +1)
+    i == 18 && return 16x*y*z*(1 -z)
+    throw(ArgumentError("no shape function $i for interpolation $ip"))
+end
+
 ###################
 # Bubble elements #
 ###################
@@ -716,18 +924,15 @@ Lagrange element with bubble stabilization.
 struct BubbleEnrichedLagrange{dim,ref_shape,order} <: Interpolation{dim,ref_shape,order} end
 getlowerdim(ip::BubbleEnrichedLagrange{dim,ref_shape,order}) where {dim,ref_shape,order} = Lagrange{dim-1,ref_shape,order}()
 
-# Vertices for all Bubble interpolations are the same
-vertices(::BubbleEnrichedLagrange{2,RefTetrahedron,order}) where {order} = ((1,),(2,),(3,))
-nvertexdofs(::BubbleEnrichedLagrange{ref_dim,ref_shape,order}) where {ref_dim,ref_shape, order} = 1
-
 ################################################
 # Lagrange-Bubble dim 2 RefTetrahedron order 1 #
 ################################################
 # Taken from https://defelement.com/elements/bubble-enriched-lagrange.html
 getnbasefunctions(::BubbleEnrichedLagrange{2,RefTetrahedron,1}) = 4
-ncelldofs(::BubbleEnrichedLagrange{2,RefTetrahedron,1}) = 1
 
-faces(::BubbleEnrichedLagrange{2,RefTetrahedron,1}) = ((1,2), (2,3), (3,1))
+vertexdof_indices(::BubbleEnrichedLagrange{2,RefTetrahedron,1}) = ((1,), (2,), (3,))
+facedof_indices(::BubbleEnrichedLagrange{2,RefTetrahedron,1}) = ((1,2), (2,3), (3,1))
+celldof_interior_indices(::BubbleEnrichedLagrange{2,RefTetrahedron,1}) = (4,)
 
 function reference_coordinates(::BubbleEnrichedLagrange{2,RefTetrahedron,1})
     return [Vec{2, Float64}((1.0, 0.0)),
@@ -752,9 +957,8 @@ end
 struct Serendipity{dim,shape,order} <: Interpolation{dim,shape,order} end
 
 # Vertices for all Serendipity interpolations are the same
-vertices(::Serendipity{2,RefCube,order}) where {order} = ((1,),(2,),(3,),(4,))
-vertices(::Serendipity{3,RefCube,order}) where {order} = ((1,),(2,),(3,),(4,),(5,),(6,),(7,),(8,))
-nvertexdofs(::Serendipity{ref_dim,ref_shape,order}) where {ref_dim,ref_shape, order} = 1
+vertexdof_indices(::Serendipity{2,RefCube,order}) where {order} = ((1,),(2,),(3,),(4,))
+vertexdof_indices(::Serendipity{3,RefCube,order}) where {order} = ((1,),(2,),(3,),(4,),(5,),(6,),(7,),(8,))
 
 #####################################
 # Serendipity dim 2 RefCube order 2 #
@@ -762,9 +966,9 @@ nvertexdofs(::Serendipity{ref_dim,ref_shape,order}) where {ref_dim,ref_shape, or
 getnbasefunctions(::Serendipity{2,RefCube,2}) = 8
 getlowerdim(::Serendipity{2,RefCube,2}) = Lagrange{1,RefCube,2}()
 getlowerorder(::Serendipity{2,RefCube,2}) = Lagrange{2,RefCube,1}()
-nfacedofs(::Serendipity{2,RefCube,2}) = 1
 
-faces(::Serendipity{2,RefCube,2}) = ((1,2,5), (2,3,6), (3,4,7), (4,1,8))
+facedof_indices(::Serendipity{2,RefCube,2}) = ((1,2,5), (2,3,6), (3,4,7), (4,1,8))
+facedof_interior_indices(::Serendipity{2,RefCube,2}) = ((5,), (6,), (7,), (8,))
 
 function reference_coordinates(::Serendipity{2,RefCube,2})
     return [Vec{2, Float64}((-1.0, -1.0)),
@@ -794,12 +998,37 @@ end
 #####################################
 # Serendipity dim 3 RefCube order 2 #
 #####################################
+# Note that second order serendipity hex has no interior face indices.
 getnbasefunctions(::Serendipity{3,RefCube,2}) = 20
 getlowerdim(::Serendipity{3,RefCube,2}) = Serendipity{2,RefCube,2}()
 getlowerorder(::Serendipity{3,RefCube,2}) = Lagrange{3,RefCube,1}()
-nedgedofs(::Serendipity{3,RefCube,2}) = 1
 
-faces(::Serendipity{3,RefCube,2}) = ((1,4,3,2,12,11,10,9), (1,2,6,5,9,18,13,17), (2,3,7,6,10,19,14,18), (3,4,8,7,11,20,15,19), (1,5,8,4,17,16,20,12), (5,6,7,8,13,14,15,16))
+facedof_indices(::Serendipity{3,RefCube,2}) = (
+    (1,4,3,2, 12,11,10,9),
+    (1,2,6,5, 9,18,13,17),
+    (2,3,7,6, 10,19,14,18),
+    (3,4,8,7, 11,20,15,19),
+    (1,5,8,4, 17,16,20,12),
+    (5,6,7,8, 13,14,15,16)
+)
+edgedof_indices(::Serendipity{3,RefCube,2}) = (
+    (1,2, 9),
+    (2,3, 10),
+    (3,4, 11),
+    (4,1, 12),
+    (1,5, 17),
+    (2,6, 18),
+    (3,7, 19),
+    (4,8, 20),
+    (5,6, 13),
+    (6,7, 14),
+    (7,8, 15),
+    (8,5, 16)
+)
+
+edgedof_interior_indices(::Serendipity{3,RefCube,2}) = (
+    (9,), (10,), (11,), (12,), (17), (18,), (19,), (20,), (13,), (14,), (15,), (16,)
+)
 
 function reference_coordinates(::Serendipity{3,RefCube,2})
     return [Vec{3, Float64}((-1.0, -1.0, -1.0)),
@@ -866,8 +1095,9 @@ and Numerical Analysis-Modélisation Mathématique et Analyse Numérique 7.R3 (1
 struct CrouzeixRaviart{dim,order} <: Interpolation{dim,RefTetrahedron,order} end
 
 getnbasefunctions(::CrouzeixRaviart{2,1}) = 3
-nfacedofs(::CrouzeixRaviart{2,1}) = 1
-faces(::CrouzeixRaviart{2,1}) = ((1,), (2,), (3,))
+
+facedof_indices(::CrouzeixRaviart{2,1}) = ((1,), (2,), (3,))
+facedof_interior_indices(::CrouzeixRaviart{2,1}) = ((1,), (2,), (3,))
 
 function reference_coordinates(::CrouzeixRaviart{2,1})
     return [Vec{2, Float64}((0.5, 0.5)),
diff --git a/src/iterators.jl b/src/iterators.jl
index 0cd0c44b68..a37a9f28c5 100644
--- a/src/iterators.jl
+++ b/src/iterators.jl
@@ -57,10 +57,10 @@ function CellCache(grid::Grid{dim,C,T}, flags::UpdateFlags=UpdateFlags()) where
     return CellCache(flags, grid, ScalarWrapper(-1), nodes, coords, nothing, Int[])
 end
 
-function CellCache(dh::Union{DofHandler{dim,T},MixedDofHandler{dim,T}}, flags::UpdateFlags=UpdateFlags()) where {dim,T}
+function CellCache(dh::Union{DofHandler{dim},MixedDofHandler{dim}}, flags::UpdateFlags=UpdateFlags()) where {dim}
     N = nnodes_per_cell(dh.grid)
     nodes = zeros(Int, N)
-    coords = zeros(Vec{dim,T}, N)
+    coords = zeros(Vec{dim, get_coordinate_eltype(dh.grid)}, N)
     n = ndofs_per_cell(dh)
     celldofs = zeros(Int, n)
     return CellCache(flags, dh.grid, ScalarWrapper(-1), nodes, coords, dh, celldofs)
diff --git a/test/test_facevalues.jl b/test/test_facevalues.jl
index eee76a9b0d..7bb0f65038 100644
--- a/test/test_facevalues.jl
+++ b/test/test_facevalues.jl
@@ -21,7 +21,7 @@ for (func_interpol, quad_rule) in  (
         fe_valtype == FaceVectorValues && @test getnbasefunctions(fv) == n_basefuncs * Ferrite.getdim(func_interpol)
 
         xs, n = valid_coordinates_and_normals(func_interpol)
-        for face in 1:getnfaces(func_interpol)
+        for face in 1:Ferrite.nfaces(func_interpol)
             reinit!(fv, xs, face)
             @test Ferrite.getcurrentface(fv) == face
 
@@ -61,7 +61,7 @@ for (func_interpol, quad_rule) in  (
             for i in 1:getnquadpoints(fv)
                 vol += getdetJdV(fv,i)
             end
-            x_face = xs[[Ferrite.faces(func_interpol)[face]...]]
+            x_face = xs[[Ferrite.facedof_indices(func_interpol)[face]...]]
             @test vol ≈ calculate_volume(Ferrite.getlowerdim(func_interpol), x_face)
 
             # Test quadrature rule after reinit! with ref. coords
diff --git a/test/test_interpolations.jl b/test/test_interpolations.jl
index 278b87be9c..e33a0695b8 100644
--- a/test/test_interpolations.jl
+++ b/test/test_interpolations.jl
@@ -1,6 +1,6 @@
 @testset "interpolations" begin
 
-for interpolation in (
+@testset "$interpolation" for interpolation in (
                       Lagrange{1, RefCube, 1}(),
                       Lagrange{1, RefCube, 2}(),
                       Lagrange{2, RefCube, 1}(),
@@ -9,13 +9,15 @@ for interpolation in (
                       Lagrange{2, RefTetrahedron, 2}(),
                       Lagrange{2, RefTetrahedron, 3}(),
                       Lagrange{2, RefTetrahedron, 4}(),
-                      Lagrange{2, RefTetrahedron, 5}(),
+                      #Lagrange{2, RefTetrahedron, 5}(),
                       Lagrange{3, RefCube, 1}(),
                       Lagrange{3, RefCube, 2}(),
                       Serendipity{2, RefCube, 2}(),
                       Serendipity{3, RefCube, 2}(),
                       Lagrange{3, RefTetrahedron, 1}(),
                       Lagrange{3, RefTetrahedron, 2}(),
+                      Lagrange{3, RefPrism, 1}(),
+                      Lagrange{3, RefPrism, 2}(),
                       #
                       DiscontinuousLagrange{1, RefCube, 0}(),
                       DiscontinuousLagrange{2, RefCube, 0}(),
@@ -34,35 +36,79 @@ for interpolation in (
     )
 
     # Test of utility functions
-    ndim = Ferrite.getdim(interpolation)
-    r_shape = Ferrite.getrefshape(interpolation)
+    ref_dim = Ferrite.getdim(interpolation)
+    ref_shape = Ferrite.getrefshape(interpolation)
     func_order = Ferrite.getorder(interpolation)
-    @test typeof(interpolation) <: Interpolation{ndim,r_shape,func_order}
+    @test typeof(interpolation) <: Interpolation{ref_dim,ref_shape,func_order}
 
     # Note that not every element formulation exists for every order and dimension.
-    applicable(Ferrite.getlowerdim, interpolation) && @test typeof(Ferrite.getlowerdim(interpolation)) <: Interpolation{ndim-1}
-    applicable(Ferrite.getlowerorder, interpolation) && @test typeof(Ferrite.getlowerorder(interpolation)) <: Interpolation{ndim,r_shape,func_order-1}
+    applicable(Ferrite.getlowerdim, interpolation) && @test typeof(Ferrite.getlowerdim(interpolation)) <: Interpolation{ref_dim-1}
+    applicable(Ferrite.getlowerorder, interpolation) && @test typeof(Ferrite.getlowerorder(interpolation)) <: Interpolation{ref_dim,ref_shape,func_order-1}
 
+    # Check partition of unity at random point.
     n_basefuncs = getnbasefunctions(interpolation)
-    x = rand(Tensor{1, ndim})
-    f = (x) -> Ferrite.value(interpolation, Tensor{1, ndim}(x))
+    x = rand(Tensor{1, ref_dim})
+    f = (x) -> Ferrite.value(interpolation, Tensor{1, ref_dim}(x))
     @test vec(ForwardDiff.jacobian(f, Array(x))') ≈
            reinterpret(Float64, Ferrite.derivative(interpolation, x))
     @test sum(Ferrite.value(interpolation, x)) ≈ 1.0
 
+    # Check if the important functions are consistens
     coords = Ferrite.reference_coordinates(interpolation)
-    for node in 1:n_basefuncs
-        N_node = Ferrite.value(interpolation, coords[node])
-        for k in 1:node
-            if k == node
-                @test N_node[k] ≈ 1.0
+    @test length(coords) == n_basefuncs
+    @test Ferrite.value(interpolation, length(coords), x) == Ferrite.value(interpolation, length(coords), x)
+    @test_throws ArgumentError Ferrite.value(interpolation, length(coords)+1, x)
+
+    # Test whether we have for each entity corresponding dof indices (possibly empty)
+    @test length(Ferrite.vertexdof_indices(interpolation)) == Ferrite.nvertices(interpolation)
+    if ref_dim > 1
+        @test length(Ferrite.facedof_indices(interpolation)) == Ferrite.nfaces(interpolation)
+        @test length(Ferrite.facedof_interior_indices(interpolation)) == Ferrite.nfaces(interpolation)
+    elseif ref_dim > 2
+        @test length(Ferrite.edgedof_indices(interpolation)) == Ferrite.nedges(interpolation)
+        @test length(Ferrite.edgedof_interior_indices(interpolation)) == Ferrite.nedges(interpolation)
+    end
+    # We have at least as many edge/face dofs as we have edge/face interior dofs
+    if ref_dim > 1
+        @test all(length.(Ferrite.facedof_interior_indices(interpolation)) .<= length.(Ferrite.facedof_indices(interpolation)))
+    elseif ref_dim > 2
+        @test all(length.(Ferrite.edgedof_interior_indices(interpolation)) .<= length.(Ferrite.edgedof_indices(interpolation)))
+    end
+    # The total number of dofs must match the number of base functions
+    totaldofs = sum(length.(Ferrite.vertexdof_indices(interpolation));init=0)
+    if ref_dim > 1
+        totaldofs += sum(length.(Ferrite.facedof_interior_indices(interpolation));init=0)
+    end
+    if ref_dim > 2
+        totaldofs += sum(length.(Ferrite.edgedof_interior_indices(interpolation));init=0) 
+    end
+    totaldofs += length(Ferrite.celldof_interior_indices(interpolation))
+    @test totaldofs == n_basefuncs
+
+    # The dof indices are valid.
+    @test all([all(0 .< i .<= n_basefuncs) for i ∈ Ferrite.vertexdof_indices(interpolation)])
+    if ref_dim > 1
+        @test all([all(0 .< i .<= n_basefuncs) for i ∈ Ferrite.facedof_indices(interpolation)])
+        @test all([all(0 .< i .<= n_basefuncs) for i ∈ Ferrite.facedof_interior_indices(interpolation)])
+    elseif ref_dim > 2
+        @test all([all(0 .< i .<= n_basefuncs) for i ∈ Ferrite.edgedof_indices(interpolation)])
+        @test all([all(0 .< i .<= n_basefuncs) for i ∈ Ferrite.edgedof_interior_indices(interpolation)])
+    end
+    @test all([all(0 .< i .<= n_basefuncs) for i ∈ Ferrite.celldof_interior_indices(interpolation)])
+
+    # Check for dirac delta property of interpolation
+    @testset "dirac delta property of dof $dof" for dof in 1:n_basefuncs
+        N_dof = Ferrite.value(interpolation, coords[dof])
+        for k in 1:n_basefuncs
+            if k == dof
+                @test N_dof[k] ≈ 1.0
             else
-                @test N_node[k] ≈ 0.0 atol=4eps(Float64)
+                @test N_dof[k] ≈ 0.0 atol=4eps(Float64) #broken=typeof(interpolation)==Lagrange{2, RefTetrahedron, 5}&&dof==4&&k==18
             end
         end
     end
 
-    # Test that faces(...) return in counter clockwise order (viewing from the outside)
+    # Test that facedof_indices(...) return in counter clockwise order (viewing from the outside)
     if interpolation isa Lagrange
         function __outward_normal(coords::Vector{<:Vec{1}}, nodes)
             n = coords[nodes[1]]
@@ -81,21 +127,22 @@ for interpolation in (
             n = (p3 - p2) × (p1 - p2)
             return n / norm(n)
         end
-        for (facenodes, normal) in zip(Ferrite.faces(interpolation), reference_normals(interpolation))
+        for (facenodes, normal) in zip(Ferrite.facedof_indices(interpolation), reference_normals(interpolation))
             @test __outward_normal(coords, facenodes) ≈ normal
         end
     end
+
     # regression for https://github.com/Ferrite-FEM/Ferrite.jl/issues/520
     interpolation_type = typeof(interpolation).name.wrapper
     if func_order > 1 && interpolation_type != Ferrite.Serendipity
-        first_order = interpolation_type{ndim,r_shape,1}() 
-        for (highorderface, firstorderface) in zip(Ferrite.faces(interpolation), Ferrite.faces(first_order))
+        first_order = interpolation_type{ref_dim,ref_shape,1}() 
+        for (highorderface, firstorderface) in zip(Ferrite.facedof_indices(interpolation), Ferrite.facedof_indices(first_order))
             for (h_node, f_node) in zip(highorderface, firstorderface)
                 @test h_node == f_node
             end
         end
-        if ndim > 2
-            for (highorderedge, firstorderedge) in zip(Ferrite.edges(interpolation), Ferrite.edges(first_order))
+        if ref_dim > 2
+            for (highorderedge, firstorderedge) in zip(Ferrite.edgedof_indices(interpolation), Ferrite.edgedof_indices(first_order))
                 for (h_node, f_node) in zip(highorderedge, firstorderedge)
                     @test h_node == f_node
                 end
diff --git a/test/test_utils.jl b/test/test_utils.jl
index 00f1327f4f..b64f0f3247 100644
--- a/test/test_utils.jl
+++ b/test/test_utils.jl
@@ -5,6 +5,7 @@
 #####################################
 reference_volume(::Interpolation{dim, RefCube}) where {dim} = 2^dim
 reference_volume(::Interpolation{dim, RefTetrahedron}) where {dim} = 1 / factorial(dim)
+reference_volume(::Interpolation{  3, RefPrism}) = 1/2
 # For faces
 reference_volume(fs::Interpolation, ::Int) = reference_volume(Ferrite.getlowerdim(fs))
 reference_volume(fs::Interpolation{2, RefTetrahedron}, face::Int) = face == 1 ? sqrt(2) : 1.0
@@ -53,6 +54,15 @@ function reference_normals(::Lagrange{3, RefCube})
             Vec{3, Float64}(( 0.0,  0.0,  1.0))]
 end
 
+# Lagrange{3, Wedge}
+function reference_normals(::Lagrange{3, RefPrism})
+    return [Vec{3, Float64}(( 0.0,  0.0, -1.0)),
+            Vec{3, Float64}(( 0.0, -1.0,  0.0)),
+            Vec{3, Float64}((-1.0,  0.0,  0.0)),
+            Vec{3, Float64}((1/√2, 1/√2,  0.0)),
+            Vec{3, Float64}(( 0.0,  0.0,  1.0))]
+end
+
 # Serendipity{2, RefCube}
 reference_normals(::Serendipity{2, RefCube, 2}) = reference_normals(Lagrange{2, RefCube,1}())
 
@@ -161,10 +171,6 @@ function calculate_volume(::Lagrange{0, RefCube, order}, ::Vector{Vec{1, T}}) wh
     return one(T)
 end
 
-getnfaces(::Interpolation{dim, RefCube}) where {dim} = 2*dim
-getnfaces(::Interpolation{2, RefTetrahedron}) = 3
-getnfaces(::Interpolation{3, RefTetrahedron}) = 4
-
 coords_on_faces(x, ::Lagrange{1, RefCube, 1}) = ([x[1]], [x[2]])
 coords_on_faces(x, ::Lagrange{1, RefCube, 2}) = ([x[1]], [x[2]])
 coords_on_faces(x, ::Lagrange{2, RefCube, 1}) =