diff --git a/Project.toml b/Project.toml
index e3c5eb3fbd..dcdb76f440 100644
--- a/Project.toml
+++ b/Project.toml
@@ -49,6 +49,7 @@ Metis = "2679e427-3c69-5b7f-982b-ece356f1e94b"
 NBInclude = "0db19996-df87-5ea3-a455-e3a50d440464"
 OhMyThreads = "67456a42-1dca-4109-a031-0a68de7e3ad5"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
+QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SHA = "ea8e919c-243c-51af-8825-aaa63cd721ce"
 TaskLocalValues = "ed4db957-447d-4319-bfb6-7fa9ae7ecf34"
@@ -56,4 +57,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 
 [targets]
-test = ["BlockArrays", "Downloads", "FerriteGmsh", "ForwardDiff", "Gmsh", "IterativeSolvers", "Metis", "NBInclude", "OhMyThreads", "ProgressMeter", "Random", "SHA", "TaskLocalValues", "Test", "TimerOutputs", "Logging"]
+test = ["BlockArrays", "Downloads", "FerriteGmsh", "ForwardDiff", "Gmsh", "IterativeSolvers", "Metis", "NBInclude", "OhMyThreads", "ProgressMeter", "Random", "SHA", "TaskLocalValues", "Test", "TimerOutputs", "Logging", "QuadGK"]
diff --git a/src/FEValues/CellValues.jl b/src/FEValues/CellValues.jl
index 18592987d2..18fa2de3a0 100644
--- a/src/FEValues/CellValues.jl
+++ b/src/FEValues/CellValues.jl
@@ -78,6 +78,7 @@ getdetJdV(cv::CellValues, q_point::Int) = cv.detJdV[q_point]
 getdetJdV(::CellValues{<:Any, <:Any, <:Any, Nothing}, ::Int) = throw(ArgumentError("detJdV is not saved in CellValues"))
 
 # Accessors for function values
+get_fun_values(cv::CellValues) = cv.fun_values
 getnbasefunctions(cv::CellValues) = getnbasefunctions(cv.fun_values)
 function_interpolation(cv::CellValues) = function_interpolation(cv.fun_values)
 function_difforder(cv::CellValues) = function_difforder(cv.fun_values)
@@ -111,7 +112,7 @@ function reinit!(cv::CellValues, cell::Union{AbstractCell, Nothing}, x::Abstract
     n_geom_basefuncs = getngeobasefunctions(geo_mapping)
 
     check_reinit_sdim_consistency(:CellValues, shape_gradient_type(cv), eltype(x))
-    if cell === nothing && !isa(mapping_type(fun_values), IdentityMapping)
+    if cell === nothing && reinit_needs_cell(cv)
         throw(ArgumentError("The cell::AbstractCell input is required to reinit! non-identity function mappings"))
     end
     if !checkbounds(Bool, x, 1:n_geom_basefuncs) || length(x) != n_geom_basefuncs
diff --git a/src/FEValues/FacetValues.jl b/src/FEValues/FacetValues.jl
index 9eb84897a3..ff6a44ffc6 100644
--- a/src/FEValues/FacetValues.jl
+++ b/src/FEValues/FacetValues.jl
@@ -135,7 +135,7 @@ function reinit!(fv::FacetValues, cell::Union{AbstractCell, Nothing}, x::Abstrac
     geo_mapping = get_geo_mapping(fv)
     fun_values = get_fun_values(fv)
 
-    if cell === nothing && !isa(mapping_type(fun_values), IdentityMapping)
+    if cell === nothing && reinit_needs_cell(fv)
         throw(ArgumentError("The cell::AbstractCell input is required to reinit! non-identity function mappings"))
     end
 
diff --git a/src/FEValues/FunctionValues.jl b/src/FEValues/FunctionValues.jl
index 006b365789..48397d9aa9 100644
--- a/src/FEValues/FunctionValues.jl
+++ b/src/FEValues/FunctionValues.jl
@@ -143,11 +143,8 @@ end
 
 # Mapping types
 struct IdentityMapping end
-# Not yet implemented:
-# struct CovariantPiolaMapping end # PR798
-# struct ContravariantPiolaMapping end # PR798
-# struct DoubleCovariantPiolaMapping end
-# struct DoubleContravariantPiolaMapping end
+struct CovariantPiolaMapping end
+struct ContravariantPiolaMapping end
 
 mapping_type(fv::FunctionValues) = mapping_type(fv.ip)
 
@@ -159,9 +156,8 @@ the function values and gradients from the reference cell
 to the physical cell geometry.
 """
 required_geo_diff_order(::IdentityMapping, fun_diff_order::Int) = fun_diff_order
-#required_geo_diff_order(::ContravariantPiolaMapping, fun_diff_order::Int) = 1 + fun_diff_order # PR798
-#required_geo_diff_order(::CovariantPiolaMapping,     fun_diff_order::Int) = 1 + fun_diff_order # PR798
-
+required_geo_diff_order(::ContravariantPiolaMapping, fun_diff_order::Int) = 1 + fun_diff_order
+required_geo_diff_order(::CovariantPiolaMapping, fun_diff_order::Int) = 1 + fun_diff_order
 
 # Support for embedded elements
 @inline calculate_Jinv(J::Tensor{2}) = inv(J)
@@ -221,6 +217,57 @@ end
     return nothing
 end
 
-# TODO in PR798, apply_mapping! for
-# * CovariantPiolaMapping
-# * ContravariantPiolaMapping
+# Covariant Piola Mapping
+@inline function apply_mapping!(funvals::FunctionValues{0}, ::CovariantPiolaMapping, q_point::Int, mapping_values, cell)
+    Jinv = inv(getjacobian(mapping_values))
+    @inbounds for j in 1:getnbasefunctions(funvals)
+        d = get_direction(funvals.ip, j, cell)
+        Nξ = funvals.Nξ[j, q_point]
+        funvals.Nx[j, q_point] = d * (Nξ ⋅ Jinv)
+    end
+    return nothing
+end
+
+@inline function apply_mapping!(funvals::FunctionValues{1}, ::CovariantPiolaMapping, q_point::Int, mapping_values, cell)
+    H = gethessian(mapping_values)
+    Jinv = inv(getjacobian(mapping_values))
+    @inbounds for j in 1:getnbasefunctions(funvals)
+        d = get_direction(funvals.ip, j, cell)
+        dNdξ = funvals.dNdξ[j, q_point]
+        Nξ = funvals.Nξ[j, q_point]
+        funvals.Nx[j, q_point] = d * (Nξ ⋅ Jinv)
+        funvals.dNdx[j, q_point] = d * (Jinv' ⋅ dNdξ ⋅ Jinv - Jinv' ⋅ (Nξ ⋅ Jinv ⋅ H ⋅ Jinv))
+    end
+    return nothing
+end
+
+# Contravariant Piola Mapping
+@inline function apply_mapping!(funvals::FunctionValues{0}, ::ContravariantPiolaMapping, q_point::Int, mapping_values, cell)
+    J = getjacobian(mapping_values)
+    detJ = det(J)
+    @inbounds for j in 1:getnbasefunctions(funvals)
+        d = get_direction(funvals.ip, j, cell)
+        Nξ = funvals.Nξ[j, q_point]
+        funvals.Nx[j, q_point] = d * (J ⋅ Nξ) / detJ
+    end
+    return nothing
+end
+
+@inline function apply_mapping!(funvals::FunctionValues{1}, ::ContravariantPiolaMapping, q_point::Int, mapping_values, cell)
+    H = gethessian(mapping_values)
+    J = getjacobian(mapping_values)
+    Jinv = inv(J)
+    detJ = det(J)
+    I2 = one(J)
+    H_Jinv = H ⋅ Jinv
+    A1 = (H_Jinv ⊡ (otimesl(I2, I2))) / detJ
+    A2 = (Jinv' ⊡ H_Jinv) / detJ
+    @inbounds for j in 1:getnbasefunctions(funvals)
+        d = get_direction(funvals.ip, j, cell)
+        dNdξ = funvals.dNdξ[j, q_point]
+        Nξ = funvals.Nξ[j, q_point]
+        funvals.Nx[j, q_point] = d * (J ⋅ Nξ) / detJ
+        funvals.dNdx[j, q_point] = d * (J ⋅ dNdξ ⋅ Jinv / detJ + A1 ⋅ Nξ - (J ⋅ Nξ) ⊗ A2)
+    end
+    return nothing
+end
diff --git a/src/FEValues/common_values.jl b/src/FEValues/common_values.jl
index 402f51cf41..08e8c82b8b 100644
--- a/src/FEValues/common_values.jl
+++ b/src/FEValues/common_values.jl
@@ -9,6 +9,13 @@ function checkquadpoint(fe_v::AbstractValues, qp::Int)
     return nothing
 end
 
+@inline function reinit_needs_cell(fe_values::AbstractValues)
+    # TODO: Might need better logic for this, but for current implementations this
+    # is ok. If someone implements a non-identity mapping that doesn't require the cell
+    # as input, this is only a slight performance issue in some cases.
+    return !isa(mapping_type(get_fun_values(fe_values)), IdentityMapping)
+end
+
 @noinline function throw_incompatible_dof_length(length_ue, n_base_funcs)
     msg = "the number of base functions ($(n_base_funcs)) does not match the length " *
         "of the vector ($(length_ue)). Perhaps you passed the global vector, " *
diff --git a/src/Ferrite.jl b/src/Ferrite.jl
index 4a6e499f53..c1abb64f9d 100644
--- a/src/Ferrite.jl
+++ b/src/Ferrite.jl
@@ -21,7 +21,7 @@ using WriteVTK:
     WriteVTK, VTKCellTypes
 using Tensors:
     Tensors, AbstractTensor, SecondOrderTensor, SymmetricTensor, Tensor, Vec, gradient,
-    rotation_tensor, symmetric, tovoigt!, hessian, otimesu
+    rotation_tensor, symmetric, tovoigt!, hessian, otimesu, otimesl
 using ForwardDiff:
     ForwardDiff
 
diff --git a/src/exports.jl b/src/exports.jl
index c710861c15..786fbf8b29 100644
--- a/src/exports.jl
+++ b/src/exports.jl
@@ -17,6 +17,9 @@ export
     Lagrange,
     DiscontinuousLagrange,
     Serendipity,
+    Nedelec,
+    RaviartThomas,
+    BrezziDouglasMarini,
     getnbasefunctions,
     getrefshape,
 
diff --git a/src/interpolations.jl b/src/interpolations.jl
index 5aaed4fdfb..3810d25291 100644
--- a/src/interpolations.jl
+++ b/src/interpolations.jl
@@ -1774,3 +1774,169 @@ function mapping_type end
 
 mapping_type(::ScalarInterpolation) = IdentityMapping()
 mapping_type(::VectorizedInterpolation) = IdentityMapping()
+
+#####################################
+# RaviartThomas (1st kind), H(div)  #
+#####################################
+struct RaviartThomas{vdim, shape, order} <: VectorInterpolation{vdim, shape, order} end
+mapping_type(::RaviartThomas) = ContravariantPiolaMapping()
+n_dbc_components(::RaviartThomas) = 1
+reference_coordinates(ip::RaviartThomas{vdim}) where {vdim} = fill(NaN * zero(Vec{vdim}), getnbasefunctions(ip))
+dirichlet_edgedof_indices(ip::RaviartThomas{2}) = edgedof_interior_indices(ip)
+dirichlet_facedof_indices(ip::RaviartThomas{3}) = facedof_interior_indices(ip)
+
+# RefTriangle
+edgedof_indices(ip::RaviartThomas{2, RefTriangle}) = edgedof_interior_indices(ip)
+facedof_indices(ip::RaviartThomas{2, RefTriangle}) = (ntuple(i -> i, getnbasefunctions(ip)),)
+
+# RefTriangle, 1st order Lagrange
+# https://defelement.com/elements/examples/triangle-raviart-thomas-lagrange-1.html
+function reference_shape_value(ip::RaviartThomas{2, RefTriangle, 1}, ξ::Vec{2}, i::Int)
+    x, y = ξ
+    i == 1 && return ξ                  # Flip sign
+    i == 2 && return Vec(x - 1, y)      # Keep sign
+    i == 3 && return Vec(x, y - 1)      # Flip sign
+    throw(ArgumentError("no shape function $i for interpolation $ip"))
+end
+
+getnbasefunctions(::RaviartThomas{2, RefTriangle, 1}) = 3
+edgedof_interior_indices(::RaviartThomas{2, RefTriangle, 1}) = ((1,), (2,), (3,))
+facedof_interior_indices(::RaviartThomas{2, RefTriangle, 1}) = ((),)
+adjust_dofs_during_distribution(::RaviartThomas) = false
+
+function get_direction(::RaviartThomas{2, RefTriangle, 1}, j, cell)
+    edge = edges(cell)[j]
+    return ifelse(edge[2] > edge[1], 1, -1)
+end
+
+# RefTriangle, 2st order Lagrange
+# https://defelement.org/elements/examples/triangle-raviart-thomas-lagrange-2.html
+function reference_shape_value(ip::RaviartThomas{2, RefTriangle, 2}, ξ::Vec{2}, i::Int)
+    x, y = ξ
+    # Face 1 (keep ordering, flip sign)
+    i == 1 && return Vec(4x * (2x - 1), 2y * (4x - 1))
+    i == 2 && return Vec(2x * (4y - 1), 4y * (2y - 1))
+    # Face 2 (flip ordering, keep signs)
+    i == 3 && return Vec(8x * y - 2x - 6y + 2, 4y * (2y - 1))
+    i == 4 && return Vec(-8x^2 - 8x * y + 12x + 6y - 4, 2y * (-4x - 4y + 3))
+    # Face 3 (keep ordering, flip sign)
+    i == 5 && return Vec(2x * (3 - 4x - 4y), -8x * y + 6x - 8y^2 + 12y - 4)
+    i == 6 && return Vec(4x * (2x - 1), 8x * y - 6x - 2y + 2)
+    # Cell
+    i == 7 && return Vec(8x * (-2x - y + 2), 8y * (-2x - y + 1))
+    i == 8 && return Vec(8x * (-2y - x + 1), 8y * (-2y - x + 2))
+    throw(ArgumentError("no shape function $i for interpolation $ip"))
+end
+
+getnbasefunctions(::RaviartThomas{2, RefTriangle, 2}) = 8
+edgedof_interior_indices(::RaviartThomas{2, RefTriangle, 2}) = ((1, 2), (3, 4), (5, 6))
+facedof_interior_indices(::RaviartThomas{2, RefTriangle, 2}) = ((7, 8),)
+adjust_dofs_during_distribution(::RaviartThomas{2, RefTriangle, 2}) = true
+
+function get_direction(::RaviartThomas{2, RefTriangle, 2}, j, cell)
+    j > 6 && return 1
+    edge = edges(cell)[(j + 1) ÷ 2]
+    return ifelse(edge[2] > edge[1], 1, -1)
+end
+
+#####################################
+# Brezzi-Douglas–Marini, H(div)     #
+#####################################
+struct BrezziDouglasMarini{vdim, shape, order} <: VectorInterpolation{vdim, shape, order} end
+mapping_type(::BrezziDouglasMarini) = ContravariantPiolaMapping()
+reference_coordinates(ip::BrezziDouglasMarini{vdim}) where {vdim} = fill(NaN * zero(Vec{vdim}), getnbasefunctions(ip))
+dirichlet_edgedof_indices(ip::BrezziDouglasMarini{2}) = edgedof_interior_indices(ip)
+n_dbc_components(::BrezziDouglasMarini) = 1
+
+# RefTriangle
+edgedof_indices(ip::BrezziDouglasMarini{2, RefTriangle}) = edgedof_interior_indices(ip)
+facedof_indices(ip::BrezziDouglasMarini{2, RefTriangle}) = (ntuple(i -> i, getnbasefunctions(ip)),)
+
+# RefTriangle, 1st order Lagrange
+# https://defelement.org/elements/examples/triangle-brezzi-douglas-marini-lagrange-1.html
+function reference_shape_value(ip::BrezziDouglasMarini{2, RefTriangle, 1}, ξ::Vec{2}, i::Int)
+    x, y = ξ
+    # Edge 1: y=1-x, n = [1, 1]/√2 (Flip sign, pos. integration outwards)
+    i == 1 && return Vec(4x, -2y) # N ⋅ n = (2√2 x - √2 (1-x)) = 3√2 x - √2
+    i == 2 && return Vec(-2x, 4y) # N ⋅ n = (-√2x + 2√2 (1-x)) = 2√2 - 3√2x
+    # Edge 2: x=0, n = [-1, 0] (reverse order to follow Ferrite convention)
+    i == 3 && return Vec(-2x - 6y + 2, 4y) # N ⋅ n = (6y - 2)
+    i == 4 && return Vec(4x + 6y - 4, -2y) # N ⋅ n = (4 - 6y)
+    # Edge 3: y=0, n = [0, -1] (Flip sign, pos. integration outwards)
+    i == 5 && return Vec(-2x, 6x + 4y - 4) # N ⋅ n = (4 - 6x)
+    i == 6 && return Vec(4x, -6x - 2y + 2) # N ⋅ n = (6x - 2)
+    throw(ArgumentError("no shape function $i for interpolation $ip"))
+end
+
+getnbasefunctions(::BrezziDouglasMarini{2, RefTriangle, 1}) = 6
+edgedof_interior_indices(::BrezziDouglasMarini{2, RefTriangle, 1}) = ((1, 2), (3, 4), (5, 6))
+adjust_dofs_during_distribution(::BrezziDouglasMarini{2, RefTriangle, 1}) = true
+
+function get_direction(::BrezziDouglasMarini{2, RefTriangle, 1}, j, cell)
+    edge = edges(cell)[(j + 1) ÷ 2]
+    return ifelse(edge[2] > edge[1], 1, -1)
+end
+
+#####################################
+# Nedelec (1st kind), H(curl)       #
+#####################################
+struct Nedelec{vdim, shape, order} <: VectorInterpolation{vdim, shape, order} end
+mapping_type(::Nedelec) = CovariantPiolaMapping()
+reference_coordinates(ip::Nedelec{vdim}) where {vdim} = fill(NaN * zero(Vec{vdim}), getnbasefunctions(ip))
+dirichlet_edgedof_indices(ip::Nedelec) = edgedof_interior_indices(ip)
+n_dbc_components(::Nedelec) = 1
+edgedof_indices(ip::Nedelec) = edgedof_interior_indices(ip)
+
+# 2D refshape (rdim == vdim for Nedelec)
+facedof_indices(ip::Nedelec{2, <:AbstractRefShape{2}}) = (ntuple(i -> i, getnbasefunctions(ip)),)
+
+# RefTriangle, 1st order Lagrange
+# https://defelement.org/elements/examples/triangle-nedelec1-lagrange-1.html
+function reference_shape_value(ip::Nedelec{2, RefTriangle, 1}, ξ::Vec{2}, i::Int)
+    x, y = ξ
+    i == 1 && return Vec(- y, x)
+    i == 2 && return Vec(- y, x - 1) # Changed signed, follow Ferrite's sign convention
+    i == 3 && return Vec(1 - y, x)
+    throw(ArgumentError("no shape function $i for interpolation $ip"))
+end
+
+getnbasefunctions(::Nedelec{2, RefTriangle, 1}) = 3
+edgedof_interior_indices(::Nedelec{2, RefTriangle, 1}) = ((1,), (2,), (3,))
+adjust_dofs_during_distribution(::Nedelec{2, RefTriangle, 1}) = false
+
+function get_direction(::Nedelec{2, RefTriangle, 1}, j, cell)
+    edge = edges(cell)[j]
+    return ifelse(edge[2] > edge[1], 1, -1)
+end
+
+# RefTriangle, 2nd order Lagrange
+# https://defelement.org/elements/examples/triangle-nedelec1-lagrange-2.html
+function reference_shape_value(ip::Nedelec{2, RefTriangle, 2}, ξ::Vec{2}, i::Int)
+    x, y = ξ
+    # Edge 1
+    i == 1 && return Vec(2 * y * (1 - 4 * x), 4 * x * (2 * x - 1))
+    i == 2 && return Vec(4 * y * (1 - 2 * y), 2 * x * (4 * y - 1))
+    # Edge 2 (flip order and sign compared to defelement)
+    i == 3 && return Vec(4 * y * (1 - 2 * y), 8 * x * y - 2 * x - 6 * y + 2)
+    i == 4 && return Vec(2 * y * (4 * x + 4 * y - 3), -8 * x^2 - 8 * x * y + 12 * x + 6 * y - 4)
+    # Edge 3
+    i == 5 && return Vec(8 * x * y - 6 * x + 8 * y^2 - 12 * y + 4, 2 * x * (-4 * x - 4 * y + 3))
+    i == 6 && return Vec(
+        -8 * x * y + 6 * x + 2 * y - 2, 4 * x * (2 * x - 1)
+    )
+    # Face
+    i == 7 && return Vec(8 * y * (-x - 2 * y + 2), 8 * x * (x + 2 * y - 1))
+    i == 8 && return Vec(8 * y * (2 * x + y - 1), 8 * x * (-2 * x - y + 2))
+    throw(ArgumentError("no shape function $i for interpolation $ip"))
+end
+
+getnbasefunctions(::Nedelec{2, RefTriangle, 2}) = 8
+edgedof_interior_indices(::Nedelec{2, RefTriangle, 2}) = ((1, 2), (3, 4), (5, 6))
+facedof_interior_indices(::Nedelec{2, RefTriangle, 2}) = ((7, 8),)
+adjust_dofs_during_distribution(::Nedelec{2, RefTriangle, 2}) = true
+
+function get_direction(::Nedelec{2, RefTriangle, 2}, j, cell)
+    j > 6 && return 1
+    edge = edges(cell)[(j + 1) ÷ 2]
+    return ifelse(edge[2] > edge[1], 1, -1)
+end
diff --git a/src/iterators.jl b/src/iterators.jl
index 84164a7138..f57bbc4084 100644
--- a/src/iterators.jl
+++ b/src/iterators.jl
@@ -84,8 +84,14 @@ function reinit!(cc::CellCache, i::Int)
 end
 
 # reinit! FEValues with CellCache
-reinit!(cv::CellValues, cc::CellCache) = reinit!(cv, cc.coords)
-reinit!(fv::FacetValues, cc::CellCache, f::Int) = reinit!(fv, cc.coords, f)
+function reinit!(cv::CellValues, cc::CellCache)
+    cell = reinit_needs_cell(cv) ? getcells(cc.grid, cellid(cc)) : nothing
+    return reinit!(cv, cell, cc.coords)
+end
+function reinit!(fv::FacetValues, cc::CellCache, f::Int)
+    cell = reinit_needs_cell(fv) ? getcells(cc.grid, cellid(cc)) : nothing
+    return reinit!(fv, cell, cc.coords, f)
+end
 
 # Accessor functions
 getnodes(cc::CellCache) = cc.nodes
diff --git a/test/runtests.jl b/test/runtests.jl
index 57cd82d8a7..88c2586b61 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -11,6 +11,7 @@ using StaticArrays
 using OrderedCollections
 using WriteVTK
 import Metis
+using QuadGK: quadgk
 
 include("test_utils.jl")
 
@@ -39,6 +40,7 @@ include("test_apply_analytical.jl")
 include("PoolAllocator.jl")
 include("test_deprecations.jl")
 include("blockarrays.jl")
+include("test_continuity.jl")
 include("test_examples.jl")
 
 @test all(x -> isdefined(Ferrite, x), names(Ferrite))  # Test that all exported symbols are defined
diff --git a/test/test_cellvalues.jl b/test/test_cellvalues.jl
index e8be79127b..84d412deb6 100644
--- a/test/test_cellvalues.jl
+++ b/test/test_cellvalues.jl
@@ -187,6 +187,69 @@
         @test Ferrite.calculate_mapping(cv2.geo_mapping, 1, cc.coords) == Ferrite.calculate_mapping(ip, ξ, cc.coords, Val(2))
     end
 
+    @testset "Non-identity mapping gradients" begin
+        function test_gradient(ip_fun, cell::CT) where {CT <: Ferrite.AbstractCell}
+            ip_geo = geometric_interpolation(CT)
+            RefShape = Ferrite.getrefshape(ip_fun)
+            x_ref = Ferrite.reference_coordinates(ip_geo)
+            # Random cell coordinates, but small pertubation to ensure 1-1 mapping.
+            cell_coords = (1 .+ rand(length(x_ref)) / 10) .* x_ref
+
+            function calculate_value_differentiable(ξ::Vec; i = 1)
+                T = eltype(ξ)
+                qr = QuadratureRule{RefShape}([zero(T)], [ξ])
+                gm = Ferrite.GeometryMapping{1}(T, ip_geo, qr)
+                mv = Ferrite.calculate_mapping(gm, 1, cell_coords)
+                fv0 = Ferrite.FunctionValues{0}(eltype(T), ip_fun, qr, ip_geo^length(ξ))
+                Ferrite.apply_mapping!(fv0, 1, mv, cell)
+                return shape_value(fv0, 1, i)
+            end
+
+            function spatial_coordinate_differentiable(ξ)
+                x = zero(ξ)
+                for i in 1:getnbasefunctions(ip_geo)
+                    x += cell_coords[i] * Ferrite.reference_shape_value(ip_geo, ξ, i)
+                end
+                return x
+            end
+
+            ξ_rand = sample_random_point(RefShape)
+            qr = QuadratureRule{RefShape}([NaN], [ξ_rand])
+            cv = CellValues(qr, ip_fun, ip_geo)
+            reinit!(cv, cell, cell_coords)
+
+            # Test jacobian calculation
+            dxdξ = gradient(spatial_coordinate_differentiable, ξ_rand) # Jacobian
+            J = Ferrite.getjacobian(Ferrite.calculate_mapping(cv.geo_mapping, 1, cell_coords))
+            @test dxdξ ≈ J
+            for i in 1:getnbasefunctions(ip_fun)
+                dNdξ, N = gradient(z -> calculate_value_differentiable(z; i), ξ_rand, :all)
+                # Test that using FunctionValues{0} and FunctionValues{1} gives same mapped value
+                @test shape_value(cv, 1, i) ≈ N
+                # Test gradient. Hard to differentiate wrt. x, easier to differentiate with ξ,
+                # and use the chain-rule to test.
+                # Note that N̂(ξ) = reference_shape_value != N(ξ) = shape_value in general.
+                #                         dNdx ⋅ dxdξ = dNdξ
+                @test shape_gradient(cv, 1, i) ⋅ dxdξ ≈ dNdξ
+            end
+        end
+        test_ips = [
+            Lagrange{RefTriangle, 2}(), Lagrange{RefQuadrilateral, 2}(), Lagrange{RefHexahedron, 2}()^3, # Test should also work for identity mapping
+            Nedelec{2, RefTriangle, 1}(), Nedelec{2, RefTriangle, 2}(),
+            RaviartThomas{2, RefTriangle, 1}(), RaviartThomas{2, RefTriangle, 2}(), BrezziDouglasMarini{2, RefTriangle, 1}(),
+        ]
+        # Makes most sense to test for nonlinear geometries, so we use such cells for testing only
+        cell_from_refshape(::Type{RefTriangle}) = QuadraticTriangle((ntuple(identity, 6)))
+        cell_from_refshape(::Type{RefQuadrilateral}) = QuadraticQuadrilateral((ntuple(identity, 9)))
+        cell_from_refshape(::Type{RefHexahedron}) = QuadraticHexahedron((ntuple(identity, 27)))
+        for ip in test_ips
+            cell = cell_from_refshape(getrefshape(ip))
+            @testset "$ip" begin
+                test_gradient(ip, cell)
+            end
+        end
+    end
+
     @testset "#265: error message for incompatible geometric interpolation" begin
         dim = 1
         deg = 1
diff --git a/test/test_continuity.jl b/test/test_continuity.jl
new file mode 100644
index 0000000000..f86add5f9a
--- /dev/null
+++ b/test/test_continuity.jl
@@ -0,0 +1,177 @@
+@testset "Field continuity" begin
+    # Integration testing to ensure that correct continuity across facets are obtained
+    # after dofdistribution with use of cell and facet values.
+
+    function find_matching_facet(grid, facet::FacetIndex)
+        cell, facetnr = facet
+        facet_vertices = Set(Ferrite.facets(getcells(grid, cell))[facetnr])
+        for cnr in 1:getncells(grid)
+            cnr == cell && continue
+            for (i, f_vert) in enumerate(Ferrite.facets(getcells(grid, cnr)))
+                facet_vertices == Set(f_vert) && return FacetIndex(cnr, i)
+            end
+        end
+        return nothing
+    end
+
+    function test_continuity(
+            dh::DofHandler, facet::FacetIndex;
+            transformation_function::Function,
+            value_function::Function = function_value
+        )
+        # transformation_function: (v,n) -> z
+        # Examples
+        # * Tangential continuity: fun(v, n) = v - (v ⋅ n)*n
+        # * Normal continuity: fun(v, n) = v ⋅ n
+        # value_function: (fe_v, q_point, ue) -> z
+
+        # Check validity of input
+        @assert length(dh.subdofhandlers) == 1
+        @assert length(Ferrite.getfieldnames(dh)) == 1
+
+        # Find the matching FaceIndex
+        cellnr, facetnr = facet
+        facet2 = find_matching_facet(dh.grid, facet)
+        facet2 === nothing && return false
+
+        # Pick "random" points on the facet
+        cell = getcells(dh.grid, cellnr)
+        RefShape = Ferrite.getrefshape(getcells(dh.grid, cellnr))
+        ip_geo = geometric_interpolation(typeof(cell))
+        ip_fun = Ferrite.getfieldinterpolation(dh, (1, 1))
+        fqr = FacetQuadratureRule{RefShape}(8)
+        fv = FacetValues(fqr, ip_fun, ip_geo)
+        cell_coords = getcoordinates(dh.grid, cellnr)
+        inds = randperm(getnquadpoints(fv))[1:min(4, getnquadpoints(fv))]
+
+        # Random dof vector to test continuity
+        u = rand(ndofs(dh))
+
+        # Calculate coordinates and function values for these
+        point_coords = zeros(eltype(cell_coords), length(inds))
+        point_normal = similar(point_coords)
+        fun_vals = zeros(typeof(shape_value(fv, 1, 1)), length(inds))
+        reinit!(fv, cell, cell_coords, facetnr)
+        ue = u[celldofs(dh, cellnr)]
+        for (i, q_point) in enumerate(inds)
+            point_coords[i] = spatial_coordinate(fv, q_point, cell_coords)
+            point_normal[i] = getnormal(fv, q_point)
+            fun_vals[i] = value_function(fv, q_point, ue)
+        end
+
+        # Calculate function values on the other cell
+        cell2 = getcells(dh.grid, facet2[1])
+        cell_coords2 = getcoordinates(dh.grid, facet2[1])
+        local_coords = map(x -> Ferrite.find_local_coordinate(ip_geo, cell_coords2, x, Ferrite.NewtonLineSearchPointFinder()), point_coords)
+        @assert all(first.(local_coords)) # check that find_local_coordinate converged
+        ξs = collect(last.(local_coords)) # Extract the local coordinates
+        qr = QuadratureRule{RefShape}(zeros(length(ξs)), ξs)
+        cv = CellValues(qr, ip_fun, ip_geo)
+        reinit!(cv, cell2, cell_coords2)
+        # fun_vals2 = similar(fun_vals)
+        ue2 = u[celldofs(dh, facet2[1])]
+        for q_point in 1:getnquadpoints(cv)
+            @assert spatial_coordinate(cv, q_point, cell_coords2) ≈ point_coords[q_point]
+            # Approximate points can contribute to the inaccuracies
+            n = point_normal[q_point]
+            v1 = fun_vals[q_point]
+            v2 = value_function(cv, q_point, ue2)
+            @test isapprox(transformation_function(v1, n), transformation_function(v2, n); rtol = 1.0e-6)
+        end
+
+        #d1 = map((v, n) -> transformation_function(v, n), fun_vals, point_normal)
+        #d2 = map((v, n) -> transformation_function(v, n), fun_vals2, point_normal)
+        # d1 = transformation_function.(fun_vals, point_normal)
+        # d1 = transformation_function.(fun_vals2, point_normal)
+        # @test all(isapprox.(d1, d2; rtol = 1.0e-6))
+        return true
+    end
+
+    tupleshift(t::NTuple{N}, shift::Int) where {N} = ntuple(i -> t[mod(i - 1 - shift, N) + 1], N)
+    #tupleshift(t::NTuple, shift::Int) = tuple(circshift(SVector(t), shift)...)
+    cell_permutations(cell::Quadrilateral) = (Quadrilateral(tupleshift(cell.nodes, shift)) for shift in 0:3)
+    cell_permutations(cell::Triangle) = (Triangle(tupleshift(cell.nodes, shift)) for shift in 0:2)
+    cell_permutations(cell::QuadraticTriangle) = (QuadraticTriangle((tupleshift(cell.nodes[1:3], shift)..., tupleshift(cell.nodes[4:6], shift)...)) for shift in 0:3)
+    cell_permutations(cell::QuadraticQuadrilateral) = (QuadraticQuadrilateral((tupleshift(cell.nodes[1:4], shift)..., tupleshift(cell.nodes[5:8], shift)..., cell.nodes[9])) for shift in 0:4)
+
+    function cell_permutations(cell::Hexahedron)
+        idx = ( #Logic on refshape: Select 1st and 2nd vertex (must be neighbours)
+            # The next follows to create inward vector with RHR, and then 4th is in same plane.
+            # The last four must be the neighbours on the other plane to the first four (same order)
+            (1, 2, 3, 4, 5, 6, 7, 8), (1, 4, 8, 5, 2, 3, 7, 6), (1, 5, 6, 2, 4, 8, 7, 3),
+            (2, 1, 5, 6, 3, 4, 8, 7), (2, 3, 4, 1, 6, 7, 8, 5), (2, 6, 7, 3, 1, 5, 8, 4),
+            (3, 2, 6, 7, 4, 1, 5, 8), (3, 4, 1, 2, 7, 8, 5, 6), (3, 7, 8, 4, 2, 6, 5, 1),
+            (4, 1, 2, 3, 8, 5, 6, 7), (4, 3, 7, 8, 1, 2, 6, 5), (4, 8, 5, 1, 3, 7, 6, 1),
+            (5, 1, 4, 8, 6, 2, 3, 7), (5, 6, 2, 1, 8, 7, 3, 4), (5, 8, 7, 6, 1, 4, 3, 2),
+            (6, 2, 1, 5, 7, 3, 4, 8), (6, 5, 8, 7, 2, 1, 4, 3), (6, 7, 3, 2, 5, 8, 4, 1),
+            (7, 3, 2, 6, 8, 4, 1, 5), (7, 6, 5, 8, 3, 2, 1, 4), (7, 8, 4, 3, 6, 5, 1, 2),
+            (8, 4, 3, 7, 5, 1, 2, 6), (8, 5, 1, 4, 7, 6, 2, 3), (8, 7, 6, 5, 4, 3, 2, 1),
+        )
+        return (Hexahedron(ntuple(i -> cell.nodes[perm[i]], 8)) for perm in idx)
+    end
+
+    function cell_permutations(cell::Tetrahedron)
+        idx = (
+            (1, 2, 3, 4), (1, 3, 4, 2), (1, 4, 2, 3),
+            (2, 1, 4, 3), (2, 3, 1, 4), (2, 4, 3, 1),
+            (3, 1, 2, 4), (3, 2, 4, 1), (3, 4, 1, 2),
+            (4, 1, 3, 2), (4, 3, 2, 1), (4, 2, 1, 3),
+        )
+        return (Tetrahedron(ntuple(i -> cell.nodes[perm[i]], 4)) for perm in idx)
+    end
+
+    continuity_function(ip::VectorizedInterpolation) = continuity_function(ip.ip)
+    continuity_function(::Lagrange) = ((v, _) -> v)
+    continuity_function(::Nedelec) = ((v, n) -> v - n * (v ⋅ n)) # Tangent continuity (H(curl))
+    continuity_function(::RaviartThomas) = ((v, n) -> v ⋅ n) # Normal continuity (H(div))
+    continuity_function(::BrezziDouglasMarini) = ((v, n) -> v ⋅ n) # Normal continuity (H(div))
+
+
+    nel = 3
+
+    cell_types = Dict(
+        RefTriangle => [Triangle, QuadraticTriangle],
+        RefQuadrilateral => [Quadrilateral, QuadraticQuadrilateral],
+        RefTetrahedron => [Tetrahedron],
+        RefHexahedron => [Hexahedron]
+    )
+
+    test_ips = [
+        Lagrange{RefTriangle, 2}(), Lagrange{RefQuadrilateral, 2}(), Lagrange{RefHexahedron, 2}()^3, # Test should also work for identity mapping
+        Nedelec{2, RefTriangle, 1}(), Nedelec{2, RefTriangle, 2}(),
+        RaviartThomas{2, RefTriangle, 1}(), RaviartThomas{2, RefTriangle, 2}(), BrezziDouglasMarini{2, RefTriangle, 1}(),
+    ]
+
+    for ip in test_ips
+        RefShape = getrefshape(ip)
+        dim = Ferrite.getrefdim(ip) # dim = sdim = rdim
+        p1, p2 = (rand(Vec{dim}), ones(Vec{dim}) + rand(Vec{dim}))
+        transfun(x) = typeof(x)(i -> sinpi(x[mod(i, length(x)) + 1] + i / 3)) / 10
+
+        for CT in cell_types[RefShape]
+            grid = generate_grid(CT, ntuple(_ -> nel, dim), p1, p2)
+            # Smoothly distort grid (to avoid spuriously badly deformed elements).
+            # A distorted grid is important to properly test the geometry mapping
+            # for 2nd order elements.
+            transform_coordinates!(grid, x -> (x + transfun(x)))
+            cellnr = getncells(grid) ÷ 2 + 1 # Should be a cell in the center
+            basecell = getcells(grid, cellnr)
+            @testset "$CT, $ip" begin
+                for testcell in cell_permutations(basecell)
+                    grid.cells[cellnr] = testcell
+                    dh = DofHandler(grid)
+                    add!(dh, :u, ip)
+                    close!(dh)
+                    cnt = 0
+                    for facetnr in 1:nfacets(RefShape)
+                        fi = FacetIndex(cellnr, facetnr)
+                        # Check continuity of function value according to continuity_function
+                        found_matching = test_continuity(dh, fi; transformation_function = continuity_function(ip))
+                        cnt += found_matching
+                    end
+                    @assert cnt > 0
+                end
+            end
+        end
+    end
+end
diff --git a/test/test_interpolations.jl b/test/test_interpolations.jl
index 3b4a497e1c..06c5158fd7 100644
--- a/test/test_interpolations.jl
+++ b/test/test_interpolations.jl
@@ -265,5 +265,81 @@ using Ferrite: reference_shape_value, reference_shape_gradient
         @test_throws ArgumentError Ferrite.facetdof_interior_indices(ip)
     end
 
+    @testset "H(curl) and H(div)" begin
+        Hcurl_interpolations = [Nedelec{2, RefTriangle, 1}(), Nedelec{2, RefTriangle, 2}()] # Nedelec{3, RefTetrahedron, 1}(), Nedelec{3, RefHexahedron, 1}()]
+        Hdiv_interpolations = [RaviartThomas{2, RefTriangle, 1}(), RaviartThomas{2, RefTriangle, 2}(), BrezziDouglasMarini{2, RefTriangle, 1}()]
+        # test_interpolation_properties.(Hcurl_interpolations)  # Requires PR1136
+        # test_interpolation_properties.(Hdiv_interpolations)   # Requires PR1136
+
+        # These reference moments define the functionals that an interpolation should fulfill
+        reference_moment(::RaviartThomas{2, RefTriangle, 1}, s, facet_shape_nr) = 1
+        reference_moment(::RaviartThomas{2, RefTriangle, 2}, s, facet_shape_nr) = facet_shape_nr == 1 ? (1 - s) : s
+        reference_moment(::BrezziDouglasMarini{2, RefTriangle, 1}, s, facet_shape_nr) = facet_shape_nr == 1 ? (1 - s) : s
+        reference_moment(::Nedelec{2, RefTriangle, 1}, s, edge_shape_nr) = 1
+        reference_moment(::Nedelec{2, RefTriangle, 2}, s, edge_shape_nr) = edge_shape_nr == 1 ? (1 - s) : s
+
+        function_space(::RaviartThomas) = Val(:Hdiv)
+        function_space(::BrezziDouglasMarini) = Val(:Hdiv)
+        function_space(::Nedelec) = Val(:Hcurl)
+        # function_space(::Lagrange) = Val(:H1)
+        # function_space(::DiscontinuousLagrange) = Val(:L2)
+
+        function test_interpolation_functionals(ip::Interpolation)
+            return test_interpolation_functionals(function_space(ip), Val(Ferrite.getrefdim(ip)), ip)
+        end
+
+        # 2D, H(div) -> facet
+        function test_interpolation_functionals(::Val{:Hdiv}, ::Val{2}, ip::Interpolation)
+            RefShape = getrefshape(ip)
+            ipg = Lagrange{RefShape, 1}()
+            for facetnr in 1:nfacets(RefShape)
+                facet_coords = getindex.((Ferrite.reference_coordinates(ipg),), Ferrite.reference_facets(RefShape)[facetnr])
+                dof_inds = Ferrite.facetdof_interior_indices(ip)[facetnr]
+                ξ(s) = facet_coords[1] + (facet_coords[2] - facet_coords[1]) * s
+                weighted_normal = norm(facet_coords[2] - facet_coords[1]) * reference_normals(ipg)[facetnr]
+                for (facet_shape_nr, shape_nr) in pairs(dof_inds)
+                    moment_fun(s) = reference_moment(ip, s, facet_shape_nr)
+                    f(s) = moment_fun(s) * (reference_shape_value(ip, ξ(s), shape_nr) ⋅ weighted_normal)
+                    val, _ = quadgk(f, 0, 1; atol = 1.0e-8)
+                    @test val ≈ 1
+                end
+
+                for shape_nr in 1:getnbasefunctions(ip)
+                    shape_nr in dof_inds && continue # Already tested
+                    for s in range(0, 1, 5)
+                        @test abs(reference_shape_value(ip, ξ(s), shape_nr) ⋅ weighted_normal) < 1.0e-10
+                    end
+                end
+            end
+        end
+
+        function test_interpolation_functionals(::Val{:Hcurl}, ::Val{2}, ip::Interpolation)
+            RefShape = getrefshape(ip)
+            ipg = Lagrange{RefShape, 1}()
+            for edgenr in 1:Ferrite.nedges(RefShape)
+                edge_coords = getindex.((Ferrite.reference_coordinates(ipg),), Ferrite.reference_edges(RefShape)[edgenr])
+                edge_vector = edge_coords[2] - edge_coords[1] # Weighted direction vector
+                dof_inds = Ferrite.edgedof_interior_indices(ip)[edgenr]
+                ξ(s) = edge_coords[1] + edge_vector * s
+                for (edge_shape_nr, shape_nr) in pairs(dof_inds)
+                    moment_fun(s) = reference_moment(ip, s, edge_shape_nr)
+                    f(s) = moment_fun(s) * (reference_shape_value(ip, ξ(s), shape_nr) ⋅ edge_vector)
+                    val, _ = quadgk(f, 0, 1; atol = 1.0e-8)
+                    @test val ≈ 1
+                end
+
+                for shape_nr in 1:getnbasefunctions(ip)
+                    shape_nr in dof_inds && continue # Already tested
+                    for s in range(0, 1, 5)
+                        @test abs(reference_shape_value(ip, ξ(s), shape_nr) ⋅ edge_vector) < 1.0e-10
+                    end
+                end
+            end
+        end
+
+        test_interpolation_functionals.(Hdiv_interpolations)
+        test_interpolation_functionals.(Hcurl_interpolations)
+    end
+
 
 end # testset