Minimize diff from main for including rotation

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]

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
 ```

src/FerriteInterfaceElements.jl

@@ -1,20 +1,21 @@
 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,
     adjust_dofs_during_distribution, default_geometric_interpolation, geometric_interpolation,
     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")
 include("interpolations.jl")
@@ -35,6 +36,7 @@ export
     function_value_jump,
     function_gradient_jump,
     getdetJdV_average,
+    midplane_rotation,
     insert_interfaces
 
 end

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,29 +14,42 @@ 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}
+    function InterfaceCellValues(ip::IP, here::CV; use_same_cv::Bool, include_R::Val) 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) )
         there = use_same_cv ? here : deepcopy(here)
-        return new{CV}(here, there, base_indices_here, base_indices_there, sides_and_baseindices)
+        R = if include_R === Val(false)
+            nothing
+        else
+            T = eltype(here.detJdV)
+            Vector{Tensor{2, Ferrite.getrefdim(ip), T}}(undef, getnquadpoints(here))
+        end
+        return new{CV,typeof(R)}(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}
+    function InterfaceCellValues(ip::IP, here::CV; use_same_cv, include_R::Val) 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)
+        R = if include_R === Val(false)
+            nothing
+        else
+            T = eltype(here.detJdV)
+            Vector{Tensor{2, Ferrite.getrefdim(ip), T}}(undef, getnquadpoints(here))
+        end
+        return new{CV,typeof(R)}(here, there, base_indices_here, base_indices_there, sides_and_baseindices, R)
     end
 end
 
@@ -43,17 +58,17 @@ InterfaceCellValues(qr::QuadratureRule, args...; kwargs...) = InterfaceCellValue
 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=use_same_cv, include_R = isa(include_R, Bool) ? Val(include_R) : 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=use_same_cv, include_R = isa(include_R, Bool) ? Val(include_R) : include_R)
 end
 
 function InterfaceCellValues(::Type{T}, qr::QuadratureRule, 
@@ -90,13 +105,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}, x::AbstractVector{Vec{sdim,T}}) where {sdim, T, CV}
+    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 cv.R !== nothing
+        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 +318,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)

test/test_cellvalues.jl

@@ -36,4 +36,38 @@
         @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
+
+    qr  = QuadratureRule{RefTriangle}(2)
+    ip = InterfaceCellInterpolation(Lagrange{RefTriangle, 1}())
+    cv = InterfaceCellValues(qr, ip; include_R=true)
+    @inferred InterfaceCellValues(qr, ip) # default type stable 
+    @inferred InterfaceCellValues(qr, ip; include_R=Val(false)) # no rotation type stable with Val
+    @inferred InterfaceCellValues(qr, ip; include_R=Val(true)) # with rotation type stable with Val.
+
+    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
+
+    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)
+    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{2}((1.0,0.0)) ≈ t
+        @test R⋅Vec{2}((0.0,1.0)) ≈ n
+    end
 end