Skip to content

Commit

Permalink
Fold Ferrite.value into methods of shape_value
Browse files Browse the repository at this point in the history 
This patch deprecates `Ferrite.value(ip::Interpolation, i::Int, ξ::Vec)`
in favor of `shape_value(ip::Interpolation, ξ::Vec, i::Int)`. Note that
the argument order is flipped to resemble the `(Cell|Face)Values` API
more. Compare e.g. `shape_value(ip, ξ, i)` with `shape_value(cv, qp,
i)`: the second argument becomes the "quadrature point" argument (`ξ` ~
`qp`), and the third argument is the shape function index, in both
cases.

Since `shape_value` is exported this introduces new (public)
functionality to users -- for some applications it is nice to be able to
evaluate the shape functions in the reference domain without having to
go through `CellValues`. For the same reason this patch also adds the
method `shape_gradient(ip::Interpolation, ξ::Vec, i::Int)` which falls
back to AD from Tensors.jl.

Fixes #609.
  • Loading branch information
@fredrikekre
fredrikekre committed May 22, 2023
1 parent 240cd94 commit 30a144c
Show file tree
Hide file tree
Showing 8 changed files with 77 additions and 64 deletions.
4 changes: 2 additions & 2 deletions src/FEValues/cell_values.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -69,10 +69,10 @@ function CellValues{IP, N_t, dNdx_t, dNdξ_t, T, dMdξ_t, QR, GIP}(qr::QR, ip::I

for (qp, ξ) in pairs(getpoints(qr))
for basefunc in 1:n_func_basefuncs
dNdξ[basefunc, qp], N[basefunc, qp] = gradient_and_value(ip, basefunc, ξ)
dNdξ[basefunc, qp], N[basefunc, qp] = shape_gradient_and_value(ip, ξ, basefunc)
end
for basefunc in 1:n_geom_basefuncs
dMdξ[basefunc, qp], M[basefunc, qp] = gradient_and_value(gip, basefunc, ξ)
dMdξ[basefunc, qp], M[basefunc, qp] = shape_gradient_and_value(gip, ξ, basefunc)
end
end

Expand Down
6 changes: 3 additions & 3 deletions src/FEValues/face_values.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -72,10 +72,10 @@ function FaceValues{IP, N_t, dNdx_t, dNdξ_t, T, dMdξ_t, QR, Normal_t, GIP}(qr:

for face in 1:n_faces, (qp, ξ) in pairs(getpoints(qr, face))
for basefunc in 1:n_func_basefuncs
dNdξ[basefunc, qp, face], N[basefunc, qp, face] = gradient_and_value(ip, basefunc, ξ)
dNdξ[basefunc, qp, face], N[basefunc, qp, face] = shape_gradient_and_value(ip, ξ, basefunc)
end
for basefunc in 1:n_geom_basefuncs
dMdξ[basefunc, qp, face], M[basefunc, qp, face] = gradient_and_value(gip, basefunc, ξ)
dMdξ[basefunc, qp, face], M[basefunc, qp, face] = shape_gradient_and_value(gip, ξ, basefunc)
end
end

Expand Down Expand Up @@ -212,7 +212,7 @@ function BCValues(::Type{T}, func_interpol::Interpolation{refshape}, geom_interp

for n_boundary_entity in 1:n_boundary_entities
for (qp, ξ) in enumerate(qrs[n_boundary_entity].points), i in 1:n_geom_basefuncs
M[i, qp, n_boundary_entity] = value(geom_interpol, i, ξ)
M[i, qp, n_boundary_entity] = shape_value(geom_interpol, ξ, i)
end
nqp[n_boundary_entity] = length(qrs[n_boundary_entity].points)
end
Expand Down
2 changes: 1 addition & 1 deletion src/PointEval/PointEvalHandler.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -126,7 +126,7 @@ function find_local_coordinate(interpolation, cell_coordinates::Vector{V}, globa
global_guess = zero(V)
J = zero(Tensor{2, dim, T})
for j in 1:n_basefuncs
dNdξ, N = gradient_and_value(interpolation, j, local_guess)
dNdξ, N = shape_gradient_and_value(interpolation, local_guess, j)
global_guess += N * cell_coordinates[j]
J += cell_coordinates[j] ⊗ dNdξ
end
Expand Down
6 changes: 3 additions & 3 deletions src/PointEval/point_values.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -60,7 +60,7 @@ function reinit!(pv::PointValues, x::AbstractVector{<:Vec{D}}, ξ::Vec{D}) where
qp = 1 # PointValues only have a single qp
# TODO: Does M need to be updated too?
for i in 1:getnbasefunctions(pv.cv.ip)
pv.cv.dNdξ[i, qp], pv.cv.N[i, qp] = gradient_and_value(pv.cv.ip, i, ξ)
pv.cv.dNdξ[i, qp], pv.cv.N[i, qp] = shape_gradient_and_value(pv.cv.ip, ξ, i)
end
reinit!(pv.cv, x)
return nothing
Expand All @@ -76,7 +76,7 @@ end

function PointValuesInternal(ξ::Vec{dim, T}, ip::IP) where {dim, T, shape <: AbstractRefShape{dim}, IP <: Interpolation{shape}}
n_func_basefuncs = getnbasefunctions(ip)
N = [value(ip, i, ξ) for i in 1:n_func_basefuncs]
N = [shape_value(ip, ξ, i) for i in 1:n_func_basefuncs]
return PointValuesInternal{IP, eltype(N)}(N, ip)
end

Expand All @@ -88,7 +88,7 @@ shape_value(pv::PointValuesInternal, qp::Int, i::Int) = (@assert qp == 1; pv.N[i
function reinit!(pv::PointValuesInternal{IP}, coord::Vec{dim}) where {dim, shape <: AbstractRefShape{dim}, IP <: Interpolation{shape}}
n_func_basefuncs = getnbasefunctions(pv.ip)
for i in 1:n_func_basefuncs
pv.N[i] = value(pv.ip, i, coord)
pv.N[i] = shape_value(pv.ip, coord, i)
end
return nothing
end
5 changes: 3 additions & 2 deletions src/deprecations.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -326,5 +326,6 @@ function Base.show(io::IO, ::CrouzeixRaviart{shape, order}) where {shape, order}
print(io, "CrouzeixRaviart{$(shape), $(order)}()")
end

@deprecate value(ip::Interpolation, ξ::Vec) [Ferrite.value(ip, i, ξ) for i in 1:getnbasefunctions(ip)] false
@deprecate derivative(ip::Interpolation, ξ::Vec) [Tensors.gradient(x -> Ferrite.value(ip, i, x), ξ) for i in 1:getnbasefunctions(ip)] false
@deprecate value(ip::Interpolation, ξ::Vec) [shape_value(ip, ξ, i) for i in 1:getnbasefunctions(ip)] false
@deprecate derivative(ip::Interpolation, ξ::Vec) [shape_gradient(ip, ξ, i) for i in 1:getnbasefunctions(ip)] false
@deprecate value(ip::Interpolation, i::Int, ξ::Vec) shape_value(ip, ξ, i) false
101 changes: 56 additions & 45 deletions src/interpolations.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -177,10 +177,6 @@ Return order of the interpolation.
@inline getorder(::Interpolation{shape,order}) where {shape,order} = order


function gradient_and_value(ip::Interpolation, i::Int, x::Vec)
return gradient(ξ -> value(ip, i, ξ), x, :all)
end

#####################
# Utility functions #
#####################
Expand All @@ -199,17 +195,32 @@ getnbasefunctions(::Interpolation)
# celldof: dof that is local to the element

"""
value(ip::Interpolation, i::Int, ξ::Vec)
shape_value(ip::Interpolation, ξ::Vec, i::Int)
Evaluates the `i`'th basis function of the interpolation `ip`
at a point `ξ` on the reference element. The index `i` must
Evaluate the value of the `i`th shape function of the interpolation `ip` in
at a point `ξ` on the reference element. The index `i` must
match the index in [`vertices(::Interpolation)`](@ref), [`faces(::Interpolation)`](@ref) and
[`edges(::Interpolation)`](@ref).
For nodal interpolations the indices also must match the
For nodal interpolations the indices also must match the
indices of [`reference_coordinates(::Interpolation)`](@ref).
"""
value(ip::Interpolation, i::Int, ξ::Vec)
shape_value(ip::Interpolation, ξ::Vec, i::Int)

"""
shape_gradient(ip::Interpolation, ξ::Vec, i::Int)
Evaluate the gradient of the `i`th shape function of the interpolation `ip` in
reference coordinate `ξ`.
"""
function shape_gradient(ip::Interpolation, ξ::Vec, i::Int)
return Tensors.gradient(x -> shape_value(ip, x, i), ξ)
end

function shape_gradient_and_value(ip::Interpolation, x::Vec, i::Int)
return gradient(ξ -> shape_value(ip, ξ, i), x, :all)
end


"""
reference_coordinates(ip::Interpolation)
Expand Down Expand Up @@ -339,8 +350,8 @@ celldof_interior_indices(ip::DiscontinuousLagrange) = ntuple(i->i, getnbasefunct
function reference_coordinates(ip::DiscontinuousLagrange{shape, order}) where {shape, order}
return reference_coordinates(Lagrange{shape,order}())
end
function value(ip::DiscontinuousLagrange{shape, order}, i::Int, ξ::Vec{dim}) where {dim, shape<:AbstractRefShape{dim}, order}
return value(Lagrange{shape, order}(), i, ξ)
function shape_value(::DiscontinuousLagrange{shape, order}, ξ::Vec{dim}, i::Int) where {dim, shape <: AbstractRefShape{dim}, order}
return shape_value(Lagrange{shape, order}(), ξ, i)
end

# Excepting the L0 element.
Expand All @@ -356,9 +367,9 @@ function reference_coordinates(ip::DiscontinuousLagrange{RefTetrahedron,0})
return [Vec{3,Float64}((1/4,1/4,1/4))]
end

function value(ip::DiscontinuousLagrange{shape,0}, i::Int, ξ::Vec{dim}) where {dim,shape<:AbstractRefShape{dim}}
function shape_value(ip::DiscontinuousLagrange{shape, 0}, ::Vec{dim, T}, i::Int) where {dim, shape <: AbstractRefShape{dim}, T}
i > 1 && throw(ArgumentError("no shape function $i for interpolation $ip"))
return 1.0
return one(T)
end

############
Expand Down Expand Up @@ -391,7 +402,7 @@ function reference_coordinates(::Lagrange{RefLine,1})
Vec{1, Float64}(( 1.0,))]
end

function value(ip::Lagrange{RefLine,1}, i::Int, ξ::Vec{1})
function shape_value(ip::Lagrange{RefLine, 1}, ξ::Vec{1}, i::Int)
ξ_x = ξ[1]
i == 1 && return (1 - ξ_x) * 0.5
i == 2 && return (1 + ξ_x) * 0.5
Expand All @@ -412,7 +423,7 @@ function reference_coordinates(::Lagrange{RefLine,2})
Vec{1, Float64}(( 0.0,))]
end

function value(ip::Lagrange{RefLine,2}, i::Int, ξ::Vec{1})
function shape_value(ip::Lagrange{RefLine, 2}, ξ::Vec{1}, i::Int)
ξ_x = ξ[1]
i == 1 && return ξ_x * (ξ_x - 1) * 0.5
i == 2 && return ξ_x * (ξ_x + 1) * 0.5
Expand All @@ -434,7 +445,7 @@ function reference_coordinates(::Lagrange{RefQuadrilateral,1})
Vec{2, Float64}((-1.0, 1.0,))]
end

function value(ip::Lagrange{RefQuadrilateral,1}, i::Int, ξ::Vec{2})
function shape_value(ip::Lagrange{RefQuadrilateral, 1}, ξ::Vec{2}, i::Int)
ξ_x = ξ[1]
ξ_y = ξ[2]
i == 1 && return (1 - ξ_x) * (1 - ξ_y) * 0.25
Expand Down Expand Up @@ -465,7 +476,7 @@ function reference_coordinates(::Lagrange{RefQuadrilateral,2})
Vec{2, Float64}(( 0.0, 0.0))]
end

function value(ip::Lagrange{RefQuadrilateral,2}, i::Int, ξ::Vec{2})
function shape_value(ip::Lagrange{RefQuadrilateral, 2}, ξ::Vec{2}, i::Int)
ξ_x = ξ[1]
ξ_y = ξ[2]
i == 1 && return (ξ_x^2 - ξ_x) * (ξ_y^2 - ξ_y) * 0.25
Expand Down Expand Up @@ -493,7 +504,7 @@ function reference_coordinates(::Lagrange{RefTriangle,1})
Vec{2, Float64}((0.0, 0.0))]
end

function value(ip::Lagrange{RefTriangle,1}, i::Int, ξ::Vec{2})
function shape_value(ip::Lagrange{RefTriangle, 1}, ξ::Vec{2}, i::Int)
ξ_x = ξ[1]
ξ_y = ξ[2]
i == 1 && return ξ_x
Expand All @@ -519,7 +530,7 @@ function reference_coordinates(::Lagrange{RefTriangle,2})
Vec{2, Float64}((0.5, 0.0))]
end

function value(ip::Lagrange{RefTriangle,2}, i::Int, ξ::Vec{2})
function shape_value(ip::Lagrange{RefTriangle, 2}, ξ::Vec{2}, i::Int)
ξ_x = ξ[1]
ξ_y = ξ[2]
γ = 1. - ξ_x - ξ_y
Expand Down Expand Up @@ -598,7 +609,7 @@ function reference_coordinates(ip::Lagrange2Tri345)
return permute!(coordpts, permdof2DLagrange2Tri345[order])
end

function value(ip::Lagrange2Tri345, i::Int, ξ::Vec{2})
function shape_value(ip::Lagrange2Tri345, ξ::Vec{2}, i::Int)
if !(0 < i <= getnbasefunctions(ip))
throw(ArgumentError("no shape function $i for interpolation $ip"))
end
Expand Down Expand Up @@ -644,7 +655,7 @@ function reference_coordinates(::Lagrange{RefTetrahedron,1})
Vec{3, Float64}((0.0, 0.0, 1.0))]
end

function value(ip::Lagrange{RefTetrahedron,1}, i::Int, ξ::Vec{3})
function shape_value(ip::Lagrange{RefTetrahedron, 1}, ξ::Vec{3}, i::Int)
ξ_x = ξ[1]
ξ_y = ξ[2]
ξ_z = ξ[3]
Expand Down Expand Up @@ -679,7 +690,7 @@ end

# http://www.colorado.edu/engineering/CAS/courses.d/AFEM.d/AFEM.Ch09.d/AFEM.Ch09.pdf
# http://www.colorado.edu/engineering/CAS/courses.d/AFEM.d/AFEM.Ch10.d/AFEM.Ch10.pdf
function value(ip::Lagrange{RefTetrahedron,2}, i::Int, ξ::Vec{3})
function shape_value(ip::Lagrange{RefTetrahedron, 2}, ξ::Vec{3}, i::Int)
ξ_x = ξ[1]
ξ_y = ξ[2]
ξ_z = ξ[3]
Expand Down Expand Up @@ -715,7 +726,7 @@ function reference_coordinates(::Lagrange{RefHexahedron,1})
Vec{3, Float64}((-1.0, 1.0, 1.0))]
end

function value(ip::Lagrange{RefHexahedron,1}, i::Int, ξ::Vec{3})
function shape_value(ip::Lagrange{RefHexahedron, 1}, ξ::Vec{3}, i::Int)
ξ_x = ξ[1]
ξ_y = ξ[2]
ξ_z = ξ[3]
Expand Down Expand Up @@ -803,7 +814,7 @@ function reference_coordinates(::Lagrange{RefHexahedron,2})
]
end

function value(ip::Lagrange{RefHexahedron,2}, i::Int, ξ::Vec{3, T}) where {T}
function shape_value(ip::Lagrange{RefHexahedron, 2}, ξ::Vec{3, T}, i::Int) where {T}
# Some local helpers.
@inline φ₁(x::T) = -0.5*x*(1-x)
@inline φ₂(x::T) = (1+x)*(1-x)
Expand Down Expand Up @@ -862,7 +873,7 @@ function reference_coordinates(::Lagrange{RefPrism,1})
Vec{3, Float64}((0.0, 1.0, 1.0))]
end

function value(ip::Lagrange{RefPrism,1}, i::Int, ξ::Vec{3})
function shape_value(ip::Lagrange{RefPrism,1}, ξ::Vec{3}, i::Int)
(x,y,z) = ξ
i == 1 && return 1-x-y -z*(1-x-y)
i == 2 && return x*(1-z)
Expand Down Expand Up @@ -942,7 +953,7 @@ function reference_coordinates(::Lagrange{RefPrism,2})
Vec{3, Float64}((1/2, 1/2, 1/2)),]
end

function value(ip::Lagrange{RefPrism,2}, i::Int, ξ::Vec{3})
function shape_value(ip::Lagrange{RefPrism, 2}, ξ::Vec{3}, i::Int)
(x,y,z) = ξ
x² = x*x
y² = y*y
Expand Down Expand Up @@ -997,7 +1008,7 @@ function reference_coordinates(::BubbleEnrichedLagrange{RefTriangle,1})
Vec{2, Float64}((1/3, 1/3)),]
end

function value(ip::BubbleEnrichedLagrange{RefTriangle,1}, i::Int, ξ::Vec{2})
function shape_value(ip::BubbleEnrichedLagrange{RefTriangle, 1}, ξ::Vec{2}, i::Int)
ξ_x = ξ[1]
ξ_y = ξ[2]
i == 1 && return ξ_x*(9ξ_y^2 + 9ξ_x*ξ_y - 9ξ_y + 1)
Expand Down Expand Up @@ -1040,7 +1051,7 @@ function reference_coordinates(::Serendipity{RefQuadrilateral,2})
Vec{2, Float64}((-1.0, 0.0))]
end

function value(ip::Serendipity{RefQuadrilateral,2}, i::Int, ξ::Vec{2})
function shape_value(ip::Serendipity{RefQuadrilateral,2}, ξ::Vec{2}, i::Int)
ξ_x = ξ[1]
ξ_y = ξ[2]
i == 1 && return (1 - ξ_x) * (1 - ξ_y) * 0.25(-ξ_x - ξ_y - 1)
Expand Down Expand Up @@ -1111,18 +1122,18 @@ function reference_coordinates(::Serendipity{RefHexahedron,2})
Vec{3, Float64}((-1.0, 1.0, 0.0)),]
end

function value(ip::Serendipity{RefHexahedron,2}, i::Int, ξ::Vec{3})
function shape_value(ip::Serendipity{RefHexahedron, 2}, ξ::Vec{3}, i::Int)
ξ_x = ξ[1]
ξ_y = ξ[2]
ξ_z = ξ[3]
i == 1 && return 0.125(1 - ξ_x) * (1 - ξ_y) * (1 - ξ_z) - 0.5(value(ip,12,ξ) + value(ip,9,ξ) + value(ip,17,ξ))
i == 2 && return 0.125(1 + ξ_x) * (1 - ξ_y) * (1 - ξ_z) - 0.5(value(ip,9,ξ) + value(ip,10,ξ) + value(ip,18,ξ))
i == 3 && return 0.125(1 + ξ_x) * (1 + ξ_y) * (1 - ξ_z) - 0.5(value(ip,10,ξ) + value(ip,11,ξ) + value(ip,19,ξ))
i == 4 && return 0.125(1 - ξ_x) * (1 + ξ_y) * (1 - ξ_z) - 0.5(value(ip,11,ξ) + value(ip,12,ξ) + value(ip,20,ξ))
i == 5 && return 0.125(1 - ξ_x) * (1 - ξ_y) * (1 + ξ_z) - 0.5(value(ip,16,ξ) + value(ip,13,ξ) + value(ip,17,ξ))
i == 6 && return 0.125(1 + ξ_x) * (1 - ξ_y) * (1 + ξ_z) - 0.5(value(ip,13,ξ) + value(ip,14,ξ) + value(ip,18,ξ))
i == 7 && return 0.125(1 + ξ_x) * (1 + ξ_y) * (1 + ξ_z) - 0.5(value(ip,14,ξ) + value(ip,15,ξ) + value(ip,19,ξ))
i == 8 && return 0.125(1 - ξ_x) * (1 + ξ_y) * (1 + ξ_z) - 0.5(value(ip,15,ξ) + value(ip,16,ξ) + value(ip,20,ξ))
i == 1 && return 0.125(1 - ξ_x) * (1 - ξ_y) * (1 - ξ_z) - 0.5(shape_value(ip, ξ, 12) + shape_value(ip, ξ, 9) + shape_value(ip, ξ, 17))
i == 2 && return 0.125(1 + ξ_x) * (1 - ξ_y) * (1 - ξ_z) - 0.5(shape_value(ip, ξ, 9) + shape_value(ip, ξ, 10) + shape_value(ip, ξ, 18))
i == 3 && return 0.125(1 + ξ_x) * (1 + ξ_y) * (1 - ξ_z) - 0.5(shape_value(ip, ξ, 10) + shape_value(ip, ξ, 11) + shape_value(ip, ξ, 19))
i == 4 && return 0.125(1 - ξ_x) * (1 + ξ_y) * (1 - ξ_z) - 0.5(shape_value(ip, ξ, 11) + shape_value(ip, ξ, 12) + shape_value(ip, ξ, 20))
i == 5 && return 0.125(1 - ξ_x) * (1 - ξ_y) * (1 + ξ_z) - 0.5(shape_value(ip, ξ, 16) + shape_value(ip, ξ, 13) + shape_value(ip, ξ, 17))
i == 6 && return 0.125(1 + ξ_x) * (1 - ξ_y) * (1 + ξ_z) - 0.5(shape_value(ip, ξ, 13) + shape_value(ip, ξ, 14) + shape_value(ip, ξ, 18))
i == 7 && return 0.125(1 + ξ_x) * (1 + ξ_y) * (1 + ξ_z) - 0.5(shape_value(ip, ξ, 14) + shape_value(ip, ξ, 15) + shape_value(ip, ξ, 19))
i == 8 && return 0.125(1 - ξ_x) * (1 + ξ_y) * (1 + ξ_z) - 0.5(shape_value(ip, ξ, 15) + shape_value(ip, ξ, 16) + shape_value(ip, ξ, 20))
i == 9 && return 0.25(1 - ξ_x^2) * (1 - ξ_y) * (1 - ξ_z)
i == 10 && return 0.25(1 + ξ_x) * (1 - ξ_y^2) * (1 - ξ_z)
i == 11 && return 0.25(1 - ξ_x^2) * (1 + ξ_y) * (1 - ξ_z)
Expand Down Expand Up @@ -1165,7 +1176,7 @@ function reference_coordinates(::CrouzeixRaviart)
Vec{2, Float64}((0.5, 0.0))]
end

function value(ip::CrouzeixRaviart, i::Int, ξ::Vec{2})
function shape_value(ip::CrouzeixRaviart, ξ::Vec{2}, i::Int)
ξ_x = ξ[1]
ξ_y = ξ[2]
i == 1 && return 2*ξ_x + 2*ξ_y - 1.0
Expand Down Expand Up @@ -1201,24 +1212,24 @@ get_n_copies(::VectorizedInterpolation{vdim}) where vdim = vdim
function getnbasefunctions(ipv::VectorizedInterpolation{vdim}) where vdim
return vdim * getnbasefunctions(ipv.ip)
end
function value(ipv::VectorizedInterpolation{vdim, shape}, I::Int, ξ::Vec{refdim, T}) where {vdim, refdim, shape <: AbstractRefShape{refdim}, T}
function shape_value(ipv::VectorizedInterpolation{vdim, shape}, ξ::Vec{refdim, T}, I::Int) where {vdim, refdim, shape <: AbstractRefShape{refdim}, T}
i0, c0 = divrem(I - 1, vdim)
i = i0 + 1
c = c0 + 1
v = value(ipv.ip, i, ξ)
v = shape_value(ipv.ip, ξ, i)
return Vec{vdim, T}(j -> j == c ? v : zero(v))
end

# vdim == refdim
function gradient_and_value(ipv::VectorizedInterpolation{dim, shape}, I::Int, ξ::Vec{dim}) where {dim, shape <: AbstractRefShape{dim}}
return invoke(gradient_and_value, Tuple{Interpolation, Int, Vec}, ipv, I, ξ)
function shape_gradient_and_value(ipv::VectorizedInterpolation{dim, shape}, ξ::Vec{dim}, I::Int) where {dim, shape <: AbstractRefShape{dim}}
return invoke(shape_gradient_and_value, Tuple{Interpolation, Vec, Int}, ipv, ξ, I)
end
# vdim != refdim
function gradient_and_value(ipv::VectorizedInterpolation{vdim, shape}, I::Int, ξ::V) where {vdim, refdim, shape <: AbstractRefShape{refdim}, T, V <: Vec{refdim, T}}
function shape_gradient_and_value(ipv::VectorizedInterpolation{vdim, shape}, ξ::V, I::Int) where {vdim, refdim, shape <: AbstractRefShape{refdim}, T, V <: Vec{refdim, T}}
# Load with dual numbers and compute the value
f = x -> value(ipv, I, x)
f = x -> shape_value(ipv, x, I)
ξd = Tensors._load(ξ, Tensors.Tag(f, V))
value_grad = value(ipv, I, ξd)
value_grad = f(ξd)
# Extract the value and gradient
val = Vec{vdim, T}(i -> Tensors.value(value_grad[i]))
grad = zero(MMatrix{vdim, refdim, T})
Expand Down
5 changes: 3 additions & 2 deletions test/test_deprecations.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -83,8 +83,9 @@ end
@testset "Ferrite.value and Ferrite.derivative" begin
ip = Lagrange{RefQuadrilateral, 1}()
ξ = zero(Vec{2})
@test (@test_deprecated Ferrite.value(ip, ξ)) == [Ferrite.value(ip, i, ξ) for i in 1:getnbasefunctions(ip)]
@test (@test_deprecated Ferrite.derivative(ip, ξ)) == [Tensors.gradient(x -> Ferrite.value(ip, i, x), ξ) for i in 1:getnbasefunctions(ip)]
@test (@test_deprecated Ferrite.value(ip, ξ)) == [shape_value(ip, ξ, i) for i in 1:getnbasefunctions(ip)]
@test (@test_deprecated Ferrite.derivative(ip, ξ)) == [shape_gradient(ip, ξ, i) for i in 1:getnbasefunctions(ip)]
@test (@test_deprecated Ferrite.value(ip, 1, ξ)) == shape_value(ip, ξ, 1)
end

end # testset deprecations
12 changes: 6 additions & 6 deletions test/test_interpolations.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,16 +47,16 @@
# Check partition of unity at random point.
n_basefuncs = getnbasefunctions(interpolation)
x = rand(Tensor{1, ref_dim})
f = (x) -> [Ferrite.value(interpolation, i, Tensor{1, ref_dim}(x)) for i in 1:n_basefuncs]
f = (x) -> [shape_value(interpolation, Tensor{1, ref_dim}(x), i) for i in 1:n_basefuncs]
@test vec(ForwardDiff.jacobian(f, Array(x))') ≈
reinterpret(Float64, [Tensors.gradient(y -> Ferrite.value(interpolation, i, y), x) for i in 1:n_basefuncs])
@test sum([Ferrite.value(interpolation, i, x) for i in 1:n_basefuncs]) ≈ 1.0
reinterpret(Float64, [shape_gradient(interpolation, x, i) for i in 1:n_basefuncs])
@test sum([shape_value(interpolation, x, i) for i in 1:n_basefuncs]) ≈ 1.0

# Check if the important functions are consistent
coords = Ferrite.reference_coordinates(interpolation)
@test length(coords) == n_basefuncs
@test Ferrite.value(interpolation, length(coords), x) == Ferrite.value(interpolation, length(coords), x)
@test_throws ArgumentError Ferrite.value(interpolation, length(coords)+1, x)
@test shape_value(interpolation, x, n_basefuncs) == shape_value(interpolation, x, n_basefuncs)
@test_throws ArgumentError shape_value(interpolation, x, n_basefuncs+1)

# Test whether we have for each entity corresponding dof indices (possibly empty)
@test length(Ferrite.vertexdof_indices(interpolation)) == Ferrite.nvertices(interpolation)
Expand Down Expand Up @@ -98,7 +98,7 @@
# Check for dirac delta property of interpolation
@testset "dirac delta property of dof $dof" for dof in 1:n_basefuncs
for k in 1:n_basefuncs
N_dof = Ferrite.value(interpolation, k, coords[dof])
N_dof = shape_value(interpolation, coords[dof], k)
if k == dof
@test N_dof ≈ 1.0
else
Expand Down

0 comments on commit 30a144c

Please sign in to comment.