Add polyorder kwarg, and test that quadrature rules fulfils this

src/Quadrature/gaussquad_tri_table.jl

@@ -108,7 +108,7 @@ function _get_dunavant_gauss_tridata(n::Int)
     return xw
 end
 
-# TheseqQuadrature rules for orders 9 to 20 have been obtained using the
+# These quadrature rules for orders 9 to 15 have been obtained using the
 # basix.make_quadrature(basix.CellType.triangle, n) calls of the
 # FEniCS / basix python package, which corresponds to Gauss-Jacobi rules.
 #

src/Quadrature/quadrature.jl

@@ -12,10 +12,14 @@ using Base.Cartesian: @nloops, @ntuple, @nexprs
 
 """
     QuadratureRule{shape}([::Type{T},] [quad_rule_type::Symbol,] order::Int)
+    QuadratureRule{shape}([::Type{T},] [quad_rule_type::Symbol,]; polyorder::Int)
     QuadratureRule{shape}(weights::AbstractVector{T}, points::AbstractVector{Vec{rdim, T}})
 
 Create a `QuadratureRule` used for integration on the refshape `shape` (of type [`AbstractRefShape`](@ref)).
-`order` is the order of the quadrature rule.
+`order` is the order of the quadrature rule, which for hypercubes corresponds to the number of quadrature points
+in each direction, and for simplices, prisms, and pyramids to the polynomial order that is fully integrated.
+To control the order of the polynomial that is exactly integrated for all reference shapes, the keyword argument
+`polyorder` can be passed instead of the `order` argument.
 `quad_rule_type` is an optional argument determining the type of quadrature rule,
 currently the `:legendre` and `:lobatto` rules are implemented for hypercubes.
 For triangles up to order 8 the default rule is the one by `:dunavant` (see [Dun:1985:hde](@cite)) and for
@@ -66,17 +70,40 @@ end
 @inline _default_quadrature_rule(::Type{RefTetrahedron}) = :keast_minimal
 
 # Fill in defaults with T=Float64
-function QuadratureRule{shape}(order::Int) where {shape <: AbstractRefShape}
-    return QuadratureRule{shape}(Float64, order)
+function QuadratureRule{shape}(args...; kwargs...) where {shape <: AbstractRefShape}
+    return QuadratureRule{shape}(Float64, args...; kwargs...)
 end
-function QuadratureRule{shape}(::Type{T}, order::Int) where {shape <: AbstractRefShape, T}
+function QuadratureRule{shape}(::Type{T}, args...; kwargs...) where {shape <: AbstractRefShape, T}
     quad_type = _default_quadrature_rule(shape)
+    return QuadratureRule{shape}(T, quad_type, args...)
+end
+function QuadratureRule{shape}(quad_type::Symbol, args...; kwargs...) where {shape <: AbstractRefShape}
+    return QuadratureRule{shape}(Float64, quad_type, args...; kwargs...)
+end
+
+# Passing polyorder kwarg:
+QuadratureRule{shape}(; kwargs...) where {shape <: AbstractRefShape} = QuadratureRule{shape}(Float64; kwargs...)
+function QuadratureRule{shape}(::Type{T}, quad_type::Symbol; polyorder::Int) where {shape <: AbstractRefShape, T}
+    order = _get_quadrature_order(shape, quad_type, polyorder)
     return QuadratureRule{shape}(T, quad_type, order)
 end
-function QuadratureRule{shape}(quad_type::Symbol, order::Int) where {shape <: AbstractRefShape}
-    return QuadratureRule{shape}(Float64, quad_type, order)
+
+function _get_quadrature_order(::Type{shape}, quad_type, polyorder) where {shape <: RefHypercube}
+    # Source: Pages pages 430 - 432 in
+    # Quarteroni, A., Sacco, R., Saleri, F. (2006). Orthogonal Polynomials in Approximation Theory.
+    # In: Numerical Mathematics. Texts in Applied Mathematics, vol 37. Springer, Berlin, Heidelberg.
+    # https://doi.org/10.1007/978-3-540-49809-4_10
+    if quad_type === :legendre
+        return ceil(Int, (polyorder + 1) / 2)
+    elseif quad_type === :lobatto
+        return ceil(Int, (polyorder + 3) / 2)
+    else
+        throw(ArgumentError("unsupported quadrature rule for $shape: $quad_type"))
+    end
 end
 
+_get_quadrature_order(::Type{shape}, quad_type, polyorder) where {shape <: Union{RefSimplex, RefPrism, RefPyramid}} = polyorder
+
 # Generate Gauss quadrature rules on hypercubes by doing an outer product
 # over all dimensions
 for dim in 1:3

test/test_quadrules.jl

@@ -237,4 +237,90 @@ using Ferrite: reference_shape_value
             end
         end
     end
+
+    @testset "polyorder" begin
+        # Define a struct that returns a complete polynomial function with random weights.
+        struct PolyFunction{dim}
+            weights::Vector{Float64}
+            exponents::Vector{NTuple{dim, Int}}
+        end
+
+        (f::PolyFunction{dim})(x::Vec{dim}) where {dim} = sum(w * sum(tuple(x...) .^ e) for (w, e) in zip(f.weights, f.exponents))
+
+        function PolyFunction{1}(order::Int)
+            return PolyFunction{1}(rand(order + 1), [(i,) for i in 0:order])
+        end
+        function PolyFunction{2}(order::Int)
+            exponents = NTuple{2, Int}[]
+            for i in 0:order
+                for j in 0:(order - i)
+                    push!(exponents, (i, j))
+                end
+            end
+            return PolyFunction{2}(rand(length(exponents)), exponents)
+        end
+        function PolyFunction{3}(order::Int)
+            exponents = NTuple{3, Int}[]
+            for i in 0:order
+                for j in 0:(order - i)
+                    for k in 0:(order - i - j)
+                        push!(exponents, (i, j, k))
+                    end
+                end
+            end
+            return PolyFunction{3}(rand(length(exponents)), exponents)
+        end
+
+        # Define functions to integrate over the difference reference shapes with hcubature.
+        function integration_reference(::Type{Ferrite.RefHypercube{dim}}, f::F; kwargs...) where {dim, F}
+            return hcubature(s -> f(Vec(s...)), -ones(Vec{dim}), ones(Vec{dim}); kwargs...)[1]
+        end
+        function integration_reference(::Type{RefTriangle}, f::F; kwargs...) where {F}
+            duffy_transform(s) = Vec((s[1] * (1 - s[2]), s[2]))
+            duffy_detJ(s) = 1 - s[2]
+            return hcubature(s -> f(duffy_transform(s)) * duffy_detJ(s), zero(Vec{2}), ones(Vec{2}); kwargs...)[1]
+        end
+        function integration_reference(::Type{RefPrism}, f::F; kwargs...) where {F}
+            duffy_transform(s) = Vec((s[1] * (1 - s[2]), s[2], s[3]))
+            duffy_detJ(s) = 1 - s[2]
+            return hcubature(s -> f(duffy_transform(s)) * duffy_detJ(s), zero(Vec{3}), ones(Vec{3}); kwargs...)[1]
+        end
+        function integration_reference(::Type{RefTetrahedron}, f::F; kwargs...) where {F}
+            duffy_transform(s) = Vec((s[1], (1 - s[1]) * s[2], (1 - s[1]) * (1 - s[2]) * s[3]))
+            duffy_detJ(s) = (1 - s[1])^2 * (1 - s[2])
+            return hcubature(s -> f(duffy_transform(s)) * duffy_detJ(s), zero(Vec{3}), ones(Vec{3}); kwargs...)[1]
+        end
+        function integration_reference(::Type{RefPyramid}, f::F; kwargs...) where {F}
+            duffy_transform(s) = Vec(s[1] * (1 - s[3]), s[2] * (1 - s[3]), s[3])
+            duffy_detJ(s) = (1 - s[3])^2
+            return hcubature(s -> f(duffy_transform(s)) * duffy_detJ(s), zero(Vec{3}), ones(Vec{3}); kwargs...)[1]
+        end
+
+        # Define function to integrate using the quadrature rule
+        function integration_check(f::F, qr::QuadratureRule) where {F}
+            s = zero(f(first(Ferrite.getpoints(qr))))
+            for (ξ, w) in zip(Ferrite.getpoints(qr), Ferrite.getweights(qr))
+                s += f(ξ) * w
+            end
+            return s
+        end
+
+        # Check that the polyorder kwarg gives the correct value
+        for ((shape, type), polyorders) in [
+                (RefLine, :legendre) => 1:10, (RefQuadrilateral, :legendre) => 1:3, (RefHexahedron, :legendre) => 1:3,
+                (RefLine, :lobatto) => 1:10, (RefQuadrilateral, :lobatto) => 1:3, (RefHexahedron, :lobatto) => 1:3,
+                (RefTriangle, :dunavant) => 1:8, (RefTriangle, :gaussjacobi) => 9:15,
+                (RefTetrahedron, :jinyun) => 1:3, (RefTetrahedron, :keast_minimal) => 4:5, (RefTetrahedron, :keast_positive) => 4,
+                (RefPrism, :polyquad) => 1:10, (RefPyramid, :polyquad) => 1:6,
+            ]
+            for polyorder in polyorders
+                @testset "QuadratureRule{$shape}($type; polyorder = $polyorder)" begin
+                    qr = QuadratureRule{shape}(type; polyorder)
+                    f = PolyFunction{Ferrite.getrefdim(shape)}(polyorder)
+                    solution = integration_reference(shape, f)
+                    @test solution ≈ integration_check(f, qr)
+                end
+            end
+        end
+    end
 end