"""
    BCValues(func_interpol::Interpolation, geom_interpol::Interpolation, boundary_type::Union{Type{<:BoundaryIndex}})

`BCValues` stores the shape values at all faces/edges/vertices (depending on `boundary_type`) for the geometric interpolation (`geom_interpol`),
for each dof-position determined by the `func_interpol`. Used mainly by the `ConstraintHandler`.
"""
struct BCValues{V_GM, FQR} <: AbstractFaceValues
    geo_mapping::V_GM # AbstractVector{GeometryMapping}
    fqr::FQR          # FaceQuadratureRule
    current_face::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{refshape}, geom_interpol::Interpolation{refshape}, boundary_type::Type{<:BoundaryIndex} = Ferrite.FaceIndex) where {T,dim,refshape <: AbstractRefShape{dim}}
    # 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{refshape,T,dim}[]
    for boundarydofs in dirichlet_boundarydof_indices(boundary_type)(func_interpol)
        dofcoords = Vec{dim,T}[]
        for boundarydof in boundarydofs
            push!(dofcoords, interpolation_coords[boundarydof])
        end
        qrf = QuadratureRule{refshape,T}(fill(T(NaN), length(dofcoords)), dofcoords) # weights will not be used
        push!(qrs, qrf)
    end
    fqr = FaceQuadratureRule(qrs)
    geo_mapping = [GeometryMapping{0}(T, geom_interpol, qr) for qr in qrs]

    BCValues(geo_mapping, fqr, ScalarWrapper(1))
end

nfaces(bcv::BCValues) = length(bcv.geo_mapping)
getcurrentface(bcv) = bcv.current_face[]

function set_current_face!(bcv::BCValues, face_nr::Int)
    # Checking face_nr before setting current_face allows us to use @inbounds 
    # when indexing by getcurrentface(fv) in other places!
    checkbounds(Bool, 1:nfaces(bcv), face_nr) || throw(ArgumentError("Face index out of range."))
    bcv.current_face[] = face_nr
end

getnquadpoints(bcv::BCValues) = @inbounds getnquadpoints(bcv.fqr, getcurrentface(bcv))