2nd order Nedelec interpolation

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
 
 

src/FEValues/FunctionValues.jl

@@ -221,35 +221,52 @@ 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
 
 @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
 
 # Contravariant Piola Mapping
 @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

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)

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")