From 4d7c140f77bfda9b7dadfe2446aa419c2e3a7248 Mon Sep 17 00:00:00 2001 From: Keita Nakamura Date: Sat, 27 Jul 2024 11:04:40 +0900 Subject: [PATCH] Add SteffenBSpline --- src/Interpolations/bspline.jl | 283 ++++++++++++++----------- src/Interpolations/gimp.jl | 2 +- src/Interpolations/kernelcorrection.jl | 14 +- src/Interpolations/mpvalue.jl | 16 +- src/Tesserae.jl | 4 + test/interpolations.jl | 45 ++-- 6 files changed, 207 insertions(+), 157 deletions(-) diff --git a/src/Interpolations/bspline.jl b/src/Interpolations/bspline.jl index 8627dfd9..9dd375af 100644 --- a/src/Interpolations/bspline.jl +++ b/src/Interpolations/bspline.jl @@ -1,42 +1,10 @@ -struct BSpline{order} <: Kernel - BSpline{order}() where {order} = new{order::Int}() -end - -""" - BSpline{1}() - LinearBSpline() - -Linear B-spline kernel. -""" -const LinearBSpline = BSpline{1} +abstract type AbstractBSpline{order} <: Kernel end -""" - BSpline{2}() - QuadraticBSpline() +gridspan(::AbstractBSpline{1}) = 2 +gridspan(::AbstractBSpline{2}) = 3 +gridspan(::AbstractBSpline{3}) = 4 -Quadratic B-spline kernel. -The peaks of this funciton are centered on the grid nodes [^Steffen]. - -[^Steffen]: [Steffen, M., Kirby, R. M., & Berzins, M. (2008). Analysis and reduction of quadrature errors in the material point method (MPM). *International journal for numerical methods in engineering*, 76(6), 922-948.](https://doi.org/10.1002/nme.2360) -""" -const QuadraticBSpline = BSpline{2} - -""" - BSpline{3}() - CubicBSpline() - -Cubic B-spline kernel. -The peaks of this funciton are centered on the grid nodes [^Steffen]. - -[^Steffen]: [Steffen, M., Kirby, R. M., & Berzins, M. (2008). Analysis and reduction of quadrature errors in the material point method (MPM). *International journal for numerical methods in engineering*, 76(6), 922-948.](https://doi.org/10.1002/nme.2360) -""" -const CubicBSpline = BSpline{3} - -gridspan(::BSpline{1}) = 2 -gridspan(::BSpline{2}) = 3 -gridspan(::BSpline{3}) = 4 - -@inline function neighboringnodes(spline::BSpline, pt, mesh::CartesianMesh{dim}) where {dim} +@inline function neighboringnodes(spline::AbstractBSpline, pt, mesh::CartesianMesh{dim}) where {dim} x = getx(pt) ξ = Tuple(normalize(x, mesh)) dims = size(mesh) @@ -49,103 +17,31 @@ gridspan(::BSpline{3}) = 4 imax = Tuple(@. min(stop, dims)) CartesianIndices(UnitRange.(imin, imax)) end -@inline _neighboringnodes_offset(::BSpline{1}) = 0.0 -@inline _neighboringnodes_offset(::BSpline{2}) = 0.5 -@inline _neighboringnodes_offset(::BSpline{3}) = 1.0 +@inline _neighboringnodes_offset(::AbstractBSpline{1}) = 0.0 +@inline _neighboringnodes_offset(::AbstractBSpline{2}) = 0.5 +@inline _neighboringnodes_offset(::AbstractBSpline{3}) = 1.0 # simple B-spline calculations -function value(::BSpline{1}, ξ::Real) +function value(::AbstractBSpline{1}, ξ::Real) ξ = abs(ξ) ξ < 1 ? 1 - ξ : zero(ξ) end -function value(::BSpline{2}, ξ::Real) +function value(::AbstractBSpline{2}, ξ::Real) ξ = abs(ξ) ξ < 0.5 ? (3 - 4ξ^2) / 4 : ξ < 1.5 ? (3 - 2ξ)^2 / 8 : zero(ξ) end -function value(::BSpline{3}, ξ::Real) +function value(::AbstractBSpline{3}, ξ::Real) ξ = abs(ξ) ξ < 1 ? (3ξ^3 - 6ξ^2 + 4) / 6 : ξ < 2 ? (2 - ξ)^3 / 6 : zero(ξ) end -@generated function value(spline::BSpline, ξ::Vec{dim}) where {dim} - quote - @_inline_meta - prod(@ntuple $dim i -> value(spline, ξ[i])) - end -end -@inline function value(spline::BSpline, pt, mesh::CartesianMesh, i) - @_propagate_inbounds_meta - xₚ = getx(pt) - xᵢ = mesh[i] - h⁻¹ = spacing_inv(mesh) - ξ = (xₚ - xᵢ) * h⁻¹ - value(spline, ξ) -end - -# Steffen, M., Kirby, R. M., & Berzins, M. (2008). -# Analysis and reduction of quadrature errors in the material point method (MPM). -# International journal for numerical methods in engineering, 76(6), 922-948. -function value(spline::BSpline{1}, ξ::Real, pos::Int)::typeof(ξ) - value(spline, ξ) -end -function value(spline::BSpline{2}, ξ::Real, pos::Int)::typeof(ξ) - if pos == 0 - ξ = abs(ξ) - ξ < 0.5 ? (3 - 4ξ^2) / 3 : - ξ < 1.5 ? (3 - 2ξ)^2 / 6 : zero(ξ) - elseif abs(pos) == 1 - ξ = sign(pos) * ξ - ξ < -1 ? zero(ξ) : - ξ < -0.5 ? 4(1 + ξ)^2 / 3 : - ξ < 0.5 ? -(28ξ^2 - 4ξ - 17) / 24 : - ξ < 1.5 ? (3 - 2ξ)^2 / 8 : zero(ξ) - else - value(spline, ξ) - end -end -function value(spline::BSpline{3}, ξ::Real, pos::Int)::typeof(ξ) - if pos == 0 - ξ = abs(ξ) - ξ < 1 ? (3ξ^3 - 6ξ^2 + 4) / 4 : - ξ < 2 ? (2 - ξ)^3 / 4 : zero(ξ) - elseif abs(pos) == 1 - ξ = sign(pos) * ξ - ξ < -1 ? zero(ξ) : - ξ < 0 ? (1 + ξ)^2 * (7 - 11ξ) / 12 : - ξ < 1 ? (7ξ^3 - 15ξ^2 + 3ξ + 7) / 12 : - ξ < 2 ? (2 - ξ)^3 / 6 : zero(ξ) - else - value(spline, ξ) - end -end -@generated function value(spline::BSpline, ξ::Vec{dim}, pos::Vec{dim}) where {dim} - quote - @_inline_meta - prod(@ntuple $dim i -> value(spline, ξ[i], pos[i])) - end -end -@inline function value(spline::BSpline, xₚ::Vec, mesh::CartesianMesh, i::CartesianIndex, ::Symbol) # last argument is pseudo argument `:steffen` - @_propagate_inbounds_meta - xᵢ = mesh[i] - h⁻¹ = spacing_inv(mesh) - ξ = (xₚ - xᵢ) * h⁻¹ - value(spline, ξ, node_position(mesh, i)) -end - -@inline function node_position(ax::AbstractVector, i::Int) - left = i - firstindex(ax) - right = lastindex(ax) - i - ifelse(left < right, left, -right) -end -node_position(mesh::CartesianMesh, index::CartesianIndex) = Vec(map(node_position, get_axes(mesh), Tuple(index))) - @inline fract(x) = x - floor(x) # Fast calculations for value, gradient and hessian # `x` must be normalized by `h` -@inline Base.values(spline::BSpline, x, args...) = only(values(identity, spline, x, args...)) -@inline function Base.values(diff, ::BSpline{1}, x::Real) +@inline Base.values(spline::AbstractBSpline, x, args...) = only(values(identity, spline, x, args...)) +@inline function Base.values(diff, ::AbstractBSpline{1}, x::Real) T = typeof(x) ξ = fract(x) vals = tuple(@. T((1-ξ, ξ))) @@ -160,7 +56,7 @@ node_position(mesh::CartesianMesh, index::CartesianIndex) = Vec(map(node_positio end vals end -@inline function Base.values(diff, ::BSpline{2}, x::Real) +@inline function Base.values(diff, ::AbstractBSpline{2}, x::Real) T = typeof(x) x′ = fract(x - T(0.5)) ξ = @. x′ - T((-0.5,0.5,1.5)) @@ -176,7 +72,7 @@ end end vals end -@inline function Base.values(diff, ::BSpline{3}, x::Real) +@inline function Base.values(diff, ::AbstractBSpline{3}, x::Real) T = typeof(x) x′ = fract(x) ξ = @. x′ - T((-1,0,1,2)) @@ -195,7 +91,7 @@ end vals end -@generated function Base.values(diff, spline::BSpline, x::Vec{dim}) where {dim} +@generated function Base.values(diff, spline::AbstractBSpline, x::Vec{dim}) where {dim} T_∇∇ws = SymmetricSecondOrderTensor{dim} ∇∇ws = Array{Expr}(undef,dim,dim) for j in 1:dim, i in 1:dim @@ -241,38 +137,171 @@ end end @inline tuple_otimes(x::Tuple) = SArray(otimes(map(Vec, x)...)) -@inline function Base.values(::typeof(identity), spline::BSpline, x::Vec, mesh::CartesianMesh) +@inline function Base.values(::typeof(identity), spline::AbstractBSpline, x::Vec, mesh::CartesianMesh) h⁻¹ = spacing_inv(mesh) (values(spline, (x - first(mesh)) * h⁻¹),) end -@inline function Base.values(::typeof(gradient), spline::BSpline, x::Vec, mesh::CartesianMesh) +@inline function Base.values(::typeof(gradient), spline::AbstractBSpline, x::Vec, mesh::CartesianMesh) xmin = first(mesh) h⁻¹ = spacing_inv(mesh) w, ∇w = values(gradient, spline, (x-xmin)*h⁻¹) (w, ∇w*h⁻¹) end -@inline function Base.values(::typeof(hessian), spline::BSpline, x::Vec, mesh::CartesianMesh) +@inline function Base.values(::typeof(hessian), spline::AbstractBSpline, x::Vec, mesh::CartesianMesh) xmin = first(mesh) h⁻¹ = spacing_inv(mesh) w, ∇w, ∇∇w = values(hessian, spline, (x-xmin)*h⁻¹) (w, ∇w*h⁻¹, ∇∇w*h⁻¹*h⁻¹) end -@inline function Base.values(::typeof(all), spline::BSpline, x::Vec, mesh::CartesianMesh) +@inline function Base.values(::typeof(all), spline::AbstractBSpline, x::Vec, mesh::CartesianMesh) xmin = first(mesh) h⁻¹ = spacing_inv(mesh) w, ∇w, ∇∇w, ∇∇∇w = values(all, spline, (x-xmin)*h⁻¹) (w, ∇w*h⁻¹, ∇∇w*h⁻¹*h⁻¹, ∇∇∇w*h⁻¹*h⁻¹*h⁻¹) end -function update_property!(mp::MPValue, it::BSpline, pt, mesh::CartesianMesh) +function update_property!(mp::MPValue, it::AbstractBSpline, pt, mesh::CartesianMesh) indices = neighboringnodes(mp) isnearbounds = size(mp.w) != size(indices) if isnearbounds @inbounds for ip in eachindex(indices) i = indices[ip] - set_shape_values!(mp, ip, value(difftype(mp), it, getx(pt), mesh, i, :steffen)) + set_kernel_values!(mp, ip, value(difftype(mp), it, getx(pt), mesh, i)) end else - set_shape_values!(mp, values(difftype(mp), it, getx(pt), mesh)) + set_kernel_values!(mp, values(difftype(mp), it, getx(pt), mesh)) end end + +struct BSpline{order} <: AbstractBSpline{order} + BSpline{order}() where {order} = new{order::Int}() +end + +""" + BSpline{1}() + LinearBSpline() + +Linear B-spline kernel. +""" +const LinearBSpline = BSpline{1} + +""" + BSpline{2}() + QuadraticBSpline() + +Quadratic B-spline kernel. +""" +const QuadraticBSpline = BSpline{2} + +""" + BSpline{3}() + CubicBSpline() + +Cubic B-spline kernel. +""" +const CubicBSpline = BSpline{3} + +@generated function value(spline::BSpline, ξ::Vec{dim}) where {dim} + quote + @_inline_meta + prod(@ntuple $dim i -> value(spline, ξ[i])) + end +end +@inline function value(spline::BSpline, pt, mesh::CartesianMesh, i) + @_propagate_inbounds_meta + xₚ = getx(pt) + xᵢ = mesh[i] + h⁻¹ = spacing_inv(mesh) + ξ = (xₚ - xᵢ) * h⁻¹ + value(spline, ξ) +end + +struct SteffenBSpline{order} <: AbstractBSpline{order} + SteffenBSpline{order}() where {order} = new{order::Int}() +end + +""" + SteffenBSpline{1}() + SteffenLinearBSpline() + +Linear B-spline kernel. +""" +const SteffenLinearBSpline = SteffenBSpline{1} + +""" + SteffenBSpline{2}() + SteffenQuadraticBSpline() + +Quadratic B-spline kernel proposed by [^Steffen]. + +[^Steffen]: [Steffen, M., Kirby, R. M., & Berzins, M. (2008). Analysis and reduction of quadrature errors in the material point method (MPM). *International journal for numerical methods in engineering*, 76(6), 922-948.](https://doi.org/10.1002/nme.2360) +""" +const SteffenQuadraticBSpline = SteffenBSpline{2} + +""" + SteffenBSpline{3}() + SteffenCubicBSpline() + +Cubic B-spline kernel proposed by [^Steffen]. + +[^Steffen]: [Steffen, M., Kirby, R. M., & Berzins, M. (2008). Analysis and reduction of quadrature errors in the material point method (MPM). *International journal for numerical methods in engineering*, 76(6), 922-948.](https://doi.org/10.1002/nme.2360) +""" +const SteffenCubicBSpline = SteffenBSpline{3} + +# Steffen, M., Kirby, R. M., & Berzins, M. (2008). +# Analysis and reduction of quadrature errors in the material point method (MPM). +# International journal for numerical methods in engineering, 76(6), 922-948. +function value(spline::SteffenBSpline{1}, ξ::Real, pos::Int)::typeof(ξ) + value(spline, ξ) +end +function value(spline::SteffenBSpline{2}, ξ::Real, pos::Int)::typeof(ξ) + if pos == 0 + ξ = abs(ξ) + ξ < 0.5 ? (3 - 4ξ^2) / 3 : + ξ < 1.5 ? (3 - 2ξ)^2 / 6 : zero(ξ) + elseif abs(pos) == 1 + ξ = sign(pos) * ξ + ξ < -1 ? zero(ξ) : + ξ < -0.5 ? 4(1 + ξ)^2 / 3 : + ξ < 0.5 ? -(28ξ^2 - 4ξ - 17) / 24 : + ξ < 1.5 ? (3 - 2ξ)^2 / 8 : zero(ξ) + else + value(spline, ξ) + end +end +function value(spline::SteffenBSpline{3}, ξ::Real, pos::Int)::typeof(ξ) + if pos == 0 + ξ = abs(ξ) + ξ < 1 ? (3ξ^3 - 6ξ^2 + 4) / 4 : + ξ < 2 ? (2 - ξ)^3 / 4 : zero(ξ) + elseif abs(pos) == 1 + ξ = sign(pos) * ξ + ξ < -1 ? zero(ξ) : + ξ < 0 ? (1 + ξ)^2 * (7 - 11ξ) / 12 : + ξ < 1 ? (7ξ^3 - 15ξ^2 + 3ξ + 7) / 12 : + ξ < 2 ? (2 - ξ)^3 / 6 : zero(ξ) + else + value(spline, ξ) + end +end +@generated function value(spline::SteffenBSpline, ξ::Vec{dim}, pos::Vec{dim}) where {dim} + quote + @_inline_meta + prod(@ntuple $dim i -> value(spline, ξ[i], pos[i])) + end +end +@inline function value(spline::SteffenBSpline, pt, mesh::CartesianMesh, i::CartesianIndex) + @_propagate_inbounds_meta + xₚ = getx(pt) + xᵢ = mesh[i] + h⁻¹ = spacing_inv(mesh) + ξ = (xₚ - xᵢ) * h⁻¹ + value(spline, ξ, node_position(mesh, i)) +end + +@inline function node_position(ax::AbstractVector, i::Int) + left = i - firstindex(ax) + right = lastindex(ax) - i + ifelse(left < right, left, -right) +end +node_position(mesh::CartesianMesh, index::CartesianIndex) = Vec(map(node_position, get_axes(mesh), Tuple(index))) diff --git a/src/Interpolations/gimp.jl b/src/Interpolations/gimp.jl index 302fe5ef..23fe77e0 100644 --- a/src/Interpolations/gimp.jl +++ b/src/Interpolations/gimp.jl @@ -50,6 +50,6 @@ end indices = neighboringnodes(mp) @inbounds @simd for ip in eachindex(indices) i = indices[ip] - set_shape_values!(mp, ip, value(difftype(mp), it, pt, mesh, i)) + set_kernel_values!(mp, ip, value(difftype(mp), it, pt, mesh, i)) end end diff --git a/src/Interpolations/kernelcorrection.jl b/src/Interpolations/kernelcorrection.jl index f855a2b6..578df72b 100644 --- a/src/Interpolations/kernelcorrection.jl +++ b/src/Interpolations/kernelcorrection.jl @@ -27,19 +27,19 @@ gridspan(kc::KernelCorrection) = gridspan(get_kernel(kc)) else @inbounds @simd for ip in eachindex(indices) i = indices[ip] - set_shape_values!(mp, ip, value(difftype(mp), get_kernel(it), pt, mesh, i)) + set_kernel_values!(mp, ip, value(difftype(mp), get_kernel(it), pt, mesh, i)) end end end # fast version for B-spline kernels -@inline function update_property!(mp::MPValue, it::KernelCorrection{<: BSpline}, pt, mesh::CartesianMesh, filter::AbstractArray{Bool} = Trues(size(mesh))) +@inline function update_property!(mp::MPValue, it::KernelCorrection{<: AbstractBSpline}, pt, mesh::CartesianMesh, filter::AbstractArray{Bool} = Trues(size(mesh))) indices = neighboringnodes(mp) isnearbounds = size(mp.w) != size(indices) || !alltrue(filter, indices) if isnearbounds update_property_nearbounds!(mp, it, pt, mesh, filter) else - set_shape_values!(mp, values(difftype(mp), get_kernel(it), getx(pt), mesh)) + set_kernel_values!(mp, values(difftype(mp), get_kernel(it), getx(pt), mesh)) end end @@ -67,10 +67,10 @@ end w = mp.w[ip] P = value(poly, xᵢ - xₚ) wq = w * (M⁻¹ ⋅ P) - hasproperty(mp, :w) && set_shape_values!(mp, ip, (wq⋅P₀,)) - hasproperty(mp, :∇w) && set_shape_values!(mp, ip, (wq⋅P₀, wq⋅∇P₀)) - hasproperty(mp, :∇∇w) && set_shape_values!(mp, ip, (wq⋅P₀, wq⋅∇P₀, wq⋅∇∇P₀)) - hasproperty(mp, :∇∇∇w) && set_shape_values!(mp, ip, (wq⋅P₀, wq⋅∇P₀, wq⋅∇∇P₀, wq⋅∇∇∇P₀)) + hasproperty(mp, :w) && set_kernel_values!(mp, ip, (wq⋅P₀,)) + hasproperty(mp, :∇w) && set_kernel_values!(mp, ip, (wq⋅P₀, wq⋅∇P₀)) + hasproperty(mp, :∇∇w) && set_kernel_values!(mp, ip, (wq⋅P₀, wq⋅∇P₀, wq⋅∇∇P₀)) + hasproperty(mp, :∇∇∇w) && set_kernel_values!(mp, ip, (wq⋅P₀, wq⋅∇P₀, wq⋅∇∇P₀, wq⋅∇∇∇P₀)) end end diff --git a/src/Interpolations/mpvalue.jl b/src/Interpolations/mpvalue.jl index 9aa17598..5e9a4bb5 100644 --- a/src/Interpolations/mpvalue.jl +++ b/src/Interpolations/mpvalue.jl @@ -103,14 +103,14 @@ end (value(hessian, f, x, args...)..., gradient(∇∇f, x)) end -@inline @propagate_inbounds set_shape_values!(mp::MPValue, ip, (w,)::Tuple{Any}) = (mp.w[ip]=w;) -@inline @propagate_inbounds set_shape_values!(mp::MPValue, ip, (w,∇w)::Tuple{Any,Any}) = (mp.w[ip]=w; mp.∇w[ip]=∇w;) -@inline @propagate_inbounds set_shape_values!(mp::MPValue, ip, (w,∇w,∇∇w)::Tuple{Any,Any,Any}) = (mp.w[ip]=w; mp.∇w[ip]=∇w; mp.∇∇w[ip]=∇∇w;) -@inline @propagate_inbounds set_shape_values!(mp::MPValue, ip, (w,∇w,∇∇w,∇∇∇w)::Tuple{Any,Any,Any,Any}) = (mp.w[ip]=w; mp.∇w[ip]=∇w; mp.∇∇w[ip]=∇∇w; mp.∇∇∇w[ip]=∇∇∇w;) -@inline set_shape_values!(mp::MPValue, (w,)::Tuple{Any}) = copyto!(mp.w, w) -@inline set_shape_values!(mp::MPValue, (w,∇w)::Tuple{Any,Any}) = (copyto!(mp.w,w); copyto!(mp.∇w,∇w);) -@inline set_shape_values!(mp::MPValue, (w,∇w,∇∇w)::Tuple{Any,Any,Any}) = (copyto!(mp.w,w); copyto!(mp.∇w,∇w); copyto!(mp.∇∇w,∇∇w);) -@inline set_shape_values!(mp::MPValue, (w,∇w,∇∇w,∇∇∇w)::Tuple{Any,Any,Any,Any}) = (copyto!(mp.w,w); copyto!(mp.∇w,∇w); copyto!(mp.∇∇w,∇∇w); copyto!(mp.∇∇∇w,∇∇∇w);) +@inline @propagate_inbounds set_kernel_values!(mp::MPValue, ip, (w,)::Tuple{Any}) = (mp.w[ip]=w;) +@inline @propagate_inbounds set_kernel_values!(mp::MPValue, ip, (w,∇w)::Tuple{Any,Any}) = (mp.w[ip]=w; mp.∇w[ip]=∇w;) +@inline @propagate_inbounds set_kernel_values!(mp::MPValue, ip, (w,∇w,∇∇w)::Tuple{Any,Any,Any}) = (mp.w[ip]=w; mp.∇w[ip]=∇w; mp.∇∇w[ip]=∇∇w;) +@inline @propagate_inbounds set_kernel_values!(mp::MPValue, ip, (w,∇w,∇∇w,∇∇∇w)::Tuple{Any,Any,Any,Any}) = (mp.w[ip]=w; mp.∇w[ip]=∇w; mp.∇∇w[ip]=∇∇w; mp.∇∇∇w[ip]=∇∇∇w;) +@inline set_kernel_values!(mp::MPValue, (w,)::Tuple{Any}) = copyto!(mp.w, w) +@inline set_kernel_values!(mp::MPValue, (w,∇w)::Tuple{Any,Any}) = (copyto!(mp.w,w); copyto!(mp.∇w,∇w);) +@inline set_kernel_values!(mp::MPValue, (w,∇w,∇∇w)::Tuple{Any,Any,Any}) = (copyto!(mp.w,w); copyto!(mp.∇w,∇w); copyto!(mp.∇∇w,∇∇w);) +@inline set_kernel_values!(mp::MPValue, (w,∇w,∇∇w,∇∇∇w)::Tuple{Any,Any,Any,Any}) = (copyto!(mp.w,w); copyto!(mp.∇w,∇w); copyto!(mp.∇∇w,∇∇w); copyto!(mp.∇∇∇w,∇∇∇w);) function Base.show(io::IO, mp::MPValue) print(io, "MPValue: \n") diff --git a/src/Tesserae.jl b/src/Tesserae.jl index 5cd687da..efb0d5d2 100644 --- a/src/Tesserae.jl +++ b/src/Tesserae.jl @@ -53,6 +53,10 @@ export LinearBSpline, QuadraticBSpline, CubicBSpline, + SteffenBSpline, + SteffenLinearBSpline, + SteffenQuadraticBSpline, + SteffenCubicBSpline, GIMP, Interpolation, KernelCorrection, diff --git a/test/interpolations.jl b/test/interpolations.jl index da15ccfb..4f8ab09c 100644 --- a/test/interpolations.jl +++ b/test/interpolations.jl @@ -3,7 +3,7 @@ @test eltype(MPValue(Vec{dim}, QuadraticBSpline())) == eltype(MPValue(Vec{dim, Float64}, QuadraticBSpline())) for T in (Float32, Float64) - for kernel in (LinearBSpline(), QuadraticBSpline(), CubicBSpline(), GIMP()) + for kernel in (LinearBSpline(), QuadraticBSpline(), CubicBSpline(), SteffenLinearBSpline(), SteffenQuadraticBSpline(), SteffenCubicBSpline(), GIMP()) for extension in (identity, KernelCorrection) it = extension(kernel) mp = @inferred MPValue(Vec{dim,T}, it) @@ -59,25 +59,42 @@ end @testset "Interpolations" begin + isapproxzero(x) = x + ones(x) ≈ ones(x) function check_partition_of_unity(mp, x; atol=sqrt(eps(eltype(mp.w)))) indices = neighboringnodes(mp) CI = CartesianIndices(indices) # local indices - isapprox(sum(mp.w[CI]), 1) && isapprox(sum(mp.∇w[CI]), zero(eltype(mp.∇w)); atol) + isapprox(sum(mp.w[CI]), 1) && isapproxzero(sum(mp.∇w[CI])) end - function check_linear_field_reproduction(mp, x, X; atol=sqrt(eps(eltype(mp.w)))) + function check_linear_field_reproduction(mp, x, X) indices = neighboringnodes(mp) CI = CartesianIndices(indices) # local indices - isapprox(mapreduce((j,i) -> X[i]*mp.w[j], +, CI, indices), x; atol) && - isapprox(mapreduce((j,i) -> X[i]⊗mp.∇w[j], +, CI, indices), I; atol) + isapprox(mapreduce((j,i) -> X[i]*mp.w[j], +, CI, indices), x) && + isapprox(mapreduce((j,i) -> X[i]⊗mp.∇w[j], +, CI, indices), I) end @testset "$it" for it in (LinearBSpline(), QuadraticBSpline(), CubicBSpline()) for T in (Float32, Float64), dim in (1,2,3) Random.seed!(1234) - mp = MPValue(Vec{dim,T}, it) + mp = MPValue(Vec{dim}, it) mesh = CartesianMesh(0.1, ntuple(i->(0,1), Val(dim))...) @test all(1:100) do _ - x = rand(Vec{dim, T}) + x = rand(Vec{dim}) + update!(mp, x, mesh) + isnearbounds = size(mp.w) != size(neighboringnodes(mp)) + PU = check_partition_of_unity(mp, x) + LFR = check_linear_field_reproduction(mp, x, mesh) + isnearbounds ? (!PU && !LFR) : (PU && LFR) + end + end + end + + @testset "$it" for it in (SteffenLinearBSpline(), SteffenQuadraticBSpline(), SteffenCubicBSpline()) + for dim in (1,2,3) + Random.seed!(1234) + mp = MPValue(Vec{dim}, it) + mesh = CartesianMesh(0.1, ntuple(i->(0,1), Val(dim))...) + @test all(1:100) do _ + x = rand(Vec{dim}) update!(mp, x, mesh) isnearbounds = size(mp.w) != size(neighboringnodes(mp)) PU = check_partition_of_unity(mp, x) @@ -89,13 +106,13 @@ end @testset "GIMP()" begin it = GIMP() - for T in (Float32, Float64), dim in (1,2,3) + for dim in (1,2,3) Random.seed!(1234) - mp = MPValue(Vec{dim,T}, it) + mp = MPValue(Vec{dim}, it) mesh = CartesianMesh(0.1, ntuple(i->(0,1), Val(dim))...) l = 0.5*spacing(mesh) / 2 @test all(1:100) do _ - x = rand(Vec{dim, T}) + x = rand(Vec{dim}) update!(mp, (;x,l), mesh) isnearbounds = any(.!(l .< x .< 1-l)) PU = check_partition_of_unity(mp, x) @@ -107,15 +124,15 @@ end end end - @testset "KernelCorrection($kernel)" for kernel in (LinearBSpline(), QuadraticBSpline(), CubicBSpline(), GIMP()) + @testset "KernelCorrection($kernel)" for kernel in (LinearBSpline(), QuadraticBSpline(), CubicBSpline(), SteffenLinearBSpline(), SteffenQuadraticBSpline(), SteffenCubicBSpline(), GIMP()) it = KernelCorrection(kernel) - for T in (Float32, Float64), dim in (1,2,3) + for dim in (1,2,3) Random.seed!(1234) - mp = MPValue(Vec{dim,T}, it) + mp = MPValue(Vec{dim}, it) mesh = CartesianMesh(0.1, ntuple(i->(0,1), Val(dim))...) l = 0.5*spacing(mesh) / 2 @test all(1:100) do _ - x = rand(Vec{dim, T}) + x = rand(Vec{dim}) update!(mp, (;x,l), mesh) PU = check_partition_of_unity(mp, x) LFR = check_linear_field_reproduction(mp, x, mesh)