Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand Down
1 change: 1 addition & 0 deletions docs/src/reference/cellvalues.md
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@ function_value_average
function_gradient_average
function_value_jump
function_gradient_jump
midplane_rotation
getdetJdV_average
get_side_and_baseindex
```
39 changes: 19 additions & 20 deletions src/FerriteInterfaceElements.jl
Original file line number Diff line number Diff line change
@@ -1,40 +1,39 @@
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")
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
106 changes: 84 additions & 22 deletions src/cellvalues.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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,
Expand Down Expand Up @@ -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]

Expand Down Expand Up @@ -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)
Expand Down
86 changes: 60 additions & 26 deletions test/test_cellvalues.jl
Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I could use some help defining a proper test case...
Not sure if the current one is fine.

Original file line number Diff line number Diff line change
@@ -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
Loading