Check quadrature point index before using at-inbounds.

Ferrite-FEM · Jan 18, 2022 · c22caa1 · c22caa1
1 parent 5994169
commit c22caa1
Showing 1 changed file with 12 additions and 0 deletions.
diff --git a/src/FEValues/common_values.jl b/src/FEValues/common_values.jl
@@ -11,6 +11,11 @@ getngeobasefunctions(cv::Values) = size(cv.M, 1)
 getn_scalarbasefunctions(cv::ScalarValues) = size(cv.N, 1)
 getn_scalarbasefunctions(cv::VectorValues{dim}) where {dim} = size(cv.N, 1) ÷ dim
 
+function checkquadpoint(cv::Values, qp::Int)
+    0 < qp <= getnquadpoints(cv) || error("quadrature point out of range")
+    return nothing
+end
+
 """
     reinit!(cv::CellValues, x::Vector)
     reinit!(bv::FaceValues, x::Vector, face::Int)
@@ -108,6 +113,7 @@ function function_value(fe_v::Values{dim}, q_point::Int, u::AbstractVector{T}, d
     isa(fe_v, VectorValues) && (n_base_funcs *= dim)
     @assert length(dof_range) == n_base_funcs
     @boundscheck checkbounds(u, dof_range)
+    @boundscheck checkquadpoint(fe_v, q_point)
     val = zero(_valuetype(fe_v, u))
     @inbounds for (i, j) in enumerate(dof_range)
         val += shape_value(fe_v, q_point, i) * u[j]
@@ -118,6 +124,7 @@ end
 function function_value(fe_v::VectorValues{dim}, q_point::Int, u::AbstractVector{Vec{dim,T}}) where {dim,T}
     n_base_funcs = getn_scalarbasefunctions(fe_v)
     @assert length(u) == n_base_funcs
+    @boundscheck checkquadpoint(fe_v, q_point)
     val = zero(Vec{dim, T})
     basefunc = 1
     @inbounds for i in 1:n_base_funcs
@@ -155,6 +162,7 @@ function function_gradient(fe_v::Values{dim}, q_point::Int, u::AbstractVector{T}
     isa(fe_v, VectorValues) && (n_base_funcs *= dim)
     @assert length(dof_range) == n_base_funcs
     @boundscheck checkbounds(u, dof_range)
+    @boundscheck checkquadpoint(fe_v, q_point)
     grad = zero(_gradienttype(fe_v, u))
     @inbounds for (i, j) in enumerate(dof_range)
         grad += shape_gradient(fe_v, q_point, i) * u[j]
@@ -168,6 +176,7 @@ Base.@pure _gradienttype(::VectorValues{dim}, ::AbstractVector{T}) where {dim,T}
 function function_gradient(fe_v::ScalarValues{dim}, q_point::Int, u::AbstractVector{Vec{dim,T}}) where {dim,T}
     n_base_funcs = getn_scalarbasefunctions(fe_v)
     @assert length(u) == n_base_funcs
+    @boundscheck checkquadpoint(fe_v, q_point)
     grad = zero(Tensor{2,dim,T})
     @inbounds for i in 1:n_base_funcs
         grad += u[i] ⊗ shape_gradient(fe_v, q_point, i)
@@ -178,6 +187,7 @@ end
 function function_gradient(fe_v::VectorValues{dim}, q_point::Int, u::AbstractVector{Vec{dim,T}}) where {dim,T}
     n_base_funcs = getn_scalarbasefunctions(fe_v)
     @assert length(u) == n_base_funcs
+    @boundscheck checkquadpoint(fe_v, q_point)
     grad = zero(Tensor{2,dim,T})
     basefunc_count = 1
     @inbounds for i in 1:n_base_funcs
@@ -224,6 +234,7 @@ where ``\\mathbf{u}_i`` are the nodal values of the function.
 function function_divergence(fe_v::ScalarValues{dim}, q_point::Int, u::AbstractVector{Vec{dim,T}}) where {dim,T}
     n_base_funcs = getn_scalarbasefunctions(fe_v)
     @assert length(u) == n_base_funcs
+    @boundscheck checkquadpoint(fe_v, q_point)
     diverg = zero(T)
     @inbounds for i in 1:n_base_funcs
         diverg += shape_gradient(fe_v, q_point, i) ⋅ u[i]
@@ -256,6 +267,7 @@ The coordinate is computed, using the geometric interpolation, as
 function spatial_coordinate(fe_v::Values{dim}, q_point::Int, x::AbstractVector{Vec{dim,T}}) where {dim,T}
     n_base_funcs = getngeobasefunctions(fe_v)
     @assert length(x) == n_base_funcs
+    @boundscheck checkquadpoint(fe_v, q_point)
     vec = zero(Vec{dim,T})
     @inbounds for i in 1:n_base_funcs
         vec += geometric_value(fe_v, q_point, i) * x[i]