diff --git a/src/interpolations.jl b/src/interpolations.jl
index 7bbda34b9c..e168f5c510 100644
--- a/src/interpolations.jl
+++ b/src/interpolations.jl
@@ -1839,6 +1839,64 @@ function get_direction(::RaviartThomas{RefTriangle, 2}, j, cell)
     return ifelse(edge[2] > edge[1], 1, -1)
 end
 
+# RefTetrahedron, 1st order Lagrange
+# https://defelement.com/elements/examples/tetrahedron-raviart-thomas-lagrange-1.html
+function reference_shape_value(ip::RaviartThomas{RefTetrahedron, 1}, ξ::Vec{3}, i::Int)
+    x, y, z = ξ
+    i == 1 && return Vec(2 * x, 2 * y, 2 * (z - 1)) # Changed sign, follow Ferrite's sign convention
+    i == 2 && return Vec(2 * x, 2 * (y - 1), 2 * z)
+    i == 3 && return Vec(2 * x, 2 * y, 2 * z)
+    i == 4 && return Vec(2 * (x - 1), 2 * y, 2 * z) # Changed sign, follow Ferrite's sign convention
+    throw(ArgumentError("no shape function $i for interpolation $ip"))
+end
+
+getnbasefunctions(::RaviartThomas{RefTetrahedron, 1}) = 4
+edgedof_indices(ip::RaviartThomas{RefTetrahedron, 1}) = edgedof_interior_indices(ip)
+facedof_indices(ip::RaviartThomas{RefTetrahedron, 1}) = facedof_interior_indices(ip)
+edgedof_interior_indices(::RaviartThomas{RefTetrahedron, 1}) = ((), (), (), (), (), ())
+facedof_interior_indices(::RaviartThomas{RefTetrahedron, 1}) = ((1), (2), (3), (4))
+adjust_dofs_during_distribution(::RaviartThomas{RefTetrahedron, 1}) = false
+
+function get_direction(::RaviartThomas{RefTetrahedron, 1}, j, cell)
+    face = faces(cell)[j]
+    _, idx_min = findmin(face)
+
+    idx_next = mod(idx_min, length(face)) + 1
+    idx_prev = mod(idx_min - 2, length(face)) + 1
+    return face[idx_next] < face[idx_prev] ? 1 : -1
+end
+
+# RefHexahedron, 1st order Lagrange
+# https://defelement.com/elements/examples/hexahedron-raviart-thomas-lagrange-1.html
+function reference_shape_value(ip::RaviartThomas{RefHexahedron, 1}, ξ::Vec{3,T}, i::Int) where {T}
+    x, y, z = ξ
+    nil = zero(T)
+
+    i == 1 && return Vec(nil, nil, (z - 1) / 2) # Changed sign, follow Ferrite's sign convention
+    i == 2 && return Vec(nil, (y - 1) / 2, nil)
+    i == 3 && return Vec((x + 1) / 2, nil, nil)
+    i == 4 && return Vec(nil, (y + 1) / 2, nil) # Changed sign, follow Ferrite's sign convention
+    i == 5 && return Vec((x - 1) / 2, nil, nil) # Changed sign, follow Ferrite's sign convention
+    i == 6 && return Vec(nil, nil, (z + 1) / 2)
+    throw(ArgumentError("no shape function $i for interpolation $ip"))
+end
+
+getnbasefunctions(::RaviartThomas{RefHexahedron, 1}) = 6
+edgedof_indices(ip::RaviartThomas{RefHexahedron, 1}) = edgedof_interior_indices(ip)
+facedof_indices(ip::RaviartThomas{RefHexahedron, 1}) = facedof_interior_indices(ip)
+edgedof_interior_indices(::RaviartThomas{RefHexahedron, 1}) = ((), (), (), (), (), (), (), (), (), (), (), ())
+facedof_interior_indices(::RaviartThomas{RefHexahedron, 1}) = ((1), (2), (3), (4), (5), (6))
+adjust_dofs_during_distribution(::RaviartThomas{RefHexahedron, 1}) = false
+
+function get_direction(::RaviartThomas{RefHexahedron, 1}, j, cell)
+    face = faces(cell)[j]
+    _, idx_min = findmin(face)
+
+    idx_next = mod(idx_min, length(face)) + 1
+    idx_prev = mod(idx_min - 2, length(face)) + 1
+    return face[idx_next] < face[idx_prev] ? 1 : -1
+end
+
 #####################################
 # Brezzi-Douglas–Marini, H(div)     #
 #####################################
@@ -1942,3 +2000,59 @@ function get_direction(::Nedelec{RefTriangle, 2}, j, cell)
     edge = edges(cell)[(j + 1) ÷ 2]
     return ifelse(edge[2] > edge[1], 1, -1)
 end
+
+# RefTetrahedron, 1st order Lagrange
+# https://defelement.org/elements/examples/tetrahedron-nedelec1-lagrange-1.html
+function reference_shape_value(ip::Nedelec{RefTetrahedron, 1}, ξ::Vec{3,T}, i::Int) where {T}
+    x, y, z = ξ
+    nil = zero(T)
+
+    i == 1 && return Vec(1 - y - z, x, x)
+    i == 2 && return Vec(-y, x, nil)
+    i == 3 && return Vec(-y, x + z - 1, -y) # Changed sign, follow Ferrite's sign convention
+    i == 4 && return Vec(z, z, 1 - x - y)
+    i == 5 && return Vec(-z, nil, x)
+    i == 6 && return Vec(nil, -z, y)
+    throw(ArgumentError("no shape function $i for interpolation $ip"))
+end
+
+getnbasefunctions(::Nedelec{RefTetrahedron, 1}) = 6
+edgedof_interior_indices(::Nedelec{RefTetrahedron, 1}) = ((1,), (2,), (3,), (4,), (5,), (6,))
+facedof_indices(::Nedelec{RefTetrahedron, 1}) = ((1, 2, 3), (1, 4, 5), (2, 5, 6), (3, 4, 6))
+adjust_dofs_during_distribution(::Nedelec{RefTetrahedron, 1}) = false
+
+function get_direction(::Nedelec{RefTetrahedron, 1}, j, cell)
+    edge = edges(cell)[j]
+    return ifelse(edge[2] > edge[1], 1, -1)
+end
+
+# RefHexahedron, 1st order Lagrange
+# https://defelement.org/elements/examples/hexahedron-nedelec1-lagrange-1.html
+function reference_shape_value(ip::Nedelec{RefHexahedron, 1}, ξ::Vec{3,T}, i::Int) where {T}
+    x, y, z = ξ
+    nil = zero(T)
+
+    i == 1 && return Vec((y * z - y - z + 1) / 8, nil, nil)
+    i == 2 && return Vec(nil, -(x * z - x + z - 1) / 8, nil)
+    i == 3 && return Vec((y * z - y + z - 1) / 8, nil, nil) # Changed sign, follow Ferrite's sign convention
+    i == 4 && return Vec(nil, -(x * z - x - z + 1) / 8, nil)  # Changed sign, follow Ferrite's sign convention
+    i == 5 && return Vec(-(z * y - z + y - 1) / 8, nil, nil)
+    i == 6 && return Vec(nil, (x * z + x + z + 1) / 8, nil)
+    i == 7 && return Vec(-(y * z + y + z + 1) / 8, nil, nil)  # Changed sign, follow Ferrite's sign convention
+    i == 8 && return Vec(nil, (z * x - z + x - 1) / 8, nil) # Changed sign, follow Ferrite's sign convention
+    i == 9 && return Vec(nil, nil, (x * y - x - y + 1) / 8)
+    i == 10 && return Vec(nil, nil, -(x * y - x + y - 1) / 8)
+    i == 11 && return Vec(nil, nil, (x * y + x + y + 1) / 8)
+    i == 12 && return Vec(nil, nil, -(y * x - y + x - 1) / 8) 
+    throw(ArgumentError("no shape function $i for interpolation $ip"))
+end
+
+getnbasefunctions(::Nedelec{RefHexahedron, 1}) = 12
+edgedof_interior_indices(::Nedelec{RefHexahedron, 1}) = ((1,), (2,), (3,), (4,), (5,), (6,), (7,), (8,), (9,), (10,), (11,), (12,))
+facedof_indices(::Nedelec{RefHexahedron, 1}) = ((1, 2, 3, 4), (1, 5, 9, 10), (2, 6, 10, 11), (3, 7, 11, 12), (4, 8, 9, 12), (5, 6, 7, 8))
+adjust_dofs_during_distribution(::Nedelec{RefHexahedron, 1}) = false
+
+function get_direction(::Nedelec{RefHexahedron, 1}, j, cell)
+    edge = edges(cell)[j]
+    return ifelse(edge[2] > edge[1], 1, -1)
+end
diff --git a/test/test_continuity.jl b/test/test_continuity.jl
index 009c659a3e..d611712077 100644
--- a/test/test_continuity.jl
+++ b/test/test_continuity.jl
@@ -130,8 +130,9 @@
 
     test_ips = [
         Lagrange{RefTriangle, 2}(), Lagrange{RefQuadrilateral, 2}(), Lagrange{RefHexahedron, 2}()^3, # Test should also work for identity mapping
-        Nedelec{RefTriangle, 1}(), Nedelec{RefTriangle, 2}(),
-        RaviartThomas{RefTriangle, 1}(), RaviartThomas{RefTriangle, 2}(), BrezziDouglasMarini{RefTriangle, 1}(),
+        Nedelec{RefTriangle, 1}(), Nedelec{RefTriangle, 2}(), Nedelec{RefTetrahedron, 1}(), Nedelec{RefHexahedron, 1}(),
+        RaviartThomas{RefTriangle, 1}(), RaviartThomas{RefTriangle, 2}(), RaviartThomas{RefTetrahedron, 1}(), RaviartThomas{RefHexahedron, 1}(),
+        BrezziDouglasMarini{RefTriangle, 1}(),
     ]
 
     for ip in test_ips
diff --git a/test/test_interpolations.jl b/test/test_interpolations.jl
index 5cb2d428f1..cfe47c00a9 100644
--- a/test/test_interpolations.jl
+++ b/test/test_interpolations.jl
@@ -351,17 +351,28 @@ end
     end
 
     @testset "H(curl) and H(div)" begin
-        Hcurl_interpolations = [Nedelec{RefTriangle, 1}(), Nedelec{RefTriangle, 2}()] # Nedelec{3, RefTetrahedron, 1}(), Nedelec{3, RefHexahedron, 1}()]
-        Hdiv_interpolations = [RaviartThomas{RefTriangle, 1}(), RaviartThomas{RefTriangle, 2}(), BrezziDouglasMarini{RefTriangle, 1}()]
+        Hcurl_interpolations = [
+            Nedelec{RefTriangle, 1}(), Nedelec{RefTriangle, 2}(), 
+            Nedelec{RefTetrahedron, 1}(), Nedelec{RefHexahedron, 1}()
+        ]
+        Hdiv_interpolations = [
+            RaviartThomas{RefTriangle, 1}(), RaviartThomas{RefTriangle, 2}(), 
+            RaviartThomas{RefTetrahedron, 1}(), RaviartThomas{RefHexahedron, 1}(),
+            BrezziDouglasMarini{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{RefTriangle, 1}, s, facet_shape_nr) = 1
         reference_moment(::RaviartThomas{RefTriangle, 2}, s, facet_shape_nr) = facet_shape_nr == 1 ? (1 - s) : s
+        reference_moment(::RaviartThomas{RefTetrahedron, 1}, s0, s1, facet_shape_nr) = 1
+        reference_moment(::RaviartThomas{RefHexahedron, 1}, s0, s1, facet_shape_nr) = 1
         reference_moment(::BrezziDouglasMarini{RefTriangle, 1}, s, facet_shape_nr) = facet_shape_nr == 1 ? (1 - s) : s
         reference_moment(::Nedelec{RefTriangle, 1}, s, edge_shape_nr) = 1
         reference_moment(::Nedelec{RefTriangle, 2}, s, edge_shape_nr) = edge_shape_nr == 1 ? (1 - s) : s
+        reference_moment(::Nedelec{RefTetrahedron, 1}, s, edge_shape_nr) = 1
+        reference_moment(::Nedelec{RefHexahedron, 1}, s, edge_shape_nr) = 1
 
         function_space(::RaviartThomas) = Val(:Hdiv)
         function_space(::BrezziDouglasMarini) = Val(:Hdiv)
@@ -373,7 +384,7 @@ end
             return test_interpolation_functionals(function_space(ip), Val(Ferrite.getrefdim(ip)), ip)
         end
 
-        # 2D, H(div) -> facet
+        # 2D, H(div)
         function test_interpolation_functionals(::Val{:Hdiv}, ::Val{2}, ip::Interpolation)
             RefShape = getrefshape(ip)
             ipg = Lagrange{RefShape, 1}()
@@ -398,7 +409,34 @@ end
             end
         end
 
-        function test_interpolation_functionals(::Val{:Hcurl}, ::Val{2}, ip::Interpolation)
+        # 3D, H(div)
+        function test_interpolation_functionals(::Val{:Hdiv}, ::Val{3}, 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]
+                ξ(s0, s1) = facet_coords[1] + (facet_coords[2] - facet_coords[1]) * s0 + (facet_coords[3] - facet_coords[1]) * s1
+                weighted_normal = reference_normals(RefShape)[facetnr] * reference_face_area(ip, facetnr)
+                for (facet_shape_nr, shape_nr) in pairs(dof_inds)
+                    moment_fun(s0, s1) = reference_moment(ip, s0, s1, facet_shape_nr)
+                    f(s0, s1) = moment_fun(s0, s1) * (reference_shape_value(ip, ξ(s0, s1), shape_nr) ⋅ weighted_normal)
+
+                    if(length(facet_coords) == 3)
+                        val, _ = quadgk(s0 -> quadgk(s1 -> f(s0, s1), 0, 1 - s0; atol = 1.0e-8)[1], 0, 1; atol = 1.0e-8) # TODO Replace quadgk by more suitable 2D cubature
+                        @test val ≈ 0.5
+                    elseif(length(facet_coords) == 4)
+                        val, _ = quadgk(s0 -> quadgk(s1 -> f(s0, s1), 0, 1; atol = 1.0e-8)[1], 0, 1; atol = 1.0e-8) # TODO Replace quadgk by more suitable 2D cubature
+                        @test val ≈ 4
+                    else
+                        error("Cubature not defined for facets that are not triangles or quadrilaterals")
+                    end             
+                end
+            end
+        end
+
+        # 2D and 3D, H(curl)
+        function test_interpolation_functionals(::Val{:Hcurl}, ::Union{Val{2}, Val{3}}, ip::Interpolation)
             RefShape = getrefshape(ip)
             ipg = Lagrange{RefShape, 1}()
             for edgenr in 1:Ferrite.nedges(RefShape)