Grammatrix fore Triangles, Makie experiments (#160)

* Update Makie code

* Update MultivariateOrthogonalPolynomialsMakieExt.jl

* tricontourf

* Add Triangle

* plotvalues no grid

* Update Project.toml

* Add grammatrix and Conjugate matrix

* Update Project.toml

* Update Project.toml

* B -> Q in Conjugate

* Update test_triangle.jl

* Update Project.toml

* Update Project.toml

* Increase coverage

* increase coverage

* Update test_triangle.jl

Author: dlfivefifty

Project.toml

@@ -1,6 +1,6 @@
 name = "MultivariateOrthogonalPolynomials"
 uuid = "4f6956fd-4f93-5457-9149-7bfc4b2ce06d"
-version = "0.5.1"
+version = "0.6"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -27,18 +27,18 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 ArrayLayouts = "1.0.9"
 BandedMatrices = "0.17.30"
 BlockArrays = "0.16.14"
-BlockBandedMatrices = "0.12"
-ClassicalOrthogonalPolynomials = "0.10, 0.11"
-ContinuumArrays = "0.13, 0.14"
+BlockBandedMatrices = "0.12.5"
+ClassicalOrthogonalPolynomials = "0.11"
+ContinuumArrays = "0.15"
 DomainSets = "0.6"
 FastTransforms = "0.15"
 FillArrays = "1.0"
-HarmonicOrthogonalPolynomials = "0.4"
+HarmonicOrthogonalPolynomials = "0.5"
 InfiniteArrays = "0.12, 0.13"
-InfiniteLinearAlgebra = "0.6.6"
+InfiniteLinearAlgebra = "0.7"
 LazyArrays = "1.0"
-LazyBandedMatrices = "0.8.4"
-QuasiArrays = "0.10, 0.11"
+LazyBandedMatrices = "0.9"
+QuasiArrays = "0.11"
 SpecialFunctions = "1, 2"
 StaticArrays = "1"
 julia = "1.9"

examples/makieplotting.jl

@@ -0,0 +1,282 @@
+# module MultivariateOrthogonalPolynomialsMakieExt
+
+using GLMakie
+using MultivariateOrthogonalPolynomials, ClassicalOrthogonalPolynomials, StaticArrays
+import Makie: mesh, mesh!
+using ContinuumArrays: plotgridvalues
+
+export contourf, contourf
+
+contourf(f::Fun; kwds...) = _mesh(meshdata(f)...; shading=false, kwds...)
+contourf!(s, f::Fun; kwds...) = _mesh!(s,  meshdata(f)...; shading=false, kwds...)
+
+
+function _mesh(p, T, v; resolution=(400,400), kwds...)
+    T_mat = Array{Int}(undef, length(T), 3)
+    for k = 1:length(T)
+        T_mat[k,:] .= T[k]
+    end
+    s = Scene(resolution=resolution)
+    mesh!(s, [first.(p) last.(p)], T_mat; color=v, kwds...)
+end
+
+
+function _surface(p, T, v; resolution=(400,400), kwds...)
+    T_mat = Array{Int}(undef, length(T), 3)
+    for k = 1:length(T)
+        T_mat[k,:] .= T[k]
+    end
+    # s = Scene(resolution=resolution)
+    mesh(first.(p), last.(p), vec(v), T_mat; kwds...)
+end
+
+
+
+function _mesh!(s, p, T, v; kwds...)
+    T_mat = Array{Int}(undef, length(T), 3)
+    for k = 1:length(T)
+        T_mat[k,:] .= T[k]
+    end
+    mesh!(s, [first.(p) last.(p)], T_mat; color=v, kwds...)
+end
+
+function meshdata(f::Fun{<:PiecewiseSpace})
+    pTv = MultivariateTriangle.meshdata.(components(f))
+    p = vcat(first.(pTv)...)
+    T = pTv[1][2]
+    cs = length(pTv[1][1])
+    for k = 2:length(pTv)
+        append!(T, (t -> (cs.+t)).(pTv[k][2]))
+        cs += length(pTv[k][1])
+    end
+
+    v = vcat(last.(pTv)...)
+
+    p, T, v
+end
+
+function meshdata(f::Fun{<:TensorSpace{<:Tuple{<:Chebyshev,<:Chebyshev}}})
+    p = points(f)
+    v = values(f)
+    n = length(p)
+    T = Vector{NTuple{3,Int}}()
+    d_x,d_y = factors(domain(f))
+    a_x,b_x = endpoints(d_x)
+    a_y,b_y = endpoints(d_y)
+    if iseven(_padua_length(n))
+        l = floor(Int, (1+sqrt(1+8n))/4)
+
+        push!(p, Vec(b_x,b_y))
+        push!(p, Vec(a_x,b_y))
+
+        push!(v, f(b_x,b_y))
+        push!(v, f(a_x,b_y))
+
+        for p = 0:l-2
+            for k = (2p*l)+1:(2p*l)+l-1
+                push!(T, (k+1, k, l+k+1))
+            end
+            for k = (2p*l)+1:(2p*l)+l-1
+                push!(T, (k, l+k, l+k+1))
+            end
+            for k = (2p*l)+l+1:(2p*l)+2l-1
+                push!(T, (k+1, k, l+k))
+            end
+            for k = (2p*l)+l+2:(2p*l)+2l
+                push!(T, (k, k+l-1, l+k))
+            end
+        end
+        for p=0:l-3
+            push!(T, ((2p+1)*l+1, (2p+2)*l+1, (2p+3)*l+1))
+        end
+        for p =0:l-2
+            push!(T, ((2p+1)*l, (2p+2)*l, (2p+3)*l))
+        end
+        push!(T, (1, n+1, l+1))
+        push!(T, (n-2l+1, n+2, n-l+1))
+    else
+        l = floor(Int, (3+sqrt(1+8n))/4)
+
+        push!(p, Vec(a_x,b_y))
+        push!(p, Vec(a_x,a_y))
+
+        push!(v, f(a_x,b_y))
+        push!(v, f(a_x,a_y))
+
+        for p = 0:l-2
+            for k = p*(2l-1)+1:p*(2l-1)+l-1
+                push!(T, (k+1, k, l+k))
+            end
+            for k = p*(2l-1)+1:p*(2l-1)+l-2
+                push!(T, (k+1, l+k, l+k+1))
+            end
+        end
+        for p = 0:l-3
+            for k = p*(2l-1)+l+1:p*(2l-1)+2l-2
+                push!(T, (k+1, k, l+k))
+            end
+            for k = p*(2l-1)+l+1:p*(2l-1)+2l-1
+                push!(T, (k, k+l-1, l+k))
+            end
+        end
+
+        for p=0:l-3
+            push!(T, (p*(2l-1) + 1, p*(2l-1) + l+1, p*(2l-1) + 2l))
+        end
+
+        for p=0:l-3
+            push!(T, (p*(2l-1) + l, p*(2l-1) + 2l-1, p*(2l-1) + 3l-1))
+        end
+
+        push!(T, (n-2l+2, n+1, n-l+2))
+        push!(T, (n-l+1, n+2, n))
+    end
+
+    p, T, v
+end
+
+
+
+meshdata(f) =
+    triangle_meshdata(points(f), values(f), (domain(f).a, domain(f).b, domain(f).c),
+                                            f.((domain(f).a, domain(f).b, domain(f).c)))
+
+function triangle_meshdata(p, v, (a, b, c), (fa, fb, fc))
+    n = length(p)
+    T = Vector{NTuple{3,Int}}()
+
+
+    if iseven(_padua_length(n))
+        l = floor(Int, (1+sqrt(1+8n))/4)
+
+        push!(p, b)
+        push!(p, c)
+
+        push!(v, fb)
+        push!(v, fc)
+
+        for p = 0:l-2
+            for k = (2p*l)+1:(2p*l)+l-1
+                push!(T, (k+1, k, l+k+1))
+            end
+            for k = (2p*l)+1:(2p*l)+l-1
+                push!(T, (k, l+k, l+k+1))
+            end
+            for k = (2p*l)+l+1:(2p*l)+2l-1
+                push!(T, (k+1, k, l+k))
+            end
+            for k = (2p*l)+l+2:(2p*l)+2l
+                push!(T, (k, k+l-1, l+k))
+            end
+        end
+        for p=0:l-3
+            push!(T, ((2p+1)*l+1, (2p+2)*l+1, (2p+3)*l+1))
+        end
+        for p =0:l-2
+            push!(T, ((2p+1)*l, (2p+2)*l, (2p+3)*l))
+        end
+        push!(T, (1, n+1, l+1))
+        push!(T, (n-2l+1, n+2, n-l+1))
+    else
+        l = floor(Int, (3+sqrt(1+8n))/4)
+
+        push!(p, Vec(c))
+        push!(p, Vec(a))
+
+        push!(v, fc)
+        push!(v, fa)
+
+        for p = 0:l-2
+            for k = p*(2l-1)+1:p*(2l-1)+l-1
+                push!(T, (k+1, k, l+k))
+            end
+            for k = p*(2l-1)+1:p*(2l-1)+l-2
+                push!(T, (k+1, l+k, l+k+1))
+            end
+        end
+        for p = 0:l-3
+            for k = p*(2l-1)+l+1:p*(2l-1)+2l-2
+                push!(T, (k+1, k, l+k))
+            end
+            for k = p*(2l-1)+l+1:p*(2l-1)+2l-1
+                push!(T, (k, k+l-1, l+k))
+            end
+        end
+
+        for p=0:l-3
+            push!(T, (p*(2l-1) + 1, p*(2l-1) + l+1, p*(2l-1) + 2l))
+        end
+
+        for p=0:l-3
+            push!(T, (p*(2l-1) + l, p*(2l-1) + 2l-1, p*(2l-1) + 3l-1))
+        end
+
+        push!(T, (n-2l+2, n+1, n-l+2))
+        push!(T, (n-l+1, n+2, n))
+    end
+
+    p, T, v
+end
+
+
+P = JacobiTriangle()
+f = expand(P, splat((x,y) -> cos(x*exp(y))))
+(a,b,c) = (SVector(0.,0.), SVector(0.,1.), SVector(1.,0.))
+xy,F = plotgridvalues(f)
+x,y = first.(xy),last.(xy)
+
+
+
+triangle_meshdata(plotgridvalues(f)..., (a,b,c), getindex.(Ref(f), (a,b,c)))
+
+
+tricontourf(vec(x), vec(y), vec(F); levels=100)
+
+using Makie, MultivariateOrthogonalPolynomials
+using ContinuumArrays: ApplyQuasiVector, plotgridvalues
+# Makie.plottype(a::ApplyQuasiVector{<:Any, typeof(*), <:Tuple{JacobiTriangle,AbstractVector}}) = Tricontourf
+
+
+function Makie.tricontourf(f::ApplyQuasiVector{<:Any, typeof(*)}; kwds...)
+    xy,F = plotgridvalues(f)
+    x,y = first.(xy),last.(xy)
+    tricontourf(vec(x), vec(y), vec(F); kwds...)
+end
+
+function Makie.tricontourf!(f::ApplyQuasiVector{<:Any, typeof(*)}; kwds...)
+    xy,F = plotgridvalues(f)
+    x,y = first.(xy),last.(xy)
+    tricontourf!(vec(x), vec(y), vec(F); kwds...)
+end
+
+
+# f = expand(Weighted(JacobiTriangle(1,1,1))[:,20])
+# tricontourf(f; levels=100)
+import Base: oneto
+KR = Block.(oneto(40))
+P = JacobiTriangle()
+f = P[:,KR] * (P \ Weighted(JacobiTriangle(1,1,1)))[KR,10]
+
+p = Figure(resolution = (3200, 2400))
+ax = Axis(p[1, 1])
+tricontourf!(f; levels=100)
+save("bubble10.png", p)
+
+
+using ContinuumArrays: AffineMap, affine
+using MultivariateOrthogonalPolynomials
+using MultivariateOrthogonalPolynomials: Triangle
+using StaticArrays
+
+a = affine(Triangle(), Triangle(SVector(1,0), SVector(0,1), SVector(1,1)))
+a[SVector(0.1,0.2)]
+
+P = JacobiTriangle()
+Q = P[affine(Triangle(SVector(1,0), SVector(0,1), SVector(1,1)), axes(P,1)), :]
+
+f = Q * [[randn(20); zeros(100)]; zeros(∞)]
+
+tricontourf(f; nlevels=300)
+
+
+# end # modul

src/MultivariateOrthogonalPolynomials.jl

@@ -6,19 +6,19 @@ using ClassicalOrthogonalPolynomials, FastTransforms, BlockBandedMatrices, Block
       LazyArrays, SpecialFunctions, LinearAlgebra, BandedMatrices, LazyBandedMatrices, ArrayLayouts,
       HarmonicOrthogonalPolynomials
 
-import Base: axes, in, ==, *, ^, \, copy, copyto!, OneTo, getindex, size, oneto, all, resize!, BroadcastStyle, similar, fill!, setindex!, convert, show, summary
+import Base: axes, in, ==, +, -, /, *, ^, \, copy, copyto!, OneTo, getindex, size, oneto, all, resize!, BroadcastStyle, similar, fill!, setindex!, convert, show, summary
 import Base.Broadcast: Broadcasted, broadcasted, DefaultArrayStyle
-import DomainSets: boundary
+import DomainSets: boundary, EuclideanDomain
 
 import QuasiArrays: LazyQuasiMatrix, LazyQuasiArrayStyle, domain
-import ContinuumArrays: @simplify, Weight, weight, grid, plotgrid, TransformFactorization, ExpansionLayout, plotvalues, unweighted, plan_grid_transform, checkpoints, transform_ldiv, AbstractBasisLayout, basis_axes
+import ContinuumArrays: @simplify, Weight, weight, grid, plotgrid, TransformFactorization, ExpansionLayout, plotvalues, unweighted, plan_grid_transform, checkpoints, transform_ldiv, AbstractBasisLayout, basis_axes, Inclusion, grammatrix
 
 import ArrayLayouts: MemoryLayout, sublayout, sub_materialize
 import BlockArrays: block, blockindex, BlockSlice, viewblock, blockcolsupport, AbstractBlockStyle, BlockStyle
 import BlockBandedMatrices: _BandedBlockBandedMatrix, AbstractBandedBlockBandedMatrix, _BandedMatrix, blockbandwidths, subblockbandwidths
 import LinearAlgebra: factorize
-import LazyArrays: arguments, paddeddata, LazyArrayStyle, LazyLayout, PaddedLayout, applylayout
-import LazyBandedMatrices: LazyBandedBlockBandedLayout, AbstractBandedBlockBandedLayout, AbstractLazyBandedBlockBandedLayout, _krontrav_axes, DiagTravLayout
+import LazyArrays: arguments, paddeddata, LazyArrayStyle, LazyLayout, PaddedLayout, applylayout, LazyMatrix, ApplyMatrix
+import LazyBandedMatrices: LazyBandedBlockBandedLayout, AbstractBandedBlockBandedLayout, AbstractLazyBandedBlockBandedLayout, _krontrav_axes, DiagTravLayout, invdiagtrav, ApplyBandedBlockBandedLayout
 import InfiniteArrays: InfiniteCardinal, OneToInf
 
 import ClassicalOrthogonalPolynomials: jacobimatrix, Weighted, orthogonalityweight, HalfWeighted, WeightedBasis, pad, recurrencecoefficients, clenshaw
@@ -32,7 +32,8 @@ export MultivariateOrthogonalPolynomial, BivariateOrthogonalPolynomial,
        DunklXuDisk, DunklXuDiskWeight, WeightedDunklXuDisk,
        Zernike, ZernikeWeight, zerniker, zernikez,
        PartialDerivative, Laplacian, AbsLaplacianPower, AngularMomentum,
-       RadialCoordinate, Weighted, Block, jacobimatrix, KronPolynomial, RectPolynomial
+       RadialCoordinate, Weighted, Block, jacobimatrix, KronPolynomial, RectPolynomial,
+       grammatrix, oneto
 
 include("ModalInterlace.jl")
 include("rect.jl")