Skip to content

Commit

Permalink
Add quadratic polynomial
Browse files Browse the repository at this point in the history
  • Loading branch information
@KeitaNakamura
KeitaNakamura committed Oct 20, 2024
1 parent 2da49de commit 438467d
Show file tree
Hide file tree
Showing 3 changed files with 33 additions and 8 deletions.
19 changes: 18 additions & 1 deletion src/Interpolations/polynomials.jl
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
struct MultiLinear end
struct MultiQuadratic end

struct Polynomial{D <: Union{Degree, MultiLinear}}
struct Polynomial{D <: Union{Degree, MultiLinear, MultiQuadratic}}
degree::D
end
Base.show(io::IO, poly::Polynomial) = print(io, Polynomial, "(", poly.degree, ")")
Expand All @@ -13,9 +14,14 @@ Base.show(io::IO, poly::Polynomial) = print(io, Polynomial, "(", poly.degree, ")
end
end

@inline Base.values(::Order{k}, p::Polynomial{<: Union{Quadratic, MultiQuadratic}}, x::Vec) where {k} = ∂ⁿ{k,:all}(x->value(p,x), x)

@inline _value(::Order{0}, ::Polynomial{Linear}, x::Vec) = vcat(one(eltype(x)), x)
@inline _value(::Order{1}, ::Polynomial{Linear}, x::Vec{dim, T}) where {dim, T} = vcat(zero(Mat{1, dim, T}), one(Mat{dim, dim, T}))
@inline _value(::Order{k}, ::Polynomial{Linear}, x::Vec{dim, T}) where {k, dim, T} = zero(Tensor{Tuple{dim+1, @Symmetry{nfill(dim,Val(k))...}}, T})
@inline _value(::Order{0}, ::Polynomial{Quadratic}, x::Vec{1}) = Vec(one(eltype(x)), x[1], x[1]^2)
@inline _value(::Order{0}, ::Polynomial{Quadratic}, x::Vec{2}) = Vec(one(eltype(x)), x[1], x[2], x[1]*x[2], x[1]^2, x[2]^2)
@inline _value(::Order{0}, ::Polynomial{Quadratic}, x::Vec{3}) = Vec(one(eltype(x)), x[1], x[2], x[3], x[1]*x[2], x[2]*x[3], x[3]*x[1], x[1]^2, x[2]^2, x[3]^2)

@inline _value(::Order{k}, ::Polynomial{MultiLinear}, x::Vec{1}) where {k} = _value(Order(k), Polynomial(Linear()), x)
@inline _value(::Order{0}, ::Polynomial{MultiLinear}, x::Vec{2}) = Vec(one(eltype(x)), x[1], x[2], x[1]*x[2])
Expand All @@ -27,3 +33,14 @@ end
@inline _value(::Order{2}, ::Polynomial{MultiLinear}, x::Vec{3, T}) where {T} = Tensor{Tuple{8, @Symmetry{3, 3}}, T}(0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,x[3],0,0,0,0,0,0,1,x[2],0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,x[1],0,0,0,0,0,0,0,0)
@inline _value(::Order{3}, ::Polynomial{MultiLinear}, x::Vec{3, T}) where {T} = Tensor{Tuple{8, @Symmetry{3, 3, 3}}, T}(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)
@inline _value(::Order{k}, ::Polynomial{MultiLinear}, x::Vec{3, T}) where {k, T} = zero(Tensor{Tuple{8, @Symmetry{nfill(3,Val(k))...}}, T})
@inline _value(::Order{0}, ::Polynomial{MultiQuadratic}, x::Vec{1}) = Vec(one(eltype(x)), x[1], x[1]^2)
@inline _value(::Order{0}, ::Polynomial{MultiQuadratic}, x::Vec{2}) = Vec(one(eltype(x)), x[1], x[2], x[1]*x[2], x[1]^2, x[2]^2, x[1]^2*x[2], x[1]*x[2]^2, x[1]^2*x[2]^2)
@inline _value(::Order{0}, ::Polynomial{MultiQuadratic}, x::Vec{3}) = Vec(one(eltype(x)),
x[1], x[2], x[3],
x[1]*x[2], x[2]*x[3], x[3]*x[1], x[1]*x[2]*x[3],
x[1]^2, x[2]^2, x[3]^2,
x[1]^2*x[2], x[1]^2*x[3], x[1]^2*x[2]*x[3],
x[2]^2*x[1], x[2]^2*x[3], x[1]*x[2]^2*x[3],
x[3]^2*x[1], x[3]^2*x[2], x[1]*x[2]*x[3]^2,
x[1]^2*x[2]^2, x[2]^2*x[3]^2, x[3]^2*x[1]^2,
x[1]^2*x[2]^2*x[3], x[1]*x[2]^2*x[3]^2, x[1]^2*x[2]*x[3]^2, x[1]^2*x[2]^2*x[3]^2)
1 change: 1 addition & 0 deletions src/Tesserae.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,6 +55,7 @@ export
Cubic,
Quartic,
MultiLinear,
MultiQuadratic,
Polynomial,
BSpline,
SteffenBSpline,
Expand Down
21 changes: 14 additions & 7 deletions test/interpolations.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -335,14 +335,21 @@ end
@testset "Polynomial" begin
for T in (Float32, Float64)
for dim in (1,2,3)
for poly in (Polynomial(Linear()), Polynomial(MultiLinear()))
for poly in (Polynomial(Linear()), Polynomial(MultiLinear()), Polynomial(Quadratic()), Polynomial(MultiQuadratic()))
xp = rand(Vec{dim, T})
vals = values(Order(4), poly, xp)
@test all(v -> eltype(v) == T, vals)
@test vals[2] ≈ gradient(x -> only(values(Order(0), poly, x)), xp)
@test vals[3] ≈ gradient(x -> gradient(x -> only(values(Order(0), poly, x)), x), xp)
@test vals[4] ≈ gradient(x -> gradient(x -> gradient(x -> only(values(Order(0), poly, x)), x), x), xp)
@test vals[5] ≈ gradient(x -> gradient(x -> gradient(x -> gradient(x -> only(values(Order(0), poly, x)), x), x), x), xp)
if poly.degree isa Union{Linear, MultiLinear}
vals = values(Order(4), poly, xp)
@test all(v -> eltype(v) == T, vals)
@test vals[2] ≈ gradient(x -> only(values(Order(0), poly, x)), xp)
@test vals[3] ≈ gradient(x -> gradient(x -> only(values(Order(0), poly, x)), x), xp)
@test vals[4] ≈ gradient(x -> gradient(x -> gradient(x -> only(values(Order(0), poly, x)), x), x), xp)
@test vals[5] ≈ gradient(x -> gradient(x -> gradient(x -> gradient(x -> only(values(Order(0), poly, x)), x), x), x), xp)
else
vals = values(Order(2), poly, xp)
@test all(v -> eltype(v) == T, vals)
@test vals[2] ≈ gradient(x -> only(values(Order(0), poly, x)), xp)
@test vals[3] ≈ gradient(x -> gradient(x -> only(values(Order(0), poly, x)), x), xp)
end
end
end
end
Expand Down

0 comments on commit 438467d

Please sign in to comment.