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..779ac2e 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,8 +12,10 @@ 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") include("interpolations.jl") @@ -35,6 +36,7 @@ export function_value_jump, function_gradient_jump, getdetJdV_average, + midplane_rotation, insert_interfaces end diff --git a/src/cellvalues.jl b/src/cellvalues.jl index 532051b..3e4c0e0 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,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) diff --git a/test/test_cellvalues.jl b/test/test_cellvalues.jl index 7f52e17..3186fd4 100644 --- a/test/test_cellvalues.jl +++ b/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 \ No newline at end of file