Move over from #798

Project.toml

@@ -49,11 +49,12 @@ Metis = "2679e427-3c69-5b7f-982b-ece356f1e94b"
 NBInclude = "0db19996-df87-5ea3-a455-e3a50d440464"
 OhMyThreads = "67456a42-1dca-4109-a031-0a68de7e3ad5"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
+QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SHA = "ea8e919c-243c-51af-8825-aaa63cd721ce"
 TaskLocalValues = "ed4db957-447d-4319-bfb6-7fa9ae7ecf34"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 
 [targets]
-test = ["BlockArrays", "Downloads", "FerriteGmsh", "ForwardDiff", "Gmsh", "IterativeSolvers", "Metis", "NBInclude", "OhMyThreads", "ProgressMeter", "Random", "SHA", "TaskLocalValues", "Test", "TimerOutputs", "Logging"]
+test = ["BlockArrays", "Downloads", "FerriteGmsh", "ForwardDiff", "Gmsh", "IterativeSolvers", "Metis", "NBInclude", "OhMyThreads", "ProgressMeter", "Random", "SHA", "TaskLocalValues", "Test", "TimerOutputs", "Logging", "QuadGK"]

src/FEValues/CellValues.jl

@@ -78,6 +78,7 @@ getdetJdV(cv::CellValues, q_point::Int) = cv.detJdV[q_point]
 getdetJdV(::CellValues{<:Any, <:Any, <:Any, Nothing}, ::Int) = throw(ArgumentError("detJdV is not saved in CellValues"))
 
 # Accessors for function values
+get_fun_values(cv::CellValues) = cv.fun_values
 getnbasefunctions(cv::CellValues) = getnbasefunctions(cv.fun_values)
 function_interpolation(cv::CellValues) = function_interpolation(cv.fun_values)
 function_difforder(cv::CellValues) = function_difforder(cv.fun_values)
@@ -111,7 +112,7 @@ function reinit!(cv::CellValues, cell::Union{AbstractCell, Nothing}, x::Abstract
     n_geom_basefuncs = getngeobasefunctions(geo_mapping)
 
     check_reinit_sdim_consistency(:CellValues, shape_gradient_type(cv), eltype(x))
-    if cell === nothing && !isa(mapping_type(fun_values), IdentityMapping)
+    if cell === nothing && reinit_needs_cell(cv)
         throw(ArgumentError("The cell::AbstractCell input is required to reinit! non-identity function mappings"))
     end
     if !checkbounds(Bool, x, 1:n_geom_basefuncs) || length(x) != n_geom_basefuncs

src/FEValues/FacetValues.jl

@@ -135,7 +135,7 @@ function reinit!(fv::FacetValues, cell::Union{AbstractCell, Nothing}, x::Abstrac
     geo_mapping = get_geo_mapping(fv)
     fun_values = get_fun_values(fv)
 
-    if cell === nothing && !isa(mapping_type(fun_values), IdentityMapping)
+    if cell === nothing && reinit_needs_cell(fv)
         throw(ArgumentError("The cell::AbstractCell input is required to reinit! non-identity function mappings"))
     end
 

src/FEValues/FunctionValues.jl

@@ -143,11 +143,8 @@ end
 
 # Mapping types
 struct IdentityMapping end
-# Not yet implemented:
-# struct CovariantPiolaMapping end # PR798
-# struct ContravariantPiolaMapping end # PR798
-# struct DoubleCovariantPiolaMapping end
-# struct DoubleContravariantPiolaMapping end
+struct CovariantPiolaMapping end
+struct ContravariantPiolaMapping end
 
 mapping_type(fv::FunctionValues) = mapping_type(fv.ip)
 
@@ -159,9 +156,8 @@ the function values and gradients from the reference cell
 to the physical cell geometry.
 """
 required_geo_diff_order(::IdentityMapping, fun_diff_order::Int) = fun_diff_order
-#required_geo_diff_order(::ContravariantPiolaMapping, fun_diff_order::Int) = 1 + fun_diff_order # PR798
-#required_geo_diff_order(::CovariantPiolaMapping,     fun_diff_order::Int) = 1 + fun_diff_order # PR798
-
+required_geo_diff_order(::ContravariantPiolaMapping, fun_diff_order::Int) = 1 + fun_diff_order
+required_geo_diff_order(::CovariantPiolaMapping, fun_diff_order::Int) = 1 + fun_diff_order
 
 # Support for embedded elements
 @inline calculate_Jinv(J::Tensor{2}) = inv(J)
@@ -221,6 +217,57 @@ end
     return nothing
 end
 
-# TODO in PR798, apply_mapping! for
-# * CovariantPiolaMapping
-# * ContravariantPiolaMapping
+# Covariant Piola Mapping
+@inline function apply_mapping!(funvals::FunctionValues{0}, ::CovariantPiolaMapping, q_point::Int, mapping_values, cell)
+    Jinv = inv(getjacobian(mapping_values))
+    @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)
+    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))
+    @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))
+    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)
+    @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
+    end
+    return nothing
+end
+
+@inline function apply_mapping!(funvals::FunctionValues{1}, ::ContravariantPiolaMapping, q_point::Int, mapping_values, cell)
+    H = gethessian(mapping_values)
+    J = getjacobian(mapping_values)
+    Jinv = inv(J)
+    detJ = det(J)
+    I2 = one(J)
+    H_Jinv = H ⋅ Jinv
+    A1 = (H_Jinv ⊡ (otimesl(I2, I2))) / detJ
+    A2 = (Jinv' ⊡ H_Jinv) / detJ
+    @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)
+    end
+    return nothing
+end

src/FEValues/common_values.jl

@@ -9,6 +9,13 @@ function checkquadpoint(fe_v::AbstractValues, qp::Int)
     return nothing
 end
 
+@inline function reinit_needs_cell(fe_values::AbstractValues)
+    # TODO: Might need better logic for this, but for current implementations this
+    # is ok. If someone implements a non-identity mapping that doesn't require the cell
+    # as input, this is only a slight performance issue in some cases.
+    return !isa(mapping_type(get_fun_values(fe_values)), IdentityMapping)
+end
+
 @noinline function throw_incompatible_dof_length(length_ue, n_base_funcs)
     msg = "the number of base functions ($(n_base_funcs)) does not match the length " *
         "of the vector ($(length_ue)). Perhaps you passed the global vector, " *

src/Ferrite.jl

@@ -21,7 +21,7 @@ using WriteVTK:
     WriteVTK, VTKCellTypes
 using Tensors:
     Tensors, AbstractTensor, SecondOrderTensor, SymmetricTensor, Tensor, Vec, gradient,
-    rotation_tensor, symmetric, tovoigt!, hessian, otimesu
+    rotation_tensor, symmetric, tovoigt!, hessian, otimesu, otimesl
 using ForwardDiff:
     ForwardDiff
 

src/exports.jl

@@ -17,6 +17,9 @@ export
     Lagrange,
     DiscontinuousLagrange,
     Serendipity,
+    Nedelec,
+    RaviartThomas,
+    BrezziDouglasMarini,
     getnbasefunctions,
     getrefshape,
 

src/interpolations.jl

@@ -1774,3 +1774,169 @@ function mapping_type end
 
 mapping_type(::ScalarInterpolation) = IdentityMapping()
 mapping_type(::VectorizedInterpolation) = IdentityMapping()
+
+#####################################
+# RaviartThomas (1st kind), H(div)  #
+#####################################
+struct RaviartThomas{vdim, shape, order} <: VectorInterpolation{vdim, shape, order} end
+mapping_type(::RaviartThomas) = ContravariantPiolaMapping()
+n_dbc_components(::RaviartThomas) = 1
+reference_coordinates(ip::RaviartThomas{vdim}) where {vdim} = fill(NaN * zero(Vec{vdim}), getnbasefunctions(ip))
+dirichlet_edgedof_indices(ip::RaviartThomas{2}) = edgedof_interior_indices(ip)
+dirichlet_facedof_indices(ip::RaviartThomas{3}) = facedof_interior_indices(ip)
+
+# RefTriangle
+edgedof_indices(ip::RaviartThomas{2, RefTriangle}) = edgedof_interior_indices(ip)
+facedof_indices(ip::RaviartThomas{2, RefTriangle}) = (ntuple(i -> i, getnbasefunctions(ip)),)
+
+# RefTriangle, 1st order Lagrange
+# https://defelement.com/elements/examples/triangle-raviart-thomas-lagrange-1.html
+function reference_shape_value(ip::RaviartThomas{2, RefTriangle, 1}, ξ::Vec{2}, i::Int)
+    x, y = ξ
+    i == 1 && return ξ                  # Flip sign
+    i == 2 && return Vec(x - 1, y)      # Keep sign
+    i == 3 && return Vec(x, y - 1)      # Flip sign
+    throw(ArgumentError("no shape function $i for interpolation $ip"))
+end
+
+getnbasefunctions(::RaviartThomas{2, RefTriangle, 1}) = 3
+edgedof_interior_indices(::RaviartThomas{2, RefTriangle, 1}) = ((1,), (2,), (3,))
+facedof_interior_indices(::RaviartThomas{2, RefTriangle, 1}) = ((),)
+adjust_dofs_during_distribution(::RaviartThomas) = false
+
+function get_direction(::RaviartThomas{2, RefTriangle, 1}, j, cell)
+    edge = edges(cell)[j]
+    return ifelse(edge[2] > edge[1], 1, -1)
+end
+
+# RefTriangle, 2st order Lagrange
+# https://defelement.org/elements/examples/triangle-raviart-thomas-lagrange-2.html
+function reference_shape_value(ip::RaviartThomas{2, RefTriangle, 2}, ξ::Vec{2}, i::Int)
+    x, y = ξ
+    # Face 1 (keep ordering, flip sign)
+    i == 1 && return Vec(4x * (2x - 1), 2y * (4x - 1))
+    i == 2 && return Vec(2x * (4y - 1), 4y * (2y - 1))
+    # Face 2 (flip ordering, keep signs)
+    i == 3 && return Vec(8x * y - 2x - 6y + 2, 4y * (2y - 1))
+    i == 4 && return Vec(-8x^2 - 8x * y + 12x + 6y - 4, 2y * (-4x - 4y + 3))
+    # Face 3 (keep ordering, flip sign)
+    i == 5 && return Vec(2x * (3 - 4x - 4y), -8x * y + 6x - 8y^2 + 12y - 4)
+    i == 6 && return Vec(4x * (2x - 1), 8x * y - 6x - 2y + 2)
+    # Cell
+    i == 7 && return Vec(8x * (-2x - y + 2), 8y * (-2x - y + 1))
+    i == 8 && return Vec(8x * (-2y - x + 1), 8y * (-2y - x + 2))
+    throw(ArgumentError("no shape function $i for interpolation $ip"))
+end
+
+getnbasefunctions(::RaviartThomas{2, RefTriangle, 2}) = 8
+edgedof_interior_indices(::RaviartThomas{2, RefTriangle, 2}) = ((1, 2), (3, 4), (5, 6))
+facedof_interior_indices(::RaviartThomas{2, RefTriangle, 2}) = ((7, 8),)
+adjust_dofs_during_distribution(::RaviartThomas{2, RefTriangle, 2}) = true
+
+function get_direction(::RaviartThomas{2, RefTriangle, 2}, j, cell)
+    j > 6 && return 1
+    edge = edges(cell)[(j + 1) ÷ 2]
+    return ifelse(edge[2] > edge[1], 1, -1)
+end
+
+#####################################
+# Brezzi-Douglas–Marini, H(div)     #
+#####################################
+struct BrezziDouglasMarini{vdim, shape, order} <: VectorInterpolation{vdim, shape, order} end
+mapping_type(::BrezziDouglasMarini) = ContravariantPiolaMapping()
+reference_coordinates(ip::BrezziDouglasMarini{vdim}) where {vdim} = fill(NaN * zero(Vec{vdim}), getnbasefunctions(ip))
+dirichlet_edgedof_indices(ip::BrezziDouglasMarini{2}) = edgedof_interior_indices(ip)
+n_dbc_components(::BrezziDouglasMarini) = 1
+
+# RefTriangle
+edgedof_indices(ip::BrezziDouglasMarini{2, RefTriangle}) = edgedof_interior_indices(ip)
+facedof_indices(ip::BrezziDouglasMarini{2, RefTriangle}) = (ntuple(i -> i, getnbasefunctions(ip)),)
+
+# RefTriangle, 1st order Lagrange
+# https://defelement.org/elements/examples/triangle-brezzi-douglas-marini-lagrange-1.html
+function reference_shape_value(ip::BrezziDouglasMarini{2, RefTriangle, 1}, ξ::Vec{2}, i::Int)
+    x, y = ξ
+    # Edge 1: y=1-x, n = [1, 1]/√2 (Flip sign, pos. integration outwards)
+    i == 1 && return Vec(4x, -2y) # N ⋅ n = (2√2 x - √2 (1-x)) = 3√2 x - √2
+    i == 2 && return Vec(-2x, 4y) # N ⋅ n = (-√2x + 2√2 (1-x)) = 2√2 - 3√2x
+    # Edge 2: x=0, n = [-1, 0] (reverse order to follow Ferrite convention)
+    i == 3 && return Vec(-2x - 6y + 2, 4y) # N ⋅ n = (6y - 2)
+    i == 4 && return Vec(4x + 6y - 4, -2y) # N ⋅ n = (4 - 6y)
+    # Edge 3: y=0, n = [0, -1] (Flip sign, pos. integration outwards)
+    i == 5 && return Vec(-2x, 6x + 4y - 4) # N ⋅ n = (4 - 6x)
+    i == 6 && return Vec(4x, -6x - 2y + 2) # N ⋅ n = (6x - 2)
+    throw(ArgumentError("no shape function $i for interpolation $ip"))
+end
+
+getnbasefunctions(::BrezziDouglasMarini{2, RefTriangle, 1}) = 6
+edgedof_interior_indices(::BrezziDouglasMarini{2, RefTriangle, 1}) = ((1, 2), (3, 4), (5, 6))
+adjust_dofs_during_distribution(::BrezziDouglasMarini{2, RefTriangle, 1}) = true
+
+function get_direction(::BrezziDouglasMarini{2, RefTriangle, 1}, j, cell)
+    edge = edges(cell)[(j + 1) ÷ 2]
+    return ifelse(edge[2] > edge[1], 1, -1)
+end
+
+#####################################
+# Nedelec (1st kind), H(curl)       #
+#####################################
+struct Nedelec{vdim, shape, order} <: VectorInterpolation{vdim, shape, order} end
+mapping_type(::Nedelec) = CovariantPiolaMapping()
+reference_coordinates(ip::Nedelec{vdim}) where {vdim} = fill(NaN * zero(Vec{vdim}), getnbasefunctions(ip))
+dirichlet_edgedof_indices(ip::Nedelec) = edgedof_interior_indices(ip)
+n_dbc_components(::Nedelec) = 1
+edgedof_indices(ip::Nedelec) = edgedof_interior_indices(ip)
+
+# 2D refshape (rdim == vdim for Nedelec)
+facedof_indices(ip::Nedelec{2, <:AbstractRefShape{2}}) = (ntuple(i -> i, getnbasefunctions(ip)),)
+
+# RefTriangle, 1st order Lagrange
+# https://defelement.org/elements/examples/triangle-nedelec1-lagrange-1.html
+function reference_shape_value(ip::Nedelec{2, RefTriangle, 1}, ξ::Vec{2}, i::Int)
+    x, y = ξ
+    i == 1 && return Vec(- y, x)
+    i == 2 && return Vec(- y, x - 1) # Changed signed, follow Ferrite's sign convention
+    i == 3 && return Vec(1 - y, x)
+    throw(ArgumentError("no shape function $i for interpolation $ip"))
+end
+
+getnbasefunctions(::Nedelec{2, RefTriangle, 1}) = 3
+edgedof_interior_indices(::Nedelec{2, RefTriangle, 1}) = ((1,), (2,), (3,))
+adjust_dofs_during_distribution(::Nedelec{2, RefTriangle, 1}) = false
+
+function get_direction(::Nedelec{2, RefTriangle, 1}, j, cell)
+    edge = edges(cell)[j]
+    return ifelse(edge[2] > edge[1], 1, -1)
+end
+
+# RefTriangle, 2nd order Lagrange
+# https://defelement.org/elements/examples/triangle-nedelec1-lagrange-2.html
+function reference_shape_value(ip::Nedelec{2, RefTriangle, 2}, ξ::Vec{2}, i::Int)
+    x, y = ξ
+    # Edge 1
+    i == 1 && return Vec(2 * y * (1 - 4 * x), 4 * x * (2 * x - 1))
+    i == 2 && return Vec(4 * y * (1 - 2 * y), 2 * x * (4 * y - 1))
+    # Edge 2 (flip order and sign compared to defelement)
+    i == 3 && return Vec(4 * y * (1 - 2 * y), 8 * x * y - 2 * x - 6 * y + 2)
+    i == 4 && return Vec(2 * y * (4 * x + 4 * y - 3), -8 * x^2 - 8 * x * y + 12 * x + 6 * y - 4)
+    # Edge 3
+    i == 5 && return Vec(8 * x * y - 6 * x + 8 * y^2 - 12 * y + 4, 2 * x * (-4 * x - 4 * y + 3))
+    i == 6 && return Vec(
+        -8 * x * y + 6 * x + 2 * y - 2, 4 * x * (2 * x - 1)
+    )
+    # Face
+    i == 7 && return Vec(8 * y * (-x - 2 * y + 2), 8 * x * (x + 2 * y - 1))
+    i == 8 && return Vec(8 * y * (2 * x + y - 1), 8 * x * (-2 * x - y + 2))
+    throw(ArgumentError("no shape function $i for interpolation $ip"))
+end
+
+getnbasefunctions(::Nedelec{2, RefTriangle, 2}) = 8
+edgedof_interior_indices(::Nedelec{2, RefTriangle, 2}) = ((1, 2), (3, 4), (5, 6))
+facedof_interior_indices(::Nedelec{2, RefTriangle, 2}) = ((7, 8),)
+adjust_dofs_during_distribution(::Nedelec{2, RefTriangle, 2}) = true
+
+function get_direction(::Nedelec{2, RefTriangle, 2}, j, cell)
+    j > 6 && return 1
+    edge = edges(cell)[(j + 1) ÷ 2]
+    return ifelse(edge[2] > edge[1], 1, -1)
+end

src/iterators.jl

@@ -84,8 +84,14 @@ function reinit!(cc::CellCache, i::Int)
 end
 
 # reinit! FEValues with CellCache
-reinit!(cv::CellValues, cc::CellCache) = reinit!(cv, cc.coords)
-reinit!(fv::FacetValues, cc::CellCache, f::Int) = reinit!(fv, cc.coords, f)
+function reinit!(cv::CellValues, cc::CellCache)
+    cell = reinit_needs_cell(cv) ? getcells(cc.grid, cellid(cc)) : nothing
+    return reinit!(cv, cell, cc.coords)
+end
+function reinit!(fv::FacetValues, cc::CellCache, f::Int)
+    cell = reinit_needs_cell(fv) ? getcells(cc.grid, cellid(cc)) : nothing
+    return reinit!(fv, cell, cc.coords, f)
+end
 
 # Accessor functions
 getnodes(cc::CellCache) = cc.nodes

test/runtests.jl

@@ -11,6 +11,7 @@ using StaticArrays
 using OrderedCollections
 using WriteVTK
 import Metis
+using QuadGK: quadgk
 
 include("test_utils.jl")
 
@@ -39,6 +40,7 @@ include("test_apply_analytical.jl")
 include("PoolAllocator.jl")
 include("test_deprecations.jl")
 include("blockarrays.jl")
+include("test_continuity.jl")
 include("test_examples.jl")
 
 @test all(x -> isdefined(Ferrite, x), names(Ferrite))  # Test that all exported symbols are defined