diff --git a/src/Dofs/DofHandler.jl b/src/Dofs/DofHandler.jl index e5ea543fd9..fbfc7d83b4 100644 --- a/src/Dofs/DofHandler.jl +++ b/src/Dofs/DofHandler.jl @@ -440,7 +440,11 @@ function _close_subdofhandler!(dh::DofHandler{sdim}, sdh::SubDofHandler, sdh_ind # TODO: More than one face dof per face in 3D are not implemented yet. This requires # keeping track of the relative rotation between the faces, not just the # direction as for faces (i.e. edges) in 2D. - sdim == 3 && @assert !any(x -> x > 1, ip_info.nfacedofs) + # + # PR #1199: + # Either the facedofs must be permuted during dof distribution (if possible), or they must be transformed during mapping. + # For now, if multiple facedofs are defined, the interpolation must have transform_facedofs_during_mapping = true. + sdim == 3 && any(x -> x > 1, ip_info.nfacedofs) && @assert transform_facedofs_during_mapping(interpolation) end # TODO: Given the InterpolationInfo it should be possible to compute ndofs_per_cell, but @@ -644,7 +648,7 @@ described therein. # typestable way. n_copies = step(dofs) @assert n_copies > 0 - if adjust_during_distribution && !orientation.regular + if adjust_during_distribution && orientation.flipped # Reverse the dofs for the path dofs = reverse(dofs) end @@ -661,12 +665,12 @@ end Returns the unique representation of an edge and its orientation. Here the unique representation is the sorted node index tuple. The -orientation is `true` if the edge is not flipped, where it is `false` +orientation is `false` if the edge is not flipped, where it is `true` if the edge is flipped. """ function sortedge(edge::Tuple{Int, Int}) a, b = edge - return a < b ? (edge, PathOrientationInfo(true)) : ((b, a), PathOrientationInfo(false)) + return a < b ? (edge, PathOrientationInfo(false)) : ((b, a), PathOrientationInfo(true)) end """ @@ -721,9 +725,10 @@ For more details we refer to [1] as we follow the methodology described therein. !!!TODO Investigate if we can somehow pass the interpolation into this function in a typestable way. """ @inline function permute_and_push!(cell_dofs::Vector{Int}, dofs::StepRange{Int, Int}, ::SurfaceOrientationInfo, adjust_during_distribution::Bool, rdim::Int) - if rdim == 3 && adjust_during_distribution && length(dofs) > 1 - error("Dof distribution for interpolations with multiple dofs per face not implemented yet.") - end + # TODO test if facedof transformation/permutation is implemented + #if rdim == 3 && adjust_during_distribution && length(dofs) > 1 + # error("Dof distribution for interpolations with multiple dofs per face not implemented yet.") + #end n_copies = step(dofs) @assert n_copies > 0 for dof in dofs @@ -739,7 +744,7 @@ function sortface(face::Tuple{Int, Int, Int}) b, c = minmax(b, c) a, c = minmax(a, c) a, b = minmax(a, b) - return (a, b, c), SurfaceOrientationInfo() # TODO fill struct + return (a, b, c), SurfaceOrientationInfo(face) end @@ -760,7 +765,7 @@ function sortface(face::Tuple{Int, Int, Int, Int}) b, c = minmax(b, c) a, c = minmax(a, c) a, b = minmax(a, b) - return (a, b, c), SurfaceOrientationInfo() # TODO fill struct + return (a, b, c), SurfaceOrientationInfo(face) end diff --git a/src/FEValues/FunctionValues.jl b/src/FEValues/FunctionValues.jl index e5ea8c3fe4..64bfb959b7 100644 --- a/src/FEValues/FunctionValues.jl +++ b/src/FEValues/FunctionValues.jl @@ -221,10 +221,15 @@ end # Covariant Piola Mapping @inline function apply_mapping!(funvals::FunctionValues{0}, ::CovariantPiolaMapping, q_point::Int, mapping_values, cell) Jinv = inv(getjacobian(mapping_values)) + + Nξ = @view funvals.Nξ[:, q_point] + ip = funvals.ip + + Mξ, _ = transform_dofs(ip, cell, Nξ, Nξ) + @inbounds for j in 1:getnbasefunctions(funvals) - d = get_direction(funvals.ip, j, cell) - Nξ = funvals.Nξ[j, q_point] - funvals.Nx[j, q_point] = d * (Nξ ⋅ Jinv) + mξ = Mξ[j] + funvals.Nx[j, q_point] = mξ ⋅ Jinv end return nothing end @@ -232,13 +237,20 @@ end @inline function apply_mapping!(funvals::FunctionValues{1}, ::CovariantPiolaMapping, q_point::Int, mapping_values, cell) H = gethessian(mapping_values) Jinv = inv(getjacobian(mapping_values)) + + Nξ = @view funvals.Nξ[:, q_point] + dNdξ = @view funvals.dNdξ[:, q_point] + ip = funvals.ip + + Mξ, dMdξ = transform_dofs(ip, cell, Nξ, dNdξ) + @inbounds for j in 1:getnbasefunctions(funvals) - d = get_direction(funvals.ip, j, cell) - dNdξ = funvals.dNdξ[j, q_point] - Nξ = funvals.Nξ[j, q_point] - funvals.Nx[j, q_point] = d * (Nξ ⋅ Jinv) - funvals.dNdx[j, q_point] = d * (Jinv' ⋅ dNdξ ⋅ Jinv - Jinv' ⋅ (Nξ ⋅ Jinv ⋅ H ⋅ Jinv)) + mξ = Mξ[j] + dmdξ = dMdξ[j] + funvals.Nx[j, q_point] = mξ ⋅ Jinv + funvals.dNdx[j, q_point] = Jinv' ⋅ dmdξ ⋅ Jinv - Jinv' ⋅ (mξ ⋅ Jinv ⋅ H ⋅ Jinv) end + return nothing end @@ -246,10 +258,15 @@ end @inline function apply_mapping!(funvals::FunctionValues{0}, ::ContravariantPiolaMapping, q_point::Int, mapping_values, cell) J = getjacobian(mapping_values) detJ = det(J) + + Nξ = @view funvals.Nξ[:, q_point] + ip = funvals.ip + + Mξ, _ = transform_dofs(ip, cell, Nξ, Nξ) + @inbounds for j in 1:getnbasefunctions(funvals) - d = get_direction(funvals.ip, j, cell) - Nξ = funvals.Nξ[j, q_point] - funvals.Nx[j, q_point] = d * (J ⋅ Nξ) / detJ + mξ = Mξ[j] + funvals.Nx[j, q_point] = J ⋅ mξ / detJ end return nothing end @@ -263,12 +280,81 @@ end H_Jinv = H ⋅ Jinv A1 = (H_Jinv ⊡ (otimesl(I2, I2))) / detJ A2 = (Jinv' ⊡ H_Jinv) / detJ + + Nξ = @view funvals.Nξ[:, q_point] + dNdξ = @view funvals.dNdξ[:, q_point] + ip = funvals.ip + + Mξ, dMdξ = transform_dofs(ip, cell, Nξ, dNdξ) + @inbounds for j in 1:getnbasefunctions(funvals) - d = get_direction(funvals.ip, j, cell) - dNdξ = funvals.dNdξ[j, q_point] - Nξ = funvals.Nξ[j, q_point] - funvals.Nx[j, q_point] = d * (J ⋅ Nξ) / detJ - funvals.dNdx[j, q_point] = d * (J ⋅ dNdξ ⋅ Jinv / detJ + A1 ⋅ Nξ - (J ⋅ Nξ) ⊗ A2) + mξ = Mξ[j] + dmdξ = dMdξ[j] + funvals.Nx[j, q_point] = J ⋅ mξ / detJ + funvals.dNdx[j, q_point] = J ⋅ dmdξ ⋅ Jinv / detJ + A1 ⋅ mξ - (J ⋅ mξ) ⊗ A2 end return nothing end + +# ============== +# Transform dofs +# ============== +function transform_dofs(ip::Interpolation, cell::AbstractCell, args...) + return transform_dofs(ip, conformity(ip), cell, args...) +end + +function transform_dofs(ip::Interpolation, ::HcurlConformity, cell::AbstractCell, Nξ, dNdξ) + Mξ = similar(Nξ) + dMdξ = similar(dNdξ) + + @inbounds for (edge, edgedofs) in enumerate(edgedof_interior_indices(ip)) + @inbounds for dof in edgedofs + d = get_direction(ip, dof, cell) + Mξ[dof] = d * Nξ[dof] + dMdξ[dof] = d * dNdξ[dof] + end + end + + @inbounds for (face, facedofs) in enumerate(facedof_interior_indices(ip)) + isempty(facedofs) && continue + + orientation = SurfaceOrientationInfo(faces(cell)[face]) + + dof_idx = collect(facedofs) + nξ = Nξ[dof_idx] + dndξ = dNdξ[dof_idx] + + # Face DOF transformation on elements that need it + # TODO it appears that the required mapping is different for each interpolation. Possible to generalize somehow? + if (transform_facedofs_during_mapping(ip)) + T = dof_transformation_matrix(ip, orientation, face) + nξ = T * nξ + dndξ = T * dndξ + end + + @inbounds for (j, dof) in enumerate(facedofs) + Mξ[dof] = nξ[j] + dMdξ[dof] = dndξ[j] + end + end + + return Mξ, dMdξ +end + +function transform_dofs(ip::Interpolation, ::HdivConformity, cell::AbstractCell, Nξ, dNdξ) + Mξ = similar(Nξ) + dMdξ = similar(dNdξ) + + @inbounds for j in 1:getnbasefunctions(ip) + d = get_direction(ip, j, cell) + Mξ[j] = d * Nξ[j] + dMdξ[j] = d * dNdξ[j] + end + + @inbounds for (face, facedofs) in enumerate(facedof_interior_indices(ip)) + isempty(facedofs) && continue + getrefdim(ip) == 3 && length(facedofs) > 1 && error("Multiple facedofs not supported for 3D H(div) interpolations: $ip") + end + + return Mξ, dMdξ +end diff --git a/src/FEValues/InterfaceValues.jl b/src/FEValues/InterfaceValues.jl index 64ea4dc3fd..e06c993832 100644 --- a/src/FEValues/InterfaceValues.jl +++ b/src/FEValues/InterfaceValues.jl @@ -419,8 +419,8 @@ Return the relative orientation info for facet B with regards to facet A. Relative orientation is computed using a [`OrientationInfo`](@ref) for each side of the interface. """ function InterfaceOrientationInfo(cell_a::AbstractCell{RefShapeA}, cell_b::AbstractCell{RefShapeB}, facet_a::Int, facet_b::Int) where {RefShapeA <: AbstractRefShape, RefShapeB <: AbstractRefShape} - OI_a = OrientationInfo(facets(cell_a)[facet_a]) - OI_b = OrientationInfo(facets(cell_b)[facet_b]) + OI_a = orientation_info(facets(cell_a)[facet_a]) + OI_b = orientation_info(facets(cell_b)[facet_b]) flipped = OI_a.flipped != OI_b.flipped shift_index = OI_b.shift_index - OI_a.shift_index return InterfaceOrientationInfo{RefShapeA, RefShapeB}(flipped, shift_index, OI_b.shift_index, facet_a, facet_b) diff --git a/src/Ferrite.jl b/src/Ferrite.jl index c1abb64f9d..1242826a52 100644 --- a/src/Ferrite.jl +++ b/src/Ferrite.jl @@ -116,6 +116,13 @@ include("PoolAllocator.jl") # Matrix/Vector utilities include("arrayutils.jl") +# Grid +include("Grid/grid.jl") +include("Grid/topology.jl") +include("Grid/utils.jl") +include("Grid/grid_generators.jl") +include("Grid/coloring.jl") + # Interpolations include("interpolations.jl") @@ -133,13 +140,6 @@ include("FEValues/PointValues.jl") include("FEValues/common_values.jl") include("FEValues/facet_integrals.jl") -# Grid -include("Grid/grid.jl") -include("Grid/topology.jl") -include("Grid/utils.jl") -include("Grid/grid_generators.jl") -include("Grid/coloring.jl") - # Dofs include("Dofs/DofHandler.jl") include("Dofs/ConstraintHandler.jl") diff --git a/src/Grid/grid.jl b/src/Grid/grid.jl index 2d3ca08552..2daa8a823c 100644 --- a/src/Grid/grid.jl +++ b/src/Grid/grid.jl @@ -664,8 +664,28 @@ end ################################# # @TODO merge this code with into the logic in `ConstraintHandler`. +function get_edge_direction(edgenodes::NTuple{2, Int}) + positive = edgenodes[2] > edgenodes[1] + return ifelse(positive, 1, -1) +end + +function get_face_direction(facenodes::NTuple{N, Int}) where {N} + N > 2 || throw(ArgumentError("A face must have at least 3 nodes")) + min_idx = argmin(facenodes) + if min_idx == 1 + positive = facenodes[2] < facenodes[end] + elseif min_idx == length(facenodes) + positive = facenodes[1] < facenodes[end - 1] + else + positive = facenodes[min_idx + 1] < facenodes[min_idx - 1] + end + return ifelse(positive, 1, -1) +end + """ PathOrientationInfo + PathOrientationInfo(flipped:Bool) + PathOrientationInfo(edgenodes::NTuple{2, Int}) Orientation information for 1D entities. @@ -674,117 +694,73 @@ associated to the vertices. To give an example, the oriented path ``` 1 ---> 2 ``` -is called *regular*, indicated by `regular=true`, while the oriented path +is called *regular*, indicated by `flipped=false`, while the oriented path ``` 2 ---> 1 ``` -is called *inverted*, indicated by `regular=false`. +is called *inverted*, indicated by `flipped=true`. """ struct PathOrientationInfo - regular::Bool # Indicator whether the orientation is regular or inverted. + flipped::Bool # Indicator whether the orientation is regular or inverted. + shift_index::Int + + function PathOrientationInfo(flipped::Bool) + return new(flipped, 0) + end +end + +function PathOrientationInfo(edgenodes::NTuple{2, Int}) + return PathOrientationInfo(get_edge_direction(edgenodes) < 0) end """ SurfaceOrientationInfo + SurfaceOrientationInfo(flipped:Bool, shift_index::Int) + SurfaceOrientationInfo(facenodes::NTuple{N, Int}) Orientation information for 2D entities. Such an entity can be possibly flipped (i.e. the defining vertex order is reverse to the spanning vertex order) and the vertices can be rotated against each other. Take for example the faces ``` -1---2 2---3 +4---3 3---2 | A | | B | -4---3 1---4 +1---2 4---1 ``` which are rotated against each other by 90° (shift index is 1) or the faces ``` -1---2 2---1 +4---3 2---3 | A | | B | -4---3 3---4 +1---2 1---4 ``` which are flipped against each other. Any combination of these can happen. The combination to map this local face to the defining face is encoded with this data structure via ``rotate \\circ flip`` where the rotation is indiced by the shift index. - !!!NOTE TODO implement me. """ struct SurfaceOrientationInfo - #flipped::Bool - #shift_index::Int -end - - -@doc raw""" - InterfaceOrientationInfo - -Orientation information for 1D and 2D entities. -The orientation is defined by the indices of the grid nodes -associated to the vertices. To give an example, the oriented path -``` -1 ---> 2 -``` -is called *regular*, indicated by `flipped=false`, while the oriented path -``` -2 ---> 1 -``` -is called *inverted*, indicated by `flipped=true`. - -2D entities can be flipped (i.e. the defining vertex order is reverse to the -spanning vertex order) and the vertices can be rotated against each other. - -The reference entity is a one with it's first node is the lowest index vertex -and its vertices span counter-clock-wise. -Take for example the faces -``` -1 2 -| \ | \ -| \ | \ -| A \ | B \ -| \ | \ -2-----3 3-----1 -``` -which are rotated against each other by 240° after tranfroming to an -equilateral triangle (shift index is 2). Or the faces -``` -3 2 -| \ | \ -| \ | \ -| A \ | B \ -| \ | \ -2-----1 3-----1 -``` -which are flipped against each other. -""" -struct OrientationInfo flipped::Bool shift_index::Int end -function OrientationInfo(edgenodes::NTuple{2, Int}) - return OrientationInfo(get_edge_direction(edgenodes) < 0, 0) -end - -function OrientationInfo(facenodes::NTuple{N, Int}) where {N} +function SurfaceOrientationInfo(facenodes::NTuple{N, Int}) where {N} min_idx = argmin(facenodes) shift_index = min_idx - 1 flipped = get_face_direction(facenodes) < 0 - return OrientationInfo(flipped, shift_index) + return SurfaceOrientationInfo(flipped, shift_index) end -function get_edge_direction(edgenodes::NTuple{2, Int}) - positive = edgenodes[2] > edgenodes[1] - return ifelse(positive, 1, -1) + +""" + orientation_info(edgenodes::NTuple{2, Int}) + orientation_info(facenodes::NTuple{N, Int}) + +Helper function that returns a `PathOrientationInfo` or `SurfaceOrientationInfo` as appropriate for the size of the node tuple. +""" +function orientation_info(edgenodes::NTuple{2, Int}) + return PathOrientationInfo(edgenodes) end -function get_face_direction(facenodes::NTuple{N, Int}) where {N} - N > 2 || throw(ArgumentError("A face must have at least 3 nodes")) - min_idx = argmin(facenodes) - if min_idx == 1 - positive = facenodes[2] < facenodes[end] - elseif min_idx == length(facenodes) - positive = facenodes[1] < facenodes[end - 1] - else - positive = facenodes[min_idx + 1] < facenodes[min_idx - 1] - end - return ifelse(positive, 1, -1) +function orientation_info(facenodes::NTuple{N, Int}) where {N} + return SurfaceOrientationInfo(facenodes) end diff --git a/src/interpolations.jl b/src/interpolations.jl index 1203e45a5d..756e3d7b34 100644 --- a/src/interpolations.jl +++ b/src/interpolations.jl @@ -99,6 +99,9 @@ interpolations, generally). """ adjust_dofs_during_distribution(::Interpolation) + +transform_facedofs_during_mapping(::Interpolation) = false + """ InterpolationInfo @@ -2114,6 +2117,77 @@ function get_direction(::Nedelec{RefTetrahedron, 1}, shape_nr, cell) return get_edge_direction(cell, shape_nr) # shape_nr = edge_nr end +# RefTetrahedron, 2nd order Lagrange +# https://defelement.org/elements/examples/tetrahedron-nedelec1-lagrange-1.html +function reference_shape_value(ip::Nedelec{RefTetrahedron, 2}, ξ::Vec{3, T}, i::Int) where {T} + x, y, z = ξ + nil = zero(T) + + # Edge 1 + i == 1 && return Vec(8(y^2 + z^2 + x * y + x * z + 2y * z) - 6(x + 2y + 2z) + 4, -2x * (4x + 4y + 4z - 3), -2x * (4x + 4y + 4z - 3)) + i == 2 && return Vec(-8x * (y + z) + 2(3x + y + z - 1), 4x * (2x - 1), 4x * (2x - 1)) + # Edge 2 + i == 3 && return Vec(2y * (1 - 4x), 4x * (2x - 1), nil) + i == 4 && return Vec(4y * (1 - 2y), 2x * (4y - 1), nil) + # Edge 3 (flip order and sign compared to defelement) + i == 5 && return Vec(-4y * (2y - 1), 8y * (x + z) - 2(x + 3y + z - 1), -4y * (2y - 1)) + i == 6 && return Vec(2y * (4x + 4y + 4z - 3), -(8(x^2 + z^2 + x * y + 2x * z + y * z) - 6(2x + y + 2z) + 4), 2y * (4x + 4y + 4z - 3)) + # Edge 4 + i == 7 && return Vec(-2z * (4x + 4y + 4z - 3), -2z * (4x + 4y + 4z - 3), 8(x^2 + y^2 + 2x * y + x * z + y * z) - 6(2x + 2y + z) + 4) + i == 8 && return Vec(4z * (2z - 1), 4z * (2z - 1), -8z * (x + y) + 2(x + y + 3z - 1)) + # Edge 5 + i == 9 && return Vec(2z * (1 - 4x), nil, 4x * (2x - 1)) + i == 10 && return Vec(4z * (1 - 2z), nil, 2x * (4z - 1)) + # Edge 6 + i == 11 && return Vec(nil, 2z * (1 - 4y), 4y * (2y - 1)) + i == 12 && return Vec(nil, 4z * (1 - 2z), 2y * (4z - 1)) + + # Face 1 (flip order compared to defelement) + i == 13 && return Vec(8y * (2x + y + z - 1), -8x * (2x + y + 2z - 2), 8x * y) + i == 14 && return Vec(-8y * (x + 2y + 2z - 2), 8x * (x + 2y + z - 1), 8x * y) + # Face 2 + i == 15 && return Vec(-8z * (x + 2y + 2z - 2), 8x * z, 8x * (x + y + 2z - 1)) + i == 16 && return Vec(8z * (2x + y + z - 1), 8x * z, -8x * (2x + 2y + z - 2)) + # Face 3 + i == 17 && return Vec(-8y * z, 16x * z, -8x * y) + i == 18 && return Vec(-8y * z, -8x * z, 16x * y) + # Face 4 (flip order compared to defelement) + i == 19 && return Vec(8y * z, 8z * (x + 2y + z - 1), -8y * (2x + 2y + z - 2)) + i == 20 && return Vec(8y * z, -8z * (2x + y + 2z - 2), 8y * (x + y + 2z - 1)) + throw(ArgumentError("no shape function $i for interpolation $ip")) +end + +getnbasefunctions(::Nedelec{RefTetrahedron, 2}) = 20 +edgedof_interior_indices(::Nedelec{RefTetrahedron, 2}) = ((1, 2), (3, 4), (5, 6), (7, 8), (9, 10), (11, 12)) +facedof_indices(::Nedelec{RefTetrahedron, 2}) = ((1, 2, 3, 4, 5, 6, 13, 14), (1, 2, 7, 8, 9, 10, 15, 16), (3, 4, 9, 10, 11, 12, 17, 18), (5, 6, 7, 8, 11, 12, 19, 20)) +facedof_interior_indices(::Nedelec{RefTetrahedron, 2}) = ((13, 14), (15, 16), (17, 18), (19, 20)) +adjust_dofs_during_distribution(::Nedelec{RefTetrahedron, 2}) = true +transform_facedofs_during_mapping(::Nedelec{RefTetrahedron, 2}) = true + +function get_direction(::Nedelec{RefTetrahedron, 2}, shape_nr, cell) + if (shape_nr < 13) + edge_nr = (shape_nr + 1) ÷ 2 + return get_edge_direction(cell, edge_nr) + else + face_nr = (shape_nr - 11) ÷ 2 + return get_face_direction(cell, face_nr) + end +end + +function dof_transformation_matrix(ip::Nedelec{RefTetrahedron, 2}, orientation::SurfaceOrientationInfo, facenr::Int) + shift_index = orientation.shift_index + if (!orientation.flipped) + shift_index == 0 && return Tensor{2, 2}((1, 0, 0, 1)) + shift_index == 1 && return Tensor{2, 2}((0, -1, 1, -1)) + shift_index == 2 && return Tensor{2, 2}((-1, 1, -1, 0)) + else + shift_index == 0 && return Tensor{2, 2}((0, 1, 1, 0)) + shift_index == 1 && return Tensor{2, 2}((-1, 0, -1, 1)) + shift_index == 2 && return Tensor{2, 2}((1, -1, 0, -1)) + end + throw(ArgumentError("no dof transformation matrix defined for $orientation for interpolation $ip")) +end + # RefHexahedron, 1st order Lagrange # https://defelement.org/elements/examples/hexahedron-nedelec1-lagrange-1.html function reference_shape_value(ip::Nedelec{RefHexahedron, 1}, ξ::Vec{3, T}, i::Int) where {T} diff --git a/test/runtests.jl b/test/runtests.jl index 6ba1a352dc..71909bcb69 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -18,6 +18,7 @@ include("test_utils.jl") # Unit tests include("test_collectionsofviews.jl") include("test_refshapes.jl") +include("test_orientations.jl") include("test_interpolations.jl") include("test_cellvalues.jl") include("test_facevalues.jl") diff --git a/test/test_cellvalues.jl b/test/test_cellvalues.jl index daffd081b4..b3e01b73df 100644 --- a/test/test_cellvalues.jl +++ b/test/test_cellvalues.jl @@ -236,12 +236,14 @@ end test_ips = [ Lagrange{RefTriangle, 2}(), Lagrange{RefQuadrilateral, 2}(), Lagrange{RefHexahedron, 2}()^3, # Test should also work for identity mapping - Nedelec{RefTriangle, 1}(), Nedelec{RefTriangle, 2}(), - RaviartThomas{RefTriangle, 1}(), RaviartThomas{RefTriangle, 2}(), BrezziDouglasMarini{RefTriangle, 1}(), + Nedelec{RefTriangle, 1}(), Nedelec{RefTriangle, 2}(), Nedelec{RefQuadrilateral, 1}(), Nedelec{RefTetrahedron, 1}(), Nedelec{RefTetrahedron, 2}(), Nedelec{RefHexahedron, 1}(), + RaviartThomas{RefTriangle, 1}(), RaviartThomas{RefTriangle, 2}(), RaviartThomas{RefQuadrilateral, 1}(), RaviartThomas{RefTetrahedron, 1}(), RaviartThomas{RefHexahedron, 1}(), + BrezziDouglasMarini{RefTriangle, 1}(), ] # Makes most sense to test for nonlinear geometries, so we use such cells for testing only cell_from_refshape(::Type{RefTriangle}) = QuadraticTriangle((ntuple(identity, 6))) cell_from_refshape(::Type{RefQuadrilateral}) = QuadraticQuadrilateral((ntuple(identity, 9))) + cell_from_refshape(::Type{RefTetrahedron}) = QuadraticTetrahedron((ntuple(identity, 10))) cell_from_refshape(::Type{RefHexahedron}) = QuadraticHexahedron((ntuple(identity, 27))) for ip in test_ips cell = cell_from_refshape(getrefshape(ip)) diff --git a/test/test_continuity.jl b/test/test_continuity.jl index 8723bd1e85..ed0625f8f9 100644 --- a/test/test_continuity.jl +++ b/test/test_continuity.jl @@ -129,7 +129,7 @@ test_ips = [ Lagrange{RefTriangle, 2}(), Lagrange{RefQuadrilateral, 2}(), Lagrange{RefHexahedron, 2}()^3, # Test should also work for identity mapping - Nedelec{RefTriangle, 1}(), Nedelec{RefTriangle, 2}(), Nedelec{RefQuadrilateral, 1}(), Nedelec{RefTetrahedron, 1}(), Nedelec{RefHexahedron, 1}(), + Nedelec{RefTriangle, 1}(), Nedelec{RefTriangle, 2}(), Nedelec{RefQuadrilateral, 1}(), Nedelec{RefTetrahedron, 1}(), Nedelec{RefTetrahedron, 2}(), Nedelec{RefHexahedron, 1}(), RaviartThomas{RefTriangle, 1}(), RaviartThomas{RefTriangle, 2}(), RaviartThomas{RefQuadrilateral, 1}(), RaviartThomas{RefTetrahedron, 1}(), RaviartThomas{RefHexahedron, 1}(), BrezziDouglasMarini{RefTriangle, 1}(), ] diff --git a/test/test_interpolations.jl b/test/test_interpolations.jl index 8145fb8067..624d6ea509 100644 --- a/test/test_interpolations.jl +++ b/test/test_interpolations.jl @@ -360,7 +360,7 @@ end @testset "H(curl) and H(div)" begin Hcurl_interpolations = [ Nedelec{RefTriangle, 1}(), Nedelec{RefTriangle, 2}(), Nedelec{RefQuadrilateral, 1}(), - Nedelec{RefTetrahedron, 1}(), Nedelec{RefHexahedron, 1}(), + Nedelec{RefTetrahedron, 1}(), Nedelec{RefTetrahedron, 2}(), Nedelec{RefHexahedron, 1}(), ] Hdiv_interpolations = [ RaviartThomas{RefTriangle, 1}(), RaviartThomas{RefTriangle, 2}(), RaviartThomas{RefQuadrilateral, 1}(), @@ -372,16 +372,16 @@ end ## Raviart-Thomas on RefTriangle reference_edge_moment(::RaviartThomas{RefTriangle, 1}, edge_shape_nr, s) = 1 reference_edge_moment(::RaviartThomas{RefTriangle, 2}, edge_shape_nr, s) = edge_shape_nr == 1 ? (1 - s) : s - reference_face_moment(::RaviartThomas{RefTriangle, 2}, face_shape_nr, s1, s2) = face_shape_nr == 1 ? Vec(1, 0) : Vec(0, 1) + reference_face_moment(::RaviartThomas{RefTriangle, 2}, face_shape_nr, shape_nr, s1, s2) = face_shape_nr == 1 ? Vec(1, 0) : Vec(0, 1) ## Raviart-Thomas on RefQuadrilateral reference_edge_moment(::RaviartThomas{RefQuadrilateral, 1}, edge_shape_nr, s) = 1 ## Raviart-Thomas on RefTetrahedron - reference_face_moment(::RaviartThomas{RefTetrahedron, 1}, face_shape_nr, s1, s2) = 1 + reference_face_moment(::RaviartThomas{RefTetrahedron, 1}, face_shape_nr, shape_nr, s1, s2) = 1 ## Raviart-Thomas on RefHexahedron - reference_face_moment(::RaviartThomas{RefHexahedron, 1}, face_shape_nr, s1, s2) = 1 + reference_face_moment(::RaviartThomas{RefHexahedron, 1}, face_shape_nr, shape_nr, s1, s2) = 1 ## Brezzi-Douglas-Marini on RefTriangle reference_edge_moment(::BrezziDouglasMarini{RefTriangle, 1}, edge_shape_nr, s) = edge_shape_nr == 1 ? (1 - s) : s @@ -389,13 +389,25 @@ end ## Nedelec on RefTriangle reference_edge_moment(::Nedelec{RefTriangle, 1}, edge_shape_nr, s) = 1 reference_edge_moment(::Nedelec{RefTriangle, 2}, edge_shape_nr, s) = edge_shape_nr == 1 ? (1 - s) : s - reference_face_moment(::Nedelec{RefTriangle, 2}, face_shape_nr, s1, s2) = face_shape_nr == 1 ? Vec(1, 0) : Vec(0, 1) + reference_face_moment(::Nedelec{RefTriangle, 2}, face_shape_nr, shape_nr, s1, s2) = face_shape_nr == 1 ? Vec(1, 0) : Vec(0, 1) ## Nedelec on RefQuadrilateral reference_edge_moment(::Nedelec{RefQuadrilateral, 1}, edge_shape_nr, s) = 1 ## Nedelec on RefTetrahedron reference_edge_moment(::Nedelec{RefTetrahedron, 1}, edge_shape_nr, s) = 1 + reference_edge_moment(::Nedelec{RefTetrahedron, 2}, edge_shape_nr, s) = edge_shape_nr == 1 ? (1 - s) : s + function reference_face_moment(::Nedelec{RefTetrahedron, 2}, face_shape_nr, shape_nr, s1, s2) + shape_nr == 13 && return Vec(0, 1, 0) + shape_nr == 14 && return Vec(1, 0, 0) + shape_nr == 15 && return Vec(1, 0, 0) + shape_nr == 16 && return Vec(0, 0, 1) + shape_nr == 17 && return Vec(-1 / √3, 1 / √3, 0.0) + shape_nr == 18 && return Vec(-1 / √3, 0.0, 1 / √3) + shape_nr == 19 && return Vec(0, 0, 1) + shape_nr == 20 && return Vec(0, 1, 0) + throw(ArgumentError("no reference edge face moment for shape_nr $shape_nr in Nedelec{RefTetrahedron, 2}")) + end ## Nedelec on RefHexahedron reference_edge_moment(::Nedelec{RefHexahedron, 1}, edge_shape_nr, s) = 1 @@ -486,13 +498,13 @@ end ξ(s1, s2) = parameterize_face(face_coords, s1, s2) for (face_shape_nr, shape_nr) in pairs(dof_inds) - f(s1, s2) = reference_face_moment(ip, face_shape_nr, s1, s2) ⋅ reference_shape_value(ip, ξ(s1, s2), shape_nr) * face_weight(ξ, s1, s2) + f(s1, s2) = reference_face_moment(ip, face_shape_nr, shape_nr, s1, s2) ⋅ reference_shape_value(ip, ξ(s1, s2), shape_nr) * face_weight(ξ, s1, s2) @test integrate_face(f) ≈ 1 # Test that the functional is zero for the other shape functions for other_shape_nr in 1:getnbasefunctions(ip) other_shape_nr == shape_nr && continue - g(s1, s2) = reference_face_moment(ip, face_shape_nr, s1, s2) ⋅ reference_shape_value(ip, ξ(s1, s2), other_shape_nr) * face_weight(ξ, s1, s2) + g(s1, s2) = reference_face_moment(ip, face_shape_nr, shape_nr, s1, s2) ⋅ reference_shape_value(ip, ξ(s1, s2), other_shape_nr) * face_weight(ξ, s1, s2) @test integrate_face(g) + 1 ≈ 1 # integrate_edge(g) ≈ 0 end end @@ -512,13 +524,13 @@ end ξ(s1, s2) = parameterize_face(face_coords, s1, s2) for (face_shape_nr, shape_nr) in pairs(dof_inds) - f(s1, s2) = reference_face_moment(ip, face_shape_nr, s1, s2) * (reference_shape_value(ip, ξ(s1, s2), shape_nr) ⋅ normal) * face_weight(ξ, s1, s2) + f(s1, s2) = reference_face_moment(ip, face_shape_nr, shape_nr, s1, s2) * (reference_shape_value(ip, ξ(s1, s2), shape_nr) ⋅ normal) * face_weight(ξ, s1, s2) @test integrate_face(f) ≈ 1 # Test that the functional is zero for the other shape functions for other_shape_nr in 1:getnbasefunctions(ip) other_shape_nr == shape_nr && continue - g(s1, s2) = reference_face_moment(ip, face_shape_nr, s1, s2) * (reference_shape_value(ip, ξ(s1, s2), other_shape_nr) ⋅ normal) * face_weight(ξ, s1, s2) + g(s1, s2) = reference_face_moment(ip, face_shape_nr, shape_nr, s1, s2) * (reference_shape_value(ip, ξ(s1, s2), other_shape_nr) ⋅ normal) * face_weight(ξ, s1, s2) @test integrate_face(g) + 1 ≈ 1 # integrate_edge(g) ≈ 0 end end @@ -569,7 +581,7 @@ end ξ(s1, s2) = parameterize_face(face_coords, s1, s2) for (face_shape_nr, shape_nr) in pairs(dof_inds) - f(s1, s2) = reference_face_moment(ip, face_shape_nr, s1, s2) ⋅ reference_shape_value(ip, ξ(s1, s2), shape_nr) * face_weight(ξ, s1, s2) + f(s1, s2) = reference_face_moment(ip, face_shape_nr, shape_nr, s1, s2) ⋅ reference_shape_value(ip, ξ(s1, s2), shape_nr) * face_weight(ξ, s1, s2) @test integrate_face(f) ≈ 1 end end diff --git a/test/test_orientations.jl b/test/test_orientations.jl new file mode 100644 index 0000000000..bede46b4c5 --- /dev/null +++ b/test/test_orientations.jl @@ -0,0 +1,75 @@ +@testset "Orientations" begin + @testset "PathOrientationInfo" begin + # Path [a, b], flipped + test_paths = [ + ((1, 2), false), # 1 -> 2, regular + ((2, 1), true), # 2 -> 1, inverted + ((5, 1), true), # 5 -> 1, inverted + ((3, 5), false), # 3 -> 5, regular + ] + + for (path, flipped) in test_paths + orientation = Ferrite.PathOrientationInfo(path) + + @test orientation.flipped == flipped + end + end + + @testset "SurfaceOrientationInfo" begin + # Surface [facenodes], flipped, shift_index + test_surfaces = [ + # 2 + # | \ + # 3 - 1, regular, rotation = 0 deg + ((1, 2, 3), false, 0), + # 1 + # | \ + # 2 - 3, regular, rotation = 120 deg + ((3, 1, 2), false, 1), + # 3 + # | \ + # 1 - 2, regular, rotation = 240 deg + ((2, 3, 1), false, 2), + # 3 + # | \ + # 2 - 1, inverted, rotation = 0 deg + ((1, 3, 2), true, 0), + # 1 + # | \ + # 3 - 2, inverted, rotation = 120 deg + ((2, 1, 3), true, 1), + # 2 + # | \ + # 1 - 3, inverted, rotation = 240 deg + ((3, 2, 1), true, 2), + # 4 - 3 + # | | + # 1 - 2, regular, rotation = 0 deg + ((1, 2, 3, 4), false, 0), + # 2 - 3 + # | | + # 1 - 4, inverted, rotation = 0 deg + ((1, 4, 3, 2), true, 0), + # 3 - 2 + # | | + # 4 - 1, regular, rotation = 90 deg + ((4, 1, 2, 3), false, 1), + # 1 - 2 + # | | + # 4 - 3, inverted, rotation = 270 deg + ((4, 3, 2, 1), true, 3), + ] + + + for (surface, flipped, shift_index) in test_surfaces + orientation = Ferrite.SurfaceOrientationInfo(surface) + + @test orientation.flipped == flipped + @test orientation.shift_index == shift_index + end + end + + @assert typeof(Ferrite.orientation_info((1, 2))) <: Ferrite.PathOrientationInfo + @assert typeof(Ferrite.orientation_info((1, 2, 3))) <: Ferrite.SurfaceOrientationInfo + @assert typeof(Ferrite.orientation_info((1, 2, 3, 4))) <: Ferrite.SurfaceOrientationInfo +end