Ensure that refshapes match when creating FEValues (#1185)

Author: KnutAM

src/FEValues/FunctionValues.jl

@@ -62,6 +62,7 @@ struct FunctionValues{DiffOrder, IP, N_t, dNdx_t, dNdξ_t, d2Ndx2_t, d2Ndξ2_t}
     end
 end
 function FunctionValues{DiffOrder}(::Type{T}, ip::Interpolation, qr::QuadratureRule, ip_geo::VectorizedInterpolation) where {DiffOrder, T}
+    assert_same_refshapes(qr, ip, ip_geo)
     n_shape = getnbasefunctions(ip)
     n_qpoints = getnquadpoints(qr)
 

src/FEValues/GeometryMapping.jl

@@ -55,13 +55,15 @@ struct GeometryMapping{DiffOrder, IP, M_t, dMdξ_t, d2Mdξ2_t}
     end
 end
 function GeometryMapping{0}(::Type{T}, ip::ScalarInterpolation, qr::QuadratureRule) where {T}
+    assert_same_refshapes(qr, ip, ip)
     n_shape = getnbasefunctions(ip)
     n_qpoints = getnquadpoints(qr)
     gm = GeometryMapping(ip, zeros(T, n_shape, n_qpoints), nothing, nothing)
     precompute_values!(gm, getpoints(qr))
     return gm
 end
 function GeometryMapping{1}(::Type{T}, ip::ScalarInterpolation, qr::QuadratureRule) where {T}
+    assert_same_refshapes(qr, ip, ip)
     n_shape = getnbasefunctions(ip)
     n_qpoints = getnquadpoints(qr)
 
@@ -73,6 +75,7 @@ function GeometryMapping{1}(::Type{T}, ip::ScalarInterpolation, qr::QuadratureRu
     return gm
 end
 function GeometryMapping{2}(::Type{T}, ip::ScalarInterpolation, qr::QuadratureRule) where {T}
+    assert_same_refshapes(qr, ip, ip)
     n_shape = getnbasefunctions(ip)
     n_qpoints = getnquadpoints(qr)
 

src/FEValues/common_values.jl

@@ -394,3 +394,16 @@ function reference_shape_hessians_gradients_and_values!(hessians::AbstractMatrix
     end
     return
 end
+
+assert_same_refshapes(::Union{QuadratureRule{RS}, FacetQuadratureRule{RS}}, ::Interpolation{RS}, ::Interpolation{RS}) where {RS} = nothing
+function assert_same_refshapes(qr::Union{QuadratureRule, FacetQuadratureRule}, ipf::Interpolation, ipg::Interpolation)
+    throw(
+        ArgumentError(
+            """
+            The reference shapes of the quadrature rule ($(getrefshape(qr))),
+            function interpolation ($(getrefshape(ipf))), and geometric interpolation ($(getrefshape(ipg))),
+            must be equal.
+            """
+        )
+    )
+end

test/test_cellvalues.jl

@@ -478,6 +478,13 @@
         end
     end
 
+    @testset "construction errors" begin
+        @test_throws ArgumentError CellValues(QuadratureRule{RefTriangle}(1), Lagrange{RefQuadrilateral, 1}())
+        @test_throws ArgumentError CellValues(QuadratureRule{RefTriangle}(1), Lagrange{RefTriangle, 1}(), Lagrange{RefQuadrilateral, 1}())
+        @test_throws ArgumentError CellValues(QuadratureRule{RefTriangle}(1), Lagrange{RefQuadrilateral, 1}(), Lagrange{RefQuadrilateral, 1}())
+        @test_throws ArgumentError CellValues(QuadratureRule{RefTriangle}(1), Lagrange{RefQuadrilateral, 1}(), Lagrange{RefTriangle, 1}())
+    end
+
     @testset "show" begin
         cv_quad = CellValues(QuadratureRule{RefQuadrilateral}(2), Lagrange{RefQuadrilateral, 2}()^2)
         showstring = sprint(show, MIME"text/plain"(), cv_quad)

test/test_facevalues.jl

@@ -177,6 +177,13 @@
         end
     end
 
+    @testset "construction errors" begin
+        @test_throws ArgumentError FacetValues(FacetQuadratureRule{RefTriangle}(1), Lagrange{RefQuadrilateral, 1}())
+        @test_throws ArgumentError FacetValues(FacetQuadratureRule{RefTriangle}(1), Lagrange{RefTriangle, 1}(), Lagrange{RefQuadrilateral, 1}())
+        @test_throws ArgumentError FacetValues(FacetQuadratureRule{RefTriangle}(1), Lagrange{RefQuadrilateral, 1}(), Lagrange{RefQuadrilateral, 1}())
+        @test_throws ArgumentError FacetValues(FacetQuadratureRule{RefTriangle}(1), Lagrange{RefQuadrilateral, 1}(), Lagrange{RefTriangle, 1}())
+    end
+
     @testset "show" begin
         # Just smoke test to make sure show doesn't error.
         fv = FacetValues(FacetQuadratureRule{RefQuadrilateral}(2), Lagrange{RefQuadrilateral, 2}())

test/test_interfacevalues.jl

@@ -155,6 +155,13 @@
             test_interfacevalues(grid, iv)
         end
     end
+    @testset "construction errors" begin
+        @test_throws ArgumentError InterfaceValues(FacetQuadratureRule{RefTriangle}(1), Lagrange{RefQuadrilateral, 1}())
+        @test_throws ArgumentError InterfaceValues(FacetQuadratureRule{RefTriangle}(1), Lagrange{RefTriangle, 1}(), Lagrange{RefQuadrilateral, 1}())
+        @test_throws ArgumentError InterfaceValues(FacetQuadratureRule{RefTriangle}(1), Lagrange{RefQuadrilateral, 1}(), Lagrange{RefQuadrilateral, 1}())
+        @test_throws ArgumentError InterfaceValues(FacetQuadratureRule{RefTriangle}(1), Lagrange{RefQuadrilateral, 1}(), Lagrange{RefTriangle, 1}())
+        @test_throws ArgumentError InterfaceValues(FacetQuadratureRule{RefTriangle}(1), Lagrange{RefTriangle, 1}(), FacetQuadratureRule{RefTriangle}(1), Lagrange{RefQuadrilateral, 1}())
+    end
     # Custom quadrature
     @testset "Custom quadrature interface values" begin
         cell_shape = Tetrahedron