Add derivatives for WeightedTriangle (#179)

* Add derivatives for WeightedTriangle

* fix bug, update CI

* FastTransforms v0.16

---------

Co-authored-by: Sheehan Olver <solver@mac.com>

Author: DanielVandH

.github/workflows/ci.yml

@@ -11,19 +11,20 @@ jobs:
       matrix:
         version:
           - '1.9'
+          - '1'
         os:
           - ubuntu-latest
           - macOS-latest
           - windows-latest
         arch:
           - x64
     steps:
-      - uses: actions/checkout@v2
-      - uses: julia-actions/setup-julia@v1
+      - uses: actions/checkout@v4
+      - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
-      - uses: actions/cache@v1
+      - uses: actions/cache@v4
         env:
           cache-name: cache-artifacts
         with:
@@ -36,6 +37,7 @@ jobs:
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
       - uses: julia-actions/julia-processcoverage@v1
-      - uses: codecov/codecov-action@v1
+      - uses: codecov/codecov-action@v4
         with:
+          token: ${{ secrets.CODECOV_TOKEN }}
           file: lcov.info

Project.toml

@@ -31,7 +31,7 @@ BlockBandedMatrices = "0.12.5"
 ClassicalOrthogonalPolynomials = "0.12"
 ContinuumArrays = "0.17"
 DomainSets = "0.6, 0.7"
-FastTransforms = "0.15.11"
+FastTransforms = "0.15.11, 0.16"
 FillArrays = "1.0"
 HarmonicOrthogonalPolynomials = "0.5"
 InfiniteArrays = "0.13"

src/rect.jl

@@ -138,7 +138,7 @@ function plan_transform(P::KronPolynomial{d}, B::Tuple{Block{1}}, dims=1:1) wher
 end
 
 applylayout(::Type{typeof(*)}, ::Lay, ::DiagTravLayout) where Lay <: AbstractBasisLayout = ExpansionLayout{Lay}()
-ContinuumArrays._mul_plotgrid(::Tuple{Any,DiagTravLayout{<:PaddedLayout}}, (P,c)) = plotgrid(P, maximum(blockcolsupport(c)))
+ContinuumArrays._mul_plotgrid(::Tuple{Any,DiagTravLayout{<:PaddedLayout}}, (P,c)::NTuple{2,Any}) = plotgrid(P, maximum(blockcolsupport(c)))
 
 pad(C::DiagTrav, ::BlockedUnitRange{RangeCumsum{Int,OneToInf{Int}}}) = DiagTrav(pad(C.array, ∞, ∞))
 

src/triangle.jl

@@ -165,6 +165,24 @@ end
     JacobiTriangle(a,b+1,c+1) * _BandedBlockBandedMatrix((k .+ convert(T, b+c))', axes(k,1), (-1,1), (-1,1))
 end
 
+# TODO: The derivatives below only hold for a, b, c > 0.
+@simplify function *(Dx::PartialDerivative{1}, P::WeightedTriangle)
+    a, b, c = P.P.a, P.P.b, P.P.c
+    n = mortar(Fill.(oneto(∞),oneto(∞)))
+    k = mortar(Base.OneTo.(oneto(∞)))    
+    scale = -(2k .+ (b + c - 1))
+    coeff1 = (k .+ (c - 1)) .* (n .- k .+ 1) ./ scale 
+    coeff2 = k .* (n .- k .+ a) ./ scale 
+    dat = BlockBroadcastArray(hcat, coeff1, coeff2)
+    WeightedTriangle(a-1, b, c-1) * _BandedBlockBandedMatrix(dat', axes(k, 1), (1, -1), (1, 0))
+end
+
+@simplify function *(Dy::PartialDerivative{2}, P::WeightedTriangle)
+    a, b, c = P.P.a, P.P.b, P.P.c 
+    k = mortar(Base.OneTo.(oneto(∞)))
+    T = promote_type(eltype(Dy), eltype(P)) # avoid bug in convert
+    WeightedTriangle(a, b-1, c-1) * _BandedBlockBandedMatrix(-one(T) * k', axes(k, 1), (1, -1), (1, -1))
+end
 
 # @simplify function *(Δ::Laplacian, P)
 #     _lap_mul(P, eltype(axes(P,1)))

test/test_triangle.jl

@@ -501,4 +501,35 @@ import MultivariateOrthogonalPolynomials: tri_forwardrecurrence, grid, TriangleR
 
         @test affine(d, axes(P,1))[affine(axes(P,1), d)[𝐱]] ≈ 𝐱
     end
+
+    @testset "Weighted derivatives" begin
+        for a in (0.001, 0.1, 0.2, 1.0)
+            for b in (0.001, 0.1, 0.2, 1.0)
+                for c in (0.001, 0.1, 0.2, 1.0)
+                    f = let a = a, b = b, c = c
+                        ((x, y),) -> x^a * y^b * (1 - x - y)^c * (x^2 + 3y^2 + x * y)
+                    end
+                    dfx = let a = a, b = b, c = c
+                        ((x, y),) -> x^a * y^b * (1 - x - y)^c * (2x + y - c * (x^2 + 3y^2 + x * y) / (1 - x - y) + a * (x + y + 3 * y^2 / x))
+                    end
+                    dfy = let a = a, b = b, c = c
+                        ((x, y),) -> x^a * y^b * (1 - x - y)^c * (x + 6y + b * (x + x^2 / y + 3y) - c * (x^2 + 3y^2 + x * y) / (1 - x - y))
+                    end
+                    P = Weighted(JacobiTriangle(a, b, c))
+                    Pf = expand(P, f)
+                    𝐱 = axes(P, 1)
+                    ∂ˣ = PartialDerivative{1}(𝐱)
+                    ∂ʸ = PartialDerivative{2}(𝐱)
+                    Pfx = ∂ˣ * Pf
+                    Pfy = ∂ʸ * Pf
+                    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