Add missing dispatches to _mass_qr

src/L2_projection.jl

@@ -101,7 +101,7 @@
   projection equation, when calling [`project`](@ref). It should match the quadrature points used when
   creating the quadrature-point variables to project.
 * The *optional* `qr_lhs` sets the quadrature rule used to integrate the left-hand-side of the projection equation,
-  and defaults to a quadrature rule that integrates the mass-matrix exactly for the given interpolation `ip`.
+  and defaults to a quadrature rule that integrates the mass-matrix accurately for the given interpolation `ip`.
 
 """
 function add!(
@@ -137,13 +137,12 @@
 end
 
 # Quadrature sufficient for integrating a mass matrix
-function _mass_qr(::Lagrange{shape, order}) where {shape <: AbstractRefShape, order}
-    return QuadratureRule{shape}(order + 1)
-end
-function _mass_qr(::Lagrange{shape, 2}) where {shape <: RefSimplex}
-    return QuadratureRule{shape}(4)
-end
-_mass_qr(ip::VectorizedInterpolation) = _mass_qr(ip.ip)
+_mass_qr(ip::Interpolation{RefShape}) where {RefShape} = QuadratureRule{RefShape}(_mass_qr_order(ip))
+_mass_qr_order(::Lagrange{<:AbstractRefShape, order}) where {order} = order + 1
+_mass_qr_order(::DiscontinuousLagrange{RefShape, order}) where {RefShape, order} = _mass_qr_order(Lagrange{RefShape, order}())
+_mass_qr_order(::Serendipity{RefShape, order}) where {RefShape, order} = _mass_qr_order(Lagrange{RefShape, order}())
+_mass_qr_order(::Lagrange{RefSimplex{refdim}, order}) where {order, refdim} = max(1, refdim * order)
+_mass_qr_order(ip::VectorizedInterpolation) = _mass_qr_order(ip.ip)
 
 function _assemble_L2_matrix(dh::DofHandler, qrs_lhs::Vector{<:QuadratureRule})
     M = Symmetric(allocate_matrix(dh))

test/test_l2_projection.jl

@@ -493,6 +493,58 @@ function test_l2proj_errorpaths()
     return
 end
 
+function test_mass_qr()
+    function _ips(shape)
+        ips = (
+            Lagrange{shape, 1}(), Lagrange{shape, 2}(),
+            Lagrange{shape, 3}(), DiscontinuousLagrange{shape, 0}(),
+            DiscontinuousLagrange{shape, 1}(), DiscontinuousLagrange{shape, 2}(), DiscontinuousLagrange{shape, 3}(),
+        )
+        if shape <: Ferrite.RefHypercube
+            ips = (ips..., Serendipity{shape, 2}())
+        end
+        return ips
+    end
+
+    function _mass_qr_test(ip::Interpolation{RefShape}) where {RefShape}
+        return QuadratureRule{RefShape}(Ferrite._mass_qr_order(ip) + 1)
+    end
+
+    function _calculate_local_me(grid, ip::Interpolation, qr::QuadratureRule)
+        cv = CellValues(qr, ip)
+        cc = CellCache(grid)
+        Ferrite.reinit!(cc, 1)
+        Ferrite.reinit!(cv, cc)
+
+        nbasef = getnbasefunctions(cv)
+        Me = zeros(nbasef, nbasef)
+        for q_point in 1:getnquadpoints(cv)
+            dV = getdetJdV(cv, q_point)
+            for i in 1:nbasef
+                Ni = shape_value(cv, q_point, i)
+                for j in 1:nbasef
+                    Nj = shape_value(cv, q_point, j)
+                    Me[i, j] += (Ni ⋅ Nj) * dV
+                end
+            end
+        end
+        return Me
+    end
+
+    for (shape, refshape) in ((Quadrilateral, RefQuadrilateral), (Triangle, RefTriangle))
+        grid = generate_grid(shape, (2, 2), Vec{2}((0.0, 0.0)), Vec{2}((1.0, 1.0)))
+        for ip in _ips(refshape)
+            me_1 = _calculate_local_me(grid, ip, Ferrite._mass_qr(ip))
+            me_2 = _calculate_local_me(grid, ip, _mass_qr_test(ip))
+            @test me_1 ≈ me_2
+            if !(me_1 ≈ me_2)
+                println("The mass matrix calculated with order+1 is not identical to the mass matrix calculated with order for $ip")
+            end
+        end
+    end
+    return
+end
+
 @testset "Test L2-Projection" begin
     test_projection(1, RefQuadrilateral)
     test_projection(1, RefTriangle)
@@ -505,4 +557,5 @@ end
     test_export(subset = true)
     test_show_l2()
     test_l2proj_errorpaths()
+    test_mass_qr()
 end