diff --git a/Project.toml b/Project.toml
index 6e1eacc..2dd2c3c 100644
--- a/Project.toml
+++ b/Project.toml
@@ -6,10 +6,12 @@ version = "1.0.0"
 [deps]
 Ferrite = "c061ca5d-56c9-439f-9c0e-210fe06d3992"
 OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
+StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
 
 [compat]
 Ferrite = "1.0"
 OrderedCollections = "1"
+StaticArraysCore = "1.4.3"
 julia = "1.10"
 
 [extras]
diff --git a/docs/src/reference/cellvalues.md b/docs/src/reference/cellvalues.md
index 5f62539..7464bd8 100644
--- a/docs/src/reference/cellvalues.md
+++ b/docs/src/reference/cellvalues.md
@@ -14,6 +14,7 @@ function_value_average
 function_gradient_average
 function_value_jump
 function_gradient_jump
+midplane_rotation
 getdetJdV_average
 get_side_and_baseindex
 ```
\ No newline at end of file
diff --git a/src/FerriteInterfaceElements.jl b/src/FerriteInterfaceElements.jl
index db4e242..a315887 100644
--- a/src/FerriteInterfaceElements.jl
+++ b/src/FerriteInterfaceElements.jl
@@ -1,11 +1,10 @@
 module FerriteInterfaceElements
 
-import Ferrite
-import Ferrite: AbstractCell, AbstractRefShape, AbstractCellValues,
+import Ferrite: Ferrite, AbstractCell, AbstractRefShape, AbstractCellValues,
     RefLine, RefQuadrilateral, RefTriangle, RefPrism, RefHexahedron,
     Line, QuadraticLine, Triangle, QuadraticTriangle, Quadrilateral, QuadraticQuadrilateral, Tetrahedron, Hexahedron,
     ScalarInterpolation, VectorizedInterpolation, Lagrange,
-    CellValues, QuadratureRule, CellCache, Grid, ExclusiveTopology, FacetIndex, Vec,
+    CellValues, QuadratureRule, CellCache, Grid, ExclusiveTopology, FacetIndex, Vec, Tensor,
     getnbasefunctions, getngeobasefunctions, getorder, n_components, getrefshape,
     vertexdof_indices, edgedof_interior_indices, facedof_interior_indices, volumedof_interior_indices,
     vertices, facets, edges,
@@ -13,28 +12,28 @@ import Ferrite: AbstractCell, AbstractRefShape, AbstractCellValues,
     checkbounds, checkquadpoint, function_value_init, function_gradient_init,
     shape_value_type, shape_gradient_type,
     reinit!, getnquadpoints, getdetJdV, shape_value, shape_gradient, function_value, function_gradient,
-    getcellset, getcells, getneighborhood
+    getcellset, getcells, getneighborhood,
+    ×, norm
 import OrderedCollections: OrderedSet
+import StaticArraysCore: SMatrix, SVector
 
 include("cells.jl")
+export InterfaceCell
+
 include("interpolations.jl")
-include("cellvalues.jl")
-include("grid.jl")
+export InterfaceCellInterpolation
 
-export
-    InterfaceCell,
-    InterfaceCellInterpolation,
-    InterfaceCellValues,
+include("cellvalues.jl")
+export InterfaceCellValues,
     get_side_and_baseindex,
-    shape_value_average,
-    shape_gradient_average,
-    shape_value_jump,
-    shape_gradient_jump,
-    function_value_average,
-    function_gradient_average,
-    function_value_jump,
-    function_gradient_jump,
-    getdetJdV_average,
-    insert_interfaces
+    shape_value_average, shape_gradient_average,
+    shape_value_jump, shape_gradient_jump,
+    function_value_average, function_gradient_average,
+    function_value_jump, function_gradient_jump,
+    midplane_rotation,
+    getdetJdV_average
+
+include("grid.jl")
+export insert_interfaces
 
 end
diff --git a/src/cellvalues.jl b/src/cellvalues.jl
index 532051b..e4eaab4 100644
--- a/src/cellvalues.jl
+++ b/src/cellvalues.jl
@@ -4,6 +4,8 @@
 An `InterfaceCellValues` is based on two `CellValues`: one for each facet of an `InterfaceCell`.
 Since one can use the same `CellValues` for both sides, be default the same object is used for better performance.
 The keyword argument `use_same_cv` can be set to `false` to disable this behavior, if needed.
+Furthermore, the rotation matrix of the midplane in the quadrature points can be computed.
+To do so, set `include_R=true`.
 
 # Fields
 - `ip::InterfaceCellInterpolation`: interpolation on the interface
@@ -12,48 +14,73 @@ The keyword argument `use_same_cv` can be set to `false` to disable this behavio
 - `base_indices_here::Vector{Int}`: base function indices on facet "here"
 - `base_indices_there::Vector{Int}`: base function indices on facet "there"
 - `sides_and_baseindices::Tuple`: side and base function for the base `CellValues` for each base function of the `InterfaceCellValues`
+- `R::Union{AbstractVector,Nothing}`: rotation matrix of the midplane in the quadrature points
 """
-struct InterfaceCellValues{CV} <: AbstractCellValues
+struct InterfaceCellValues{CV,TR} <: AbstractCellValues
     here::CV
     there::CV
     base_indices_here::Vector{Int}
     base_indices_there::Vector{Int}
     sides_and_baseindices::Tuple
+    R::TR # Union{AbstractVector,Nothing} Rotation matrix in quadrature points
 
-    function InterfaceCellValues(ip::IP, here::CV; use_same_cv) where {IP<:InterfaceCellInterpolation, CV<:CellValues}
-        sides_and_baseindices = Tuple( get_side_and_baseindex(ip, i) for i in 1:getnbasefunctions(ip) )
-        base_indices_here  = collect( get_interface_index(ip, :here,  i) for i in 1:getnbasefunctions(ip.base) )
-        base_indices_there = collect( get_interface_index(ip, :there, i) for i in 1:getnbasefunctions(ip.base) )
+    function InterfaceCellValues(ip::IP, here::CV, use_same_cv::Bool, include_R::Val{false}) where {
+            sdim, shape<:AbstractRefShape{sdim}, 
+            IP<:Union{InterfaceCellInterpolation{shape}, VectorizedInterpolation{<:Any,<:Any,<:Any,<:InterfaceCellInterpolation{shape}}},
+            CV<:CellValues}
+        sides_and_baseindices, base_indices_here, base_indices_there = _prepare_baseindices(ip)
         there = use_same_cv ? here : deepcopy(here)
-        return new{CV}(here, there, base_indices_here, base_indices_there, sides_and_baseindices)
+        
+        TR = Nothing
+        R  = nothing
+        return new{CV,TR}(here, there, base_indices_here, base_indices_there, sides_and_baseindices, R)
     end
 
-    function InterfaceCellValues(ip::IP, here::CV; use_same_cv) where {IP<:VectorizedInterpolation{<:Any,<:Any,<:Any,<:InterfaceCellInterpolation}, CV<:CellValues}
-        sides_and_baseindices = Tuple( get_side_and_baseindex(ip, i) for i in 1:getnbasefunctions(ip) )
-        ip = ip.ip
-        base_indices_here  = collect( get_interface_index(ip, :here,  i) for i in 1:getnbasefunctions(ip.base) )
-        base_indices_there = collect( get_interface_index(ip, :there, i) for i in 1:getnbasefunctions(ip.base) )
-        there = use_same_cv ? here : deepcopy(here)
-        return new{CV}(here, there, base_indices_here, base_indices_there, sides_and_baseindices)
-    end
+    function InterfaceCellValues(ip::IP, here::CV, use_same_cv::Bool, include_R::Val{true}) where {
+        sdim, shape<:AbstractRefShape{sdim}, 
+        IP<:Union{InterfaceCellInterpolation{shape}, VectorizedInterpolation{<:Any,<:Any,<:Any,<:InterfaceCellInterpolation{shape}}},
+        CV<:CellValues}
+    sides_and_baseindices, base_indices_here, base_indices_there = _prepare_baseindices(ip)
+    there = use_same_cv ? here : deepcopy(here)
+    
+    T = eltype(here.detJdV)
+    TR = Vector{Tensor{2, sdim, T}}
+    R = TR(undef, getnquadpoints(here))
+    return new{CV,TR}(here, there, base_indices_here, base_indices_there, sides_and_baseindices, R)
+end
+end
+
+function _prepare_baseindices(ip::IP) where {IP<:InterfaceCellInterpolation}
+    sides_and_baseindices = Tuple( get_side_and_baseindex(ip, i) for i in 1:getnbasefunctions(ip) )
+    base_indices_here  = collect( get_interface_index(ip, :here,  i) for i in 1:getnbasefunctions(ip.base) )
+    base_indices_there = collect( get_interface_index(ip, :there, i) for i in 1:getnbasefunctions(ip.base) )
+    return sides_and_baseindices, base_indices_here, base_indices_there
+end
+function _prepare_baseindices(ip::IP) where {IP<:VectorizedInterpolation{<:Any,<:Any,<:Any,<:InterfaceCellInterpolation}}
+    sides_and_baseindices = Tuple( get_side_and_baseindex(ip, i) for i in 1:getnbasefunctions(ip) )
+    ip = ip.ip
+    base_indices_here  = collect( get_interface_index(ip, :here,  i) for i in 1:getnbasefunctions(ip.base) )
+    base_indices_there = collect( get_interface_index(ip, :there, i) for i in 1:getnbasefunctions(ip.base) )
+    return sides_and_baseindices, base_indices_here, base_indices_there
 end
 
 InterfaceCellValues(qr::QuadratureRule, args...; kwargs...) = InterfaceCellValues(Float64, qr, args...; kwargs...)
+InterfaceCellValues(ip, here, use_same_cv, include_R::Bool) = InterfaceCellValues(ip, here, use_same_cv, Val(include_R))
 
 function InterfaceCellValues(::Type{T}, qr::QuadratureRule, 
             ip::InterfaceCellInterpolation, 
             ip_geo::VectorizedInterpolation{sdim,<:Any,<:Any,<:InterfaceCellInterpolation} = default_geometric_interpolation(ip); 
-            use_same_cv=true, kwargs...) where {T, sdim}
+            use_same_cv=true, include_R=Val(false), kwargs...) where {T, sdim}
     cv = CellValues(T, qr, ip.base, VectorizedInterpolation{sdim}(ip_geo.ip.base); kwargs...)
-    return InterfaceCellValues(ip, cv; use_same_cv=use_same_cv)
+    return InterfaceCellValues(ip, cv, use_same_cv, include_R)
 end
 function InterfaceCellValues(::Type{T}, qr::QuadratureRule, 
             ip::VectorizedInterpolation{vdim,<:Any,<:Any,<:InterfaceCellInterpolation}, 
             ip_geo::VectorizedInterpolation{sdim,<:Any,<:Any,<:InterfaceCellInterpolation} = default_geometric_interpolation(ip); 
-            use_same_cv=true, kwargs...) where {T, vdim, sdim}
+            use_same_cv=true, include_R=Val(false), kwargs...) where {T, vdim, sdim}
     cv = CellValues(T, qr, VectorizedInterpolation{vdim}(ip.ip.base), 
                            VectorizedInterpolation{sdim}(ip_geo.ip.base); kwargs...)
-    return InterfaceCellValues(ip, cv; use_same_cv=use_same_cv)
+    return InterfaceCellValues(ip, cv, use_same_cv, include_R)
 end
 
 function InterfaceCellValues(::Type{T}, qr::QuadratureRule, 
@@ -90,13 +117,39 @@ Ferrite.shape_gradient_type(cv::InterfaceCellValues) = shape_gradient_type(cv.he
 
 Ferrite.reinit!(cv::InterfaceCellValues, cc::CellCache) = reinit!(cv, cc.coords) # TODO: Needed?
 
-function Ferrite.reinit!(cv::InterfaceCellValues, x::AbstractVector{Vec{sdim,T}}) where {sdim, T}
-    reinit!(cv.here, @view x[cv.base_indices_here])
-    if cv.here === cv.there
-        reinit!(cv.there, @view x[cv.base_indices_there])
+function Ferrite.reinit!(cv::InterfaceCellValues{CV,TR}, x::AbstractVector{Vec{sdim,T}}) where {sdim, T, CV, TR}
+    x_here  = @view x[cv.base_indices_here]
+    reinit!(cv.here, x_here)
+
+    if ! (cv.here === cv.there)
+        x_there = @view x[cv.base_indices_there]
+        reinit!(cv.there, x_there)
+    end
+
+    if ! isnothing(cv.R)
+        x_there = @view x[cv.base_indices_there]
+        for qp in 1:getnquadpoints(cv.here)
+            mapping_here  = Ferrite.calculate_mapping(cv.here.geo_mapping,  qp, x_here)
+            mapping_there = Ferrite.calculate_mapping(cv.there.geo_mapping, qp, x_there)
+            J_here  = Ferrite.getjacobian(mapping_here)
+            J_there = Ferrite.getjacobian(mapping_there)
+            J = (J_here + J_there) / 2
+            cv.R[qp] = _get_R_from_J(J)
+        end
     end
     return nothing
 end
+function _get_R_from_J(J::SMatrix{2,1,T}) where T 
+    v1 = J[:, 1]
+    v2 = SVector{2,T}((-v1[2], v1[1]))
+    return Tensor{2,2,T}(((v1/norm(v1))..., (v2/norm(v2))...))
+end
+function _get_R_from_J(J::SMatrix{3,2,T}) where T 
+    v1 = J[:, 1]
+    v3 = v1 × J[:, 2]
+    v2 = v3 × v1
+    return Tensor{2,3,T}(((v1/norm(v1))..., (v2/norm(v2))..., (v3/norm(v3))...))
+end
 
 get_side_and_baseindex(cv::InterfaceCellValues, i::Integer) = cv.sides_and_baseindices[i]
 
@@ -277,6 +330,15 @@ function Ferrite.function_gradient_jump(cv::InterfaceCellValues, qp::Int, u::Abs
     return function_gradient(cv, qp, u, false) - function_gradient(cv, qp, u, true)
 end
 
+"""
+    midplane_rotation(cv::InterfaceCellValues, qp::Int)
+
+Compute the rotation matrix of the midplane in quadrature point.
+"""
+function midplane_rotation(cv::InterfaceCellValues, qp::Int)
+    return cv.R[qp]
+end
+
 function Base.show(io::IO, d::MIME"text/plain", cv::InterfaceCellValues)
     if cv.here === cv.there
         print(io, "InterfaceCellValues based on a single CellValues:\n"); show(io, d, cv.here)
diff --git a/test/test_cellvalues.jl b/test/test_cellvalues.jl
index 7f52e17..1157308 100644
--- a/test/test_cellvalues.jl
+++ b/test/test_cellvalues.jl
@@ -1,39 +1,73 @@
 @testset "InterfaceCellValues" begin
-    qr = QuadratureRule{RefTriangle}(1)
-    for (fip, gip) in ( (InterfaceCellInterpolation(Lagrange{RefTriangle, 1}()), InterfaceCellInterpolation(Lagrange{RefTriangle, 1}())),
-                        (InterfaceCellInterpolation(Lagrange{RefTriangle, 2}()), InterfaceCellInterpolation(Lagrange{RefTriangle, 1}())),
-                        (InterfaceCellInterpolation(Lagrange{RefTriangle, 1}()), InterfaceCellInterpolation(Lagrange{RefTriangle, 2}())),
-                        (InterfaceCellInterpolation(Lagrange{RefTriangle, 2}()), InterfaceCellInterpolation(Lagrange{RefTriangle, 2}())),
-                      )
-        @test InterfaceCellValues(qr, fip, gip^3) isa InterfaceCellValues
-        @test InterfaceCellValues(qr, fip^3, gip^3) isa InterfaceCellValues
-        @test InterfaceCellValues(qr, fip^3, gip) isa InterfaceCellValues
+    qr = 
+    for (refshape, dim, iporder, giporder) in ( 
+            (RefTriangle, 3, 1, 1), (RefTriangle, 3, 1, 2), (RefTriangle, 3, 2, 1), (RefTriangle, 3, 2, 2),
+            (RefLine, 2, 1, 1), (RefLine, 2, 1, 2), (RefLine, 2, 2, 1), (RefLine, 2, 2, 2))
+        qr  = QuadratureRule{refshape}(2)
+        fip = InterfaceCellInterpolation(Lagrange{refshape, iporder}())
+        gip = InterfaceCellInterpolation(Lagrange{refshape, giporder}())
+        @test InterfaceCellValues(qr, fip, gip^dim) isa InterfaceCellValues
+        @test InterfaceCellValues(qr, fip^dim, gip^dim) isa InterfaceCellValues
+        @test InterfaceCellValues(qr, fip^dim, gip) isa InterfaceCellValues
         @test InterfaceCellValues(qr, fip, gip) isa InterfaceCellValues
         @test InterfaceCellValues(Float32, qr, fip) isa InterfaceCellValues
         @test InterfaceCellValues(qr, fip) isa InterfaceCellValues
         @test InterfaceCellValues(qr, fip; use_same_cv=false) isa InterfaceCellValues
     end
 
-    ip = InterfaceCellInterpolation(Lagrange{RefTriangle, 1}())
-    cv = InterfaceCellValues(qr, ip)
+    for (refshape, dim) in ((RefTriangle,3), (RefLine,2))
+        qr  = QuadratureRule{refshape}(2)
+        ip = InterfaceCellInterpolation(Lagrange{refshape, 1}())
+        cv = InterfaceCellValues(qr, ip)
+
+        @test getnbasefunctions(cv) == 2*dim
+        @test Ferrite.getngeobasefunctions(cv) == 2*dim
+
+        x = repeat([rand(Vec{dim}) for _ in 1:dim], 2)
+        reinit!(cv, x)
+        nbf = getnbasefunctions(cv)
+        here, there = rand(2)
+        u = vcat(ones(nbf÷2).*here, ones(nbf÷2).*there)
+        for qp in 1:getnquadpoints(cv)
+            @test function_value(cv, qp, u, true) ≈ here
+            @test function_value(cv, qp, u, false) ≈ there
+            @test all(abs.(function_gradient(cv, qp, u, true)) .≤ 1e-14)
+            @test all(abs.(function_gradient(cv, qp, u, false)) .≤ 1e-14)
+            @test function_value_average(cv, qp, u) ≈ (here + there)/2
+            @test function_value_jump(cv, qp, u) ≈ there - here
+            @test all(abs.(function_gradient_average(cv, qp, u)) .≤ 1e-14)
+            @test all(abs.(function_gradient_jump(cv, qp, u)) .≤ 1e-14)
+            @test getdetJdV_average(cv, qp) == (getdetJdV(cv.here, qp) + getdetJdV(cv.there, qp)) / 2
+        end
+    end
 
-    @test getnbasefunctions(cv) == 6
-    @test Ferrite.getngeobasefunctions(cv) == 6
+    qr  = QuadratureRule{RefTriangle}(2)
+    ip = InterfaceCellInterpolation(Lagrange{RefTriangle, 1}())
+    cv = InterfaceCellValues(qr, ip; include_R=true)
+    x  = repeat([Vec{3}((1.0,0.0,0.0)), Vec{3}((0.0,1.0,0.0)), Vec{3}((0.0,0.0,1.0))], 2)
+    n  = Vec{3}(( sqrt(1/3),  sqrt(1/3),   sqrt(1/3)))
+    t₁ = Vec{3}(( sqrt(1/2),  0.0,        -sqrt(1/2)))
+    t₂ = Vec{3}((-sqrt(1/6), 2*sqrt(1/6), -sqrt(1/6)))
+    reinit!(cv, x)
+    for qp in 1:getnquadpoints(cv)
+        R = midplane_rotation(cv, qp)
+        @test tdot(R) ≈ one(R) # Fails e.g. when dx/dξ₁ not perpendicular to dx/ξ₂
+        @test R⋅Vec{3}((1.0,0.0,0.0)) ≈ t₁
+        @test R⋅Vec{3}((0.0,1.0,0.0)) ≈ t₂
+        @test R⋅Vec{3}((0.0,0.0,1.0)) ≈ n
+    end
 
-    x = repeat([rand(Vec{3}), rand(Vec{3}), rand(Vec{3})], 2)
+    qr  = QuadratureRule{RefLine}(2)
+    ip = InterfaceCellInterpolation(Lagrange{RefLine, 1}())
+    cv = InterfaceCellValues(qr, ip; include_R=true)
+    x  = repeat([Vec{2}((1.0,0.0)), Vec{2}((0.0,1.0))], 2)
+    n  = Vec{2}((-sqrt(1/2), -sqrt(1/2)))
+    t  = Vec{2}((-sqrt(1/2),  sqrt(1/2)))
     reinit!(cv, x)
-    nbf = getnbasefunctions(cv)
-    here, there = rand(2)
-    u = vcat(ones(nbf÷2).*here, ones(nbf÷2).*there)
     for qp in 1:getnquadpoints(cv)
-        @test function_value(cv, qp, u, true) ≈ here
-        @test function_value(cv, qp, u, false) ≈ there
-        @test all(abs.(function_gradient(cv, qp, u, true)) .≤ 1e-14)
-        @test all(abs.(function_gradient(cv, qp, u, false)) .≤ 1e-14)
-        @test function_value_average(cv, qp, u) ≈ (here + there)/2
-        @test function_value_jump(cv, qp, u) ≈ there - here
-        @test all(abs.(function_gradient_average(cv, qp, u)) .≤ 1e-14)
-        @test all(abs.(function_gradient_jump(cv, qp, u)) .≤ 1e-14)
-        @test getdetJdV_average(cv, qp) == (getdetJdV(cv.here, qp) + getdetJdV(cv.there, qp)) / 2
+        R = midplane_rotation(cv, qp)
+        @test tdot(R) ≈ one(R) # Fails e.g. when dx/dξ₁ not perpendicular to dx/ξ₂
+        @test R⋅Vec{2}((1.0,0.0)) ≈ t
+        @test R⋅Vec{2}((0.0,1.0)) ≈ n
     end
 end
\ No newline at end of file