Support more triangle gram matrices (#189)

* Support more triangle gram matrices

* Update ci.yml

* Add weightedgrammatrix overload

* Weak formulation tests

* tests should pass

* Update test_triangle.jl

Author: dlfivefifty

.github/workflows/ci.yml

@@ -13,7 +13,6 @@ jobs:
           - '1.10'
         os:
           - ubuntu-latest
-          - macOS-latest
           - windows-latest
         arch:
           - x64

Project.toml

@@ -1,6 +1,6 @@
 name = "MultivariateOrthogonalPolynomials"
 uuid = "4f6956fd-4f93-5457-9149-7bfc4b2ce06d"
-version = "0.7"
+version = "0.7.1"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -37,7 +37,7 @@ HarmonicOrthogonalPolynomials = "0.6"
 InfiniteArrays = "0.14"
 InfiniteLinearAlgebra = "0.8"
 LazyArrays = "2.0"
-LazyBandedMatrices = "0.10"
+LazyBandedMatrices = "0.10.3"
 QuasiArrays = "0.11"
 SpecialFunctions = "1, 2"
 StaticArrays = "1"

src/MultivariateOrthogonalPolynomials.jl

@@ -21,7 +21,7 @@ import LazyArrays: arguments, paddeddata, LazyArrayStyle, LazyLayout, PaddedLayo
 import LazyBandedMatrices: LazyBandedBlockBandedLayout, AbstractBandedBlockBandedLayout, AbstractLazyBandedBlockBandedLayout, _krontrav_axes, DiagTravLayout, invdiagtrav, ApplyBandedBlockBandedLayout
 import InfiniteArrays: InfiniteCardinal, OneToInf
 
-import ClassicalOrthogonalPolynomials: jacobimatrix, Weighted, orthogonalityweight, HalfWeighted, WeightedBasis, pad, recurrencecoefficients, clenshaw
+import ClassicalOrthogonalPolynomials: jacobimatrix, Weighted, orthogonalityweight, HalfWeighted, WeightedBasis, pad, recurrencecoefficients, clenshaw, weightedgrammatrix
 import HarmonicOrthogonalPolynomials: BivariateOrthogonalPolynomial, MultivariateOrthogonalPolynomial, Plan,
                                           PartialDerivative, AngularMomentum, BlockOneTo, BlockRange1, interlace,
                                           MultivariateOPLayout, MAX_PLOT_BLOCKS

src/triangle.jl

@@ -188,15 +188,22 @@ end
 #     _lap_mul(P, eltype(axes(P,1)))
 # end
 
-# _BandedBlockBandedMatrix((@. exp(loggamma(n+k+b+c)+loggamma(n-k+a+1)+loggamma(k+b)+loggamma(k+c)-loggamma(n+k+a+b+c)-loggamma(k+b+c)-loggamma(n-k+1)-loggamma(k))/((2n+a+b+c)*(2k+b+c-1)))',
-# axes(k,1), (0,0), (0,0))
+
 
 function grammatrix(A::JacobiTriangle)
     @assert A == JacobiTriangle()
     n = mortar(Fill.(oneto(∞),oneto(∞)))
     k = mortar(Base.OneTo.(oneto(∞)))
     _BandedBlockBandedMatrix(BroadcastVector{eltype(A)}((n,k) -> exp(loggamma(n+k)+loggamma(n-k+1)+loggamma(k)+loggamma(k)-loggamma(n+k)-loggamma(k)-loggamma(n-k+1)-loggamma(k))/((2n)*(2k-1)), n, k)',
-                                axes(k,1), (0,0), (0,0))
+    axes(k,1), (0,0), (0,0))
+end
+
+function weightedgrammatrix(A::JacobiTriangle)
+    n = mortar(Fill.(oneto(∞),oneto(∞)))
+    a,b,c = A.a,A.b,A.c
+    k = mortar(Base.OneTo.(oneto(∞)))
+    _BandedBlockBandedMatrix(BroadcastVector{eltype(A)}((n,k,a,b,c) -> exp(loggamma(n+k+b+c)+loggamma(n-k+a+1)+loggamma(k+b)+loggamma(k+c)-loggamma(n+k+a+b+c)-loggamma(k+b+c)-loggamma(n-k+1)-loggamma(k))/((2n+a+b+c)*(2k+b+c-1)), n, k, a, b, c)',
+    axes(k,1), (0,0), (0,0))
 end
 
 """
@@ -399,7 +406,7 @@ function \(w_A::WeightedTriangle, B::JacobiTriangle)
     w_A \ (TriangleWeight(0,0,0) .* B)
 end
 
-function \(A::JacobiTriangle, B::JacobiTriangle)
+@simplify function \(A::JacobiTriangle, B::JacobiTriangle)
     if A == B
         Eye{promote_type(eltype(A),eltype(B))}((axes(B,2),))
     elseif A.a == B.a + 1 && A.b == B.b && A.c == B.c
@@ -742,6 +749,23 @@ end
 getindex(f::ApplyQuasiVector{T,typeof(*),<:Tuple{JacobiTriangle,AbstractVector}}, xys::AbstractArray{<:SVector{2}}) where T =
     [f[xy] for xy in xys]
 
+####
+# Grammatrices
+####
+
+"""
+    legendre_triangle_grammatrix
+
+computes the grammatrix by first re-expanding in Legendre
+"""
+function legendre_triangle_grammatrix(A, B)
+    P = JacobiTriangle{eltype(B)}()
+    (P\A)'*grammatrix(P)*(P\B)
+end
+
+@simplify *(Ac::QuasiAdjoint{<:Any,<:JacobiTriangle}, B::Weighted{<:Any,<:JacobiTriangle}) = legendre_triangle_grammatrix(parent(Ac),B)
+
+
 # getindex(f::Expansion{T,<:JacobiTriangle}, x::AbstractVector{<:Number}) where T =
 #     copyto!(Vector{T}(undef, length(x)), view(f, x))
 
@@ -836,4 +860,12 @@ function plotvalues(u::ApplyQuasiVector{T,typeof(*),<:Tuple{JacobiTriangle, Abst
     F = TriIPlan{T}(Block(N), P.a, P.b, P.c)
     C = F * DiagTrav(invdiagtrav(c)[1:N,1:N]) # transform to grid
     C
+end
+
+
+
+function Base._sum(P::JacobiTriangle{T}, dims) where T
+    @assert dims == 1
+    @assert P.a == P.b == P.c == 0
+    Hcat(one(T)/2, Zeros{T}(1,∞))
 end

test/test_triangle.jl

@@ -1,6 +1,6 @@
 using MultivariateOrthogonalPolynomials, StaticArrays, BlockArrays, BlockBandedMatrices, ArrayLayouts,
         QuasiArrays, Test, ClassicalOrthogonalPolynomials, BandedMatrices, FastTransforms, LinearAlgebra, ContinuumArrays
-import MultivariateOrthogonalPolynomials: tri_forwardrecurrence, grid, TriangleRecurrenceA, TriangleRecurrenceB, TriangleRecurrenceC, xy_muladd!, ExpansionLayout, Triangle, ApplyBandedBlockBandedLayout
+import MultivariateOrthogonalPolynomials: tri_forwardrecurrence, grid, TriangleRecurrenceA, TriangleRecurrenceB, TriangleRecurrenceC, xy_muladd!, ExpansionLayout, Triangle, ApplyBandedBlockBandedLayout, weightedgrammatrix
 
 @testset "Triangle" begin
     @testset "basics" begin
@@ -522,14 +522,71 @@ import MultivariateOrthogonalPolynomials: tri_forwardrecurrence, grid, TriangleR
                     ∂ʸ = PartialDerivative{2}(𝐱)
                     Pfx = ∂ˣ * Pf
                     Pfy = ∂ʸ * Pf
-                    for _ in 1:100 
+
+                    @test (coefficients(∂ˣ * P) * coefficients(Pf))[1:50] ≈ coefficients(Pfx)[1:50] ≈ coefficients(∂ˣ * P)[1:50,1:50] * coefficients(Pf)[1:50]
+                    for _ in 1:100
                         u, v = minmax(rand(2)...)
                         x, y = v - u, 1 - v # random point in the triangle
                         @test Pfx[SVector(x, y)] ≈ dfx((x, y))
                         @test Pfy[SVector(x, y)] ≈ dfy((x, y))
                     end
-                end 
+                end
             end
         end
     end
+
+    @testset "Weighted grammatrix" begin
+        P = JacobiTriangle()
+        D = weightedgrammatrix(P)
+
+        KR = Block.(1:10)
+        @test D[KR,KR] == grammatrix(P)[KR,KR]
+        for k = 1:5
+            @test sum(P[:,k].^2) ≈ D[k,k]
+        end
+
+        Q = JacobiTriangle(1,1,1)
+        D = weightedgrammatrix(Q)
+
+
+
+        𝐱 = axes(P,1)
+        x,y = first.(𝐱),last.(𝐱)
+        for k = 1:5
+            @test sum(x .* y .* (1 .- x .- y) .* Q[:,k].^2) ≈ D[k,k]
+        end
+
+
+        W = Weighted(Q)
+
+        PW = P'W
+        WP = W'P
+        for k = 1:5, j=1:5
+            @test sum(W[:,j] .* P[:,k]) ≈ PW[k,j] atol=1E-14
+            @test PW[k,j] ≈ WP[j,k] atol=1E-14
+        end
+
+        f = expand(P, splat((x,y) -> exp(x*cos(y))))
+
+        c = W'f
+
+        for k = 1:5
+            @test c[k] ≈ sum(expand(P, splat((x,y) -> (W[SVector(x,y),k]::Float64) * exp(x*cos(y)))))
+        end
+    end
+
+    @testset "Weak formulation" begin
+        P = JacobiTriangle()
+        W = Weighted(JacobiTriangle(1,1,1))
+        𝐱 = axes(W,1)
+        ∂_x = PartialDerivative{1}(𝐱)
+        ∂_y = PartialDerivative{2}(𝐱)
+        Δ = -((∂_x*W)'*(∂_x*W) + (∂_y*W)'*(∂_y*W))
+        M = W'W
+        f = expand(P, splat((x,y) -> exp(x*cos(y))))
+        κ = 10
+        A = Δ + κ^2*M
+        c = \(A, W'f; tolerance=1E-4)
+        @test (W*c)[SVector(0.1,0.2)] ≈ -0.005927539175184257 # empirical
+    end
 end