Store interpolations and quadrature rules in FEValues.

docs/src/literate/linear_shell.jl

@@ -304,11 +304,11 @@ function integrate_shell!(ke, cv, qr_ooplane, X, data)
             B = ForwardDiff.jacobian(
                 (a) -> strain(a, N, dNdx, ζ, dζdx, q, ef1, ef2, h), zeros(Float64, ndofs) )
 
-            dV = det(J) * cv.qr_weights[iqp] * qr_ooplane.weights[oqp]
+            dV = det(J) * cv.qr.weights[iqp] * qr_ooplane.weights[oqp]
             ke .+= B'*data.C*B * dV
         end
     end
 end;
 
 # Run everything:
-main()
+main()

src/FEValues/cell_values.jl

@@ -41,7 +41,11 @@ struct CellScalarValues{dim,T<:Real,refshape<:AbstractRefShape} <: CellValues{di
     detJdV::Vector{T}
     M::Matrix{T}
     dMdξ::Matrix{Vec{dim,T}}
-    qr_weights::Vector{T}
+    qr::QuadratureRule{dim,refshape,T}
+    # The following fields are deliberately abstract -- they are never used in
+    # performance critical code, just stored here for convenience.
+    func_interp::Interpolation{dim,refshape}
+    geo_interp::Interpolation{dim,refshape}
 end
 
 function CellScalarValues(quad_rule::QuadratureRule, func_interpol::Interpolation,
@@ -78,7 +82,7 @@ function CellScalarValues(::Type{T}, quad_rule::QuadratureRule{dim,shape}, func_
 
     detJdV = fill(T(NaN), n_qpoints)
 
-    CellScalarValues{dim,T,shape}(N, dNdx, dNdξ, detJdV, M, dMdξ, quad_rule.weights)
+    CellScalarValues{dim,T,shape}(N, dNdx, dNdξ, detJdV, M, dMdξ, quad_rule, func_interpol, geom_interpol)
 end
 
 # CellVectorValues
@@ -89,7 +93,11 @@ struct CellVectorValues{dim,T<:Real,refshape<:AbstractRefShape,M} <: CellValues{
     detJdV::Vector{T}
     M::Matrix{T}
     dMdξ::Matrix{Vec{dim,T}}
-    qr_weights::Vector{T}
+    qr::QuadratureRule{dim,refshape,T}
+    # The following fields are deliberately abstract -- they are never used in
+    # performance critical code, just stored here for convenience.
+    func_interp::Interpolation{dim,refshape}
+    geo_interp::Interpolation{dim,refshape}
 end
 
 function CellVectorValues(quad_rule::QuadratureRule, func_interpol::Interpolation, geom_interpol::Interpolation=func_interpol)
@@ -137,22 +145,22 @@ function CellVectorValues(::Type{T}, quad_rule::QuadratureRule{dim,shape}, func_
     detJdV = fill(T(NaN), n_qpoints)
     MM = Tensors.n_components(Tensors.get_base(eltype(dNdx)))
 
-    CellVectorValues{dim,T,shape,MM}(N, dNdx, dNdξ, detJdV, M, dMdξ, quad_rule.weights)
+    CellVectorValues{dim,T,shape,MM}(N, dNdx, dNdξ, detJdV, M, dMdξ, quad_rule, func_interpol, geom_interpol)
 end
 
 function reinit!(cv::CellValues{dim}, x::AbstractVector{Vec{dim,T}}) where {dim,T}
     n_geom_basefuncs = getngeobasefunctions(cv)
     n_func_basefuncs = getnbasefunctions(cv)
     @assert length(x) == n_geom_basefuncs
 
-    @inbounds for i in 1:length(cv.qr_weights)
-        w = cv.qr_weights[i]
+    @inbounds for i in 1:length(cv.qr.weights)
+        w = cv.qr.weights[i]
         fecv_J = zero(Tensor{2,dim})
         for j in 1:n_geom_basefuncs
             fecv_J += x[j] ⊗ cv.dMdξ[j, i]
         end
         detJ = det(fecv_J)
-        detJ > 0.0 || throw(ArgumentError("det(J) is not positive: det(J) = $(detJ)"))
+        detJ > 0.0 || throw_detJ_not_pos(detJ)
         cv.detJdV[i] = detJ * w
         Jinv = inv(fecv_J)
         for j in 1:n_func_basefuncs

src/FEValues/common_values.jl

@@ -12,6 +12,7 @@ FieldTrait(::Type{<:PointScalarValues}) = ScalarValued()
 FieldTrait(::Type{<:CellVectorValues}) = VectorValued()
 FieldTrait(::Type{<:FaceVectorValues}) = VectorValued()
 
+@noinline throw_detJ_not_pos(detJ) = throw(ArgumentError("det(J) is not positive: det(J) = $(detJ)"))
 
 getnbasefunctions(cv::Values) = size(cv.N, 1)
 getngeobasefunctions(cv::Values) = size(cv.M, 1)
@@ -40,7 +41,7 @@ reinit!
 
 Return the number of quadrature points for the `Values` object.
 """
-getnquadpoints(fe::Values) = length(fe.qr_weights)
+getnquadpoints(fe::Values) = length(fe.qr.weights)
 
 """
     getdetJdV(fe_v::Values, q_point::Int)

src/FEValues/face_values.jl

@@ -48,8 +48,16 @@ struct FaceScalarValues{dim,T<:Real,refshape<:AbstractRefShape} <: FaceValues{di
     normals::Vector{Vec{dim,T}}
     M::Array{T,3}
     dMdξ::Array{Vec{dim,T},3}
-    qr_weights::Vector{T}
+    # 'Any' is 'dim-1' here -- this is deliberately abstractly typed. Only qr.weights is
+    # accessed in performance critical code so this doesn't seem to be a problem in
+    # practice since qr.weights is correctly inferred as Vector{T}, and T is a parameter
+    # of the struct.
+    qr::QuadratureRule{<:Any,refshape,T}
     current_face::ScalarWrapper{Int}
+    # The following fields are deliberately abstract -- they are never used in
+    # performance critical code, just stored here for convenience.
+    func_interp::Interpolation{dim,refshape}
+    geo_interp::Interpolation{dim,refshape}
 end
 
 function FaceScalarValues(quad_rule::QuadratureRule, func_interpol::Interpolation,
@@ -93,7 +101,7 @@ function FaceScalarValues(::Type{T}, quad_rule::QuadratureRule{dim_qr,shape}, fu
 
     detJdV = fill(T(NaN), n_qpoints, n_faces)
 
-    FaceScalarValues{dim,T,shape}(N, dNdx, dNdξ, detJdV, normals, M, dMdξ, quad_rule.weights, ScalarWrapper(0))
+    FaceScalarValues{dim,T,shape}(N, dNdx, dNdξ, detJdV, normals, M, dMdξ, quad_rule, ScalarWrapper(0), func_interpol, geom_interpol)
 end
 
 # FaceVectorValues
@@ -105,8 +113,16 @@ struct FaceVectorValues{dim,T<:Real,refshape<:AbstractRefShape,M} <: FaceValues{
     normals::Vector{Vec{dim,T}}
     M::Array{T,3}
     dMdξ::Array{Vec{dim,T},3}
-    qr_weights::Vector{T}
+    # 'Any' is 'dim-1' here -- this is deliberately abstractly typed. Only qr.weights is
+    # accessed in performance critical code so this doesn't seem to be a problem in
+    # practice since qr.weights is correctly inferred as Vector{T}, and T is a parameter
+    # of the struct.
+    qr::QuadratureRule{<:Any,refshape,T}
     current_face::ScalarWrapper{Int}
+    # The following fields are deliberately abstract -- they are never used in
+    # performance critical code, just stored here for convenience.
+    func_interp::Interpolation{dim,refshape}
+    geo_interp::Interpolation{dim,refshape}
 end
 
 function FaceVectorValues(quad_rule::QuadratureRule, func_interpol::Interpolation, geom_interpol::Interpolation=func_interpol)
@@ -161,7 +177,7 @@ function FaceVectorValues(::Type{T}, quad_rule::QuadratureRule{dim_qr,shape}, fu
     detJdV = fill(T(NaN), n_qpoints, n_faces)
     MM = Tensors.n_components(Tensors.get_base(eltype(dNdx)))
 
-    FaceVectorValues{dim,T,shape,MM}(N, dNdx, dNdξ, detJdV, normals, M, dMdξ, quad_rule.weights, ScalarWrapper(0))
+    FaceVectorValues{dim,T,shape,MM}(N, dNdx, dNdξ, detJdV, normals, M, dMdξ, quad_rule, ScalarWrapper(0), func_interpol, geom_interpol)
 end
 
 function reinit!(fv::FaceValues{dim}, x::AbstractVector{Vec{dim,T}}, face::Int) where {dim,T}
@@ -173,8 +189,8 @@ function reinit!(fv::FaceValues{dim}, x::AbstractVector{Vec{dim,T}}, face::Int)
     fv.current_face[] = face
     cb = getcurrentface(fv)
 
-    @inbounds for i in 1:length(fv.qr_weights)
-        w = fv.qr_weights[i]
+    @inbounds for i in 1:length(fv.qr.weights)
+        w = fv.qr.weights[i]
         fefv_J = zero(Tensor{2,dim})
         for j in 1:n_geom_basefuncs
             fefv_J += x[j] ⊗ fv.dMdξ[j, i, cb]
@@ -183,7 +199,7 @@ function reinit!(fv::FaceValues{dim}, x::AbstractVector{Vec{dim,T}}, face::Int)
         fv.normals[i] = weight_norm / norm(weight_norm)
         detJ = norm(weight_norm)
 
-        detJ > 0.0 || throw(ArgumentError("det(J) is not positive: det(J) = $(detJ)"))
+        detJ > 0.0 || throw_detJ_not_pos(detJ)
         fv.detJdV[i, cb] = detJ * w
         Jinv = inv(fefv_J)
         for j in 1:n_func_basefuncs

src/Quadrature/quadrature.jl

@@ -40,6 +40,8 @@ struct QuadratureRule{dim,shape,T}
     points::Vector{Vec{dim,T}}
 end
 
+Base.copy(qr::QuadratureRule) = qr # TODO: Is it ever useful to get an actual copy?
+
 """
     getweights(qr::QuadratureRule)
 

src/interpolations.jl

@@ -32,6 +32,8 @@ julia> getnbasefunctions(ip)
 """
 abstract type Interpolation{dim,shape,order} end
 
+Base.copy(ip::Interpolation) = ip
+
 """
 Return the dimension of an `Interpolation`
 """
@@ -693,4 +695,4 @@ function value(ip::CrouzeixRaviart{2,1}, i::Int, ξ::Vec{2})
     i == 2 && return 1.0 - 2*ξ_x
     i == 3 && return 1.0 - 2*ξ_y
     throw(ArgumentError("no shape function $i for interpolation $ip"))
-end
+end

test/test_cellvalues.jl

@@ -76,10 +76,14 @@ for (func_interpol, quad_rule) in  (
 
         # test copy
         cvc = copy(cv)
+        @test typeof(cv) == typeof(cvc)
         for fname in fieldnames(typeof(cv))
-            @test typeof(cv) == typeof(cvc)
-            @test pointer(getfield(cv, fname)) != pointer(getfield(cvc, fname))
-            @test getfield(cv, fname) == getfield(cvc, fname)
+            v = getfield(cv, fname)
+            vc = getfield(cvc, fname)
+            if hasmethod(pointer, Tuple{typeof(v)})
+                @test pointer(getfield(cv, fname)) != pointer(getfield(cvc, fname))
+            end
+            @test v == vc
         end
     end
 end

test/test_facevalues.jl

@@ -83,12 +83,15 @@ for (func_interpol, quad_rule) in  (
 
         # test copy
         fvc = copy(fv)
+        @test typeof(fv) == typeof(fvc)
         for fname in fieldnames(typeof(fv))
-            @test typeof(fv) == typeof(fvc)
-            if !isa(getfield(fv, fname), Ferrite.ScalarWrapper)
-                @test pointer(getfield(fv, fname)) != pointer(getfield(fvc, fname))
-                @test getfield(fv, fname) == getfield(fvc, fname)
+            v = getfield(fv, fname)
+            v isa Ferrite.ScalarWrapper && continue
+            vc = getfield(fvc, fname)
+            if hasmethod(pointer, Tuple{typeof(v)})
+                @test pointer(v) != pointer(vc)
             end
+            @test v == vc
         end
     end
 end