piecewise multiplication (#10)

* piecewise multiplication

* Update runtests.jl

examples/heat.jl

@@ -1,27 +1,55 @@
-using PiecewiseOrthogonalPolynomials, ClassicalOrthogonalPolynomials, DifferentialEquations, BlockArrays, Plots
+using PiecewiseOrthogonalPolynomials, ClassicalOrthogonalPolynomials, DifferentialEquations, BlockArrays, BlockBandedMatrices, Plots
+
+heat(u, (Δ, M), t) = M \ (Δ * u)
 
 r = range(-1,1; length=4)
 C = ContinuousPolynomial{1}(r)
 P = ContinuousPolynomial{0}(r)
 x = axes(C,1)
 D = Derivative(x)
 
-K = Block.(1:20)
+K = Block.(1:40)
 M = (C'C)[K,K]
 Δ = BlockBandedMatrix((-(D*C)'*(D*C))[K,K])
 
 
 c₀ = M \ ((C'P)*(P \ exp.(-30x.^2)))[K]
 
-heat(u, (Δ, M), t) = M \ (Δ * u)
-
 u = solve(ODEProblem(heat, c₀, (0.,1.), (Δ, M)), Tsit5(), reltol=1e-8, abstol=1e-8)
 
 g = range(-1,1;length=100)
 plot(g, C[g,K]*u(0.0))
 plot!(g, C[g,K]*u(0.5))
 plot!(g, C[g,K]*u(1))
 
-s = P / P \ broadcast(x -> abs(x) ≤ 1/3 ? 1.0 : 0.5, x)
+# u(t,x) = C[x,:] * c(t)
+
+# a(x) * u_tt = Δ*u
+# C' * (a. * C) * c_tt = -((D*C') * D*C) * c
+
+K = Block.(1:20)
+a = P / P \ broadcast(x -> abs(x) ≤ 1/3 ? 0.5 : 1.0, x)
+M̃ = (C'*(a .* C))[K,K]
+M = (C'C)[K,K]
+Δ = BlockBandedMatrix((-(D*C)'*(D*C))[K,K])
 
-s .* P
+c₀ = M \ ((C'P)*(P \ exp.(-30x.^2)))[K]
+
+u = solve(ODEProblem(heat, c₀, (0.,1.), (Δ, M̃)), Tsit5(), reltol=1e-8, abstol=1e-8)
+
+g = range(-1,1;length=100)
+p = plot()
+for t in range(0,1; length=4)
+    plot!(g, C[g,K]*u(t))
+end; p
+
+
+wave(up, u, (Δ, M), t) = M \ (Δ * u)
+
+u = solve(SecondOrderODEProblem(wave, zero(c₀), c₀, (0.,1.), (Δ, M̃)), Tsit5(), reltol=1e-8, abstol=1e-8)
+
+g = range(-1,1;length=100)
+p = plot()
+for t in range(0,1; length=10)
+    plot!(g, C[g,K]*u(t).x[2])
+end; p

src/PiecewiseOrthogonalPolynomials.jl

@@ -4,13 +4,18 @@ using ClassicalOrthogonalPolynomials, LinearAlgebra, BlockArrays, BlockBandedMat
 import BlockArrays: BlockSlice, block, blockindex, blockvec
 import BlockBandedMatrices: _BandedBlockBandedMatrix
 import ClassicalOrthogonalPolynomials: grid, massmatrix
-import ContinuumArrays: @simplify, factorize, TransformFactorization
+import ContinuumArrays: @simplify, factorize, TransformFactorization, AbstractBasisLayout, MemoryLayout, layout_broadcasted, ExpansionLayout, basis
+import LazyArrays: paddeddata
+import LazyBandedMatrices: BlockBroadcastMatrix
 import Base: axes, getindex, ==, \, OneTo
 
 export PiecewisePolynomial, ContinuousPolynomial, Derivative, Block
 
 abstract type AbstractPiecewisePolynomial{order,T,P<:AbstractVector} <: Basis{T} end
 
+struct PiecewisePolynomialLayout{order} <: AbstractBasisLayout end
+MemoryLayout(::Type{<:AbstractPiecewisePolynomial{order}}) where order = PiecewisePolynomialLayout{order}()
+
 struct PiecewisePolynomial{T,Bas,P<:AbstractVector} <: AbstractPiecewisePolynomial{0,T,P}
     basis::Bas
     points::P
@@ -179,5 +184,27 @@ end
 end
 
 
+### multiplication
+
+function layout_broadcasted(::Tuple{ExpansionLayout{PiecewisePolynomialLayout{0}},PiecewisePolynomialLayout{0}}, ::typeof(*), a, P)
+    @assert basis(a) == P
+    _,c = a.args
+    T = eltype(c)
+    m = length(P.points)
+    B = PiecewisePolynomial(P).basis
+
+    ops = [(cₖ = [paddeddata(c)[k:m-1:end]; Zeros{T}(∞)];
+            aₖ = B  * cₖ;
+            unitblocks(B \ (aₖ .* B))) for k = 1:m-1]
+
+    P * BlockBroadcastMatrix{T}(Diagonal, ops...)
+end
+
+
+function layout_broadcasted(::Tuple{ExpansionLayout{PiecewisePolynomialLayout{0}},PiecewisePolynomialLayout{1}}, ::typeof(*), a, C)
+    P = ContinuousPolynomial{0}(C)
+    (a .* P) * (P \ C)
+end
+
 
 end # module

test/runtests.jl

@@ -1,4 +1,4 @@
-using PiecewiseOrthogonalPolynomials, ClassicalOrthogonalPolynomials, BlockArrays, Test
+using PiecewiseOrthogonalPolynomials, ClassicalOrthogonalPolynomials, BlockArrays, Test, FillArrays, LinearAlgebra
 
 
 @testset "transform" begin
@@ -56,6 +56,19 @@ end
     P\D*C
 end
 
+@testset "multiplication" begin
+    r = range(-1,1; length=4)
+    P = ContinuousPolynomial{0}(r)
+    C = ContinuousPolynomial{1}(r)
+    x = axes(P,1)
+    a = P / P \ broadcast(x -> abs(x) ≤ 1/3 ? 1.0 : 0.5, x)
+    @test (P \ (a .* P))[Block.(1:2), Block.(1:2)] ≈ Diagonal([0.5,1,0.5,0.5,1,0.5])
+
+    c = [randn(4); Zeros(∞)]
+    @test ((a .* C) * c)[0.1] ≈ (C*c)[0.1]
+    @test ((a .* C) * c)[0.4] ≈ (C*c)[0.4]/2
+end
+
 @testset "weak Laplacian" begin
     r = range(0,1; length=5)
     C = ContinuousPolynomial{1}(r)