Skip to content

Commit

Permalink
Use Tensor for interface transformation matrices (#1100)
Browse files Browse the repository at this point in the history 
* Fix type stability issue of InterfaceValues
* SMatrix -> Tensor internally for InterfaceValues
* Regression test coverage
  • Loading branch information
@AbdAlazezAhmed
AbdAlazezAhmed authored Nov 26, 2024
1 parent aa35522 commit e7a671d
Show file tree
Hide file tree
Showing 3 changed files with 84 additions and 29 deletions.
20 changes: 10 additions & 10 deletions src/FEValues/InterfaceValues.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -454,16 +454,16 @@ end
θ = 2 * shift_index / 3
θpre = 2 * lowest_node_shift_index / 3

flipping = SMatrix{3, 3}(1.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, 1.0)
flipping = Tensor{2, 3}((1.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, 1.0))

translate_1 = SMatrix{3, 3}(1.0, 0.0, 0.0, 0.0, 1.0, 0.0, -sinpi(2 / 3) / 3, -0.5, 1.0)
stretch_1 = SMatrix{3, 3}(sinpi(2 / 3), 0.5, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0)
translate_1 = Tensor{2, 3}((1.0, 0.0, 0.0, 0.0, 1.0, 0.0, -sinpi(2 / 3) / 3, -0.5, 1.0))
stretch_1 = Tensor{2, 3}((sinpi(2 / 3), 0.5, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0))

translate_2 = SMatrix{3, 3}(1.0, 0.0, 0.0, 0.0, 1.0, 0.0, sinpi(2 / 3) / 3, 0.5, 1.0)
stretch_2 = SMatrix{3, 3}(1 / sinpi(2 / 3), -1 / 2 / sinpi(2 / 3), 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0)
translate_2 = Tensor{2, 3}((1.0, 0.0, 0.0, 0.0, 1.0, 0.0, sinpi(2 / 3) / 3, 0.5, 1.0))
stretch_2 = Tensor{2, 3}((1 / sinpi(2 / 3), -1 / 2 / sinpi(2 / 3), 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0))

return flipped ? stretch_2 * translate_2 * rotation_tensor(0, 0, θpre * pi) * flipping * rotation_tensor(0, 0, (θ - θpre) * pi) * translate_1 * stretch_1 :
stretch_2 * translate_2 * rotation_tensor(0, 0, θ * pi) * translate_1 * stretch_1
return flipped ? stretch_2 ⋅ translate_2 ⋅ rotation_tensor(0, 0, θpre * pi) ⋅ flipping ⋅ rotation_tensor(0, 0, (θ - θpre) * pi) ⋅ translate_1 ⋅ stretch_1 :
stretch_2 ⋅ translate_2 ⋅ rotation_tensor(0, 0, θ * pi) ⋅ translate_1 ⋅ stretch_1
end

@inline function _get_transformation_matrix(::NTuple{4, Int}, interface_transformation::InterfaceOrientationInfo)
Expand All @@ -474,8 +474,8 @@ end
θ = 2 * shift_index / 4
θpre = 2 * lowest_node_shift_index / 4

flipping = SMatrix{3, 3}(0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0)
return flipped ? rotation_tensor(0, 0, θpre * pi) * flipping * rotation_tensor(0, 0, (θ - θpre) * pi) : rotation_tensor(0, 0, θ * pi)
flipping = Tensor{2, 3}((0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0))
return flipped ? rotation_tensor(0, 0, θpre * pi) ⋅ flipping ⋅ rotation_tensor(0, 0, (θ - θpre) * pi) : rotation_tensor(0, 0, θ * pi)
end

@inline function _get_transformation_matrix(::NTuple{N, Int}, ::InterfaceOrientationInfo) where {N}
Expand Down Expand Up @@ -553,7 +553,7 @@ function transform_interface_points!(dst::AbstractVector{Vec{3, Float64}}, point
M = get_transformation_matrix(interface_transformation)
for (idx, point) in pairs(points)
face_point = element_to_facet_transformation(point, RefShapeA, facet_a)
result = M * Vec(face_point[1], face_point[2], 1.0)
result = M ⋅ Vec(face_point[1], face_point[2], 1.0)
dst[idx] = facet_to_element_transformation(Vec(result[1], result[2]), RefShapeB, facet_b)
end
return nothing
Expand Down
4 changes: 2 additions & 2 deletions src/iterators.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -341,10 +341,10 @@ end
`InterfaceIterator` is stateful and should not be used for things other than `for`-looping
(e.g. broadcasting over, or collecting the iterator may yield unexpected results).
"""
struct InterfaceIterator{IC <: InterfaceCache, G <: Grid}
struct InterfaceIterator{IC <: InterfaceCache, G <: AbstractGrid, TopologyType <: AbstractTopology}
cache::IC
grid::G
topology::ExclusiveTopology
topology::TopologyType
end

function InterfaceIterator(
Expand Down
89 changes: 72 additions & 17 deletions test/test_interfacevalues.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -200,23 +200,78 @@
# DiscontinuousLagrange{RefTriangle, 1}(), FacetQuadratureRule{RefTriangle}(2))
# end
@testset "Unordered nodes 3D" begin
dim = 2
nodes = [
Node((-1.0, 0.0, 0.0)), Node((0.0, 0.0, 0.0)), Node((1.0, 0.0, 0.0)),
Node((-1.0, 1.0, 0.0)), Node((0.0, 1.0, 0.0)), Node((1.0, 1.0, 0.0)),
Node((-1.0, 0.0, 1.0)), Node((0.0, 0.0, 1.0)), Node((1.0, 0.0, 1.0)),
Node((-1.0, 1.0, 1.0)), Node((0.0, 1.0, 1.0)), Node((1.0, 1.0, 1.0)),
]
cells = [
Hexahedron((1, 2, 5, 4, 7, 8, 11, 10)),
Hexahedron((5, 6, 12, 11, 2, 3, 9, 8)),
]

grid = Grid(cells, nodes)
test_interfacevalues(
grid,
InterfaceValues(FacetQuadratureRule{RefHexahedron}(2), DiscontinuousLagrange{RefHexahedron, 1}())
)
@testset "Hexahedron" begin
nodes = [
Node((-1.0, 0.0, 0.0)), Node((0.0, 0.0, 0.0)), Node((1.0, 0.0, 0.0)),
Node((-1.0, 1.0, 0.0)), Node((0.0, 1.0, 0.0)), Node((1.0, 1.0, 0.0)),
Node((-1.0, 0.0, 1.0)), Node((0.0, 0.0, 1.0)), Node((1.0, 0.0, 1.0)),
Node((-1.0, 1.0, 1.0)), Node((0.0, 1.0, 1.0)), Node((1.0, 1.0, 1.0)),
]
cells = [
Hexahedron((1, 2, 5, 4, 7, 8, 11, 10)),
Hexahedron((5, 6, 12, 11, 2, 3, 9, 8)),
]
grid = Grid(cells, nodes)
test_interfacevalues(
grid,
InterfaceValues(FacetQuadratureRule{RefHexahedron}(2), DiscontinuousLagrange{RefHexahedron, 1}())
)
orientation_info = Ferrite.InterfaceOrientationInfo(getcells(grid, 1), getcells(grid, 2), 3, 5)
@testset "Interface Orientation" begin
@test orientation_info.flipped == true
@test Ferrite.get_transformation_matrix(orientation_info) isa Tensor{2, 3}
end
@testset "Flipped normal Interface Orientation" begin
nodes = [
Node((-1.0, 0.0, 0.0)), Node((0.0, 0.0, 0.0)), Node((1.0, 0.0, 0.0)),
Node((-1.0, 1.0, 0.0)), Node((0.0, 1.0, 0.0)), Node((1.0, 1.0, 0.0)),
Node((-1.0, 0.0, 1.0)), Node((0.0, 0.0, 1.0)), Node((1.0, 0.0, 1.0)),
Node((-1.0, 1.0, 1.0)), Node((0.0, 1.0, 1.0)), Node((1.0, 1.0, 1.0)),
]
cells = [
Hexahedron((1, 4, 5, 2, 7, 10, 11, 8)),
Hexahedron((5, 6, 12, 11, 2, 3, 9, 8)),
]
grid = Grid(cells, nodes)
orientation_info = Ferrite.InterfaceOrientationInfo(getcells(grid, 1), getcells(grid, 2), 4, 5)
@test orientation_info.flipped == false
@test Ferrite.get_transformation_matrix(orientation_info) isa Tensor{2, 3}
end
end
@testset "Tetrahedron" begin
nodes = [
Node((0.0, 0.0, 0.0)), Node((1.0, 0.0, 0.0)), Node((0.0, 1.0, 0.0)),
Node((0.0, 0.0, 1.0)), Node((-1.0, 0.0, 0.0)),
]
cells = [
Tetrahedron((1, 2, 3, 4)),
Tetrahedron((1, 3, 5, 4)),
]
grid = Grid(cells, nodes)
test_interfacevalues(
grid,
InterfaceValues(FacetQuadratureRule{RefTetrahedron}(2), DiscontinuousLagrange{RefTetrahedron, 1}())
)
orientation_info = Ferrite.InterfaceOrientationInfo(getcells(grid, 1), getcells(grid, 2), 4, 2)
@testset "Interface Orientation" begin
@test orientation_info.flipped == true
@test Ferrite.get_transformation_matrix(orientation_info) isa Tensor{2, 3}
end
@testset "Flipped normal Interface Orientation" begin
nodes = [
Node((0.0, 0.0, 0.0)), Node((1.0, 0.0, 0.0)), Node((0.0, 1.0, 0.0)),
Node((0.0, 0.0, 1.0)), Node((-1.0, 0.0, 0.0)),
]
cells = [
Tetrahedron((1, 2, 4, 3)),
Tetrahedron((1, 3, 5, 4)),
]
grid = Grid(cells, nodes)
orientation_info = Ferrite.InterfaceOrientationInfo(getcells(grid, 1), getcells(grid, 2), 4, 2)
@test orientation_info.flipped == false
@test Ferrite.get_transformation_matrix(orientation_info) isa Tensor{2, 3}
end
end
end
@testset "Interface dof_range" begin
grid = generate_grid(Quadrilateral, (3, 3))
Expand Down

0 comments on commit e7a671d

Please sign in to comment.