Skip to content

Commit

Permalink
Add SteffenBSpline
Browse files Browse the repository at this point in the history
  • Loading branch information
@KeitaNakamura
KeitaNakamura committed Jul 27, 2024
1 parent 22d77f5 commit 4d7c140
Show file tree
Hide file tree
Showing 6 changed files with 207 additions and 157 deletions.
283 changes: 156 additions & 127 deletions src/Interpolations/bspline.jl
View file
Original file line number Diff line number Diff line change
@@ -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)
Expand All @@ -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-ξ, ξ)))
Expand All @@ -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))
Expand All @@ -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))
Expand All @@ -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
Expand Down Expand Up @@ -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)))
2 changes: 1 addition & 1 deletion src/Interpolations/gimp.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Loading

0 comments on commit 4d7c140

Please sign in to comment.