Triangle examples in examples/triangleexamples.jl fail

[DanielVandH]

The examples here fail (below is Lines 1-12 from that file)

julia> using ApproxFun, MultivariateOrthogonalPolynomials, BlockArrays, SpecialFunctions, FillArrays, Plots
       import ApproxFun: blockbandwidths, Vec, PiecewiseSegment
       import MultivariateOrthogonalPolynomials: DirichletTriangle
       S = TriangleWeight(1,1,1,JacobiTriangle(1,1,1))
       Δ = Laplacian(S)
WARNING: could not import ApproxFun.blockbandwidths into Main
WARNING: could not import ApproxFun.Vec into Main
WARNING: could not import MultivariateOrthogonalPolynomials.DirichletTriangle into Main
ERROR: MethodError: no method matching TriangleWeight(::Int64, ::Int64, ::Int64, ::JacobiTriangle{Float64, Int64})

Closest candidates are:
  TriangleWeight(::T, ::T, ::T) where T
   @ MultivariateOrthogonalPolynomials C:\Users\User\.julia\packages\MultivariateOrthogonalPolynomials\dZw1Q\src\triangle.jl:111

Stacktrace:
 [1] top-level scope
   @ Untitled-1:4

I believe in this case the fix is to do S = WeightedTriangle(1, 1, 1)? This error appears later in the script as well (might be some others too - Laplacian(S) also errors -, but I haven't run it that far). I could try and fix it up later once I've looked into it all a bit more.

[DanielVandH]

Perhaps as a simple first step it would be useful to just run these scripts as parts of the CI to just see if they run completely, e.g. in runtests.jl

dir = readdir("./examples")
filter!(file -> endswith(file, ".jl"), dir)
for example_path in dir
    script = joinpath("./examples", example_path)
    mod = @eval module $(gensym()) end
    Base.include(mod, script) # make sure each script is self-contained
end

unless you have something else you usually do instead for this.

[dlfivefifty]

That’s a great idea. Can you add a PR?

[dlfivefifty]

Note that’s using the old ApproxFun version. We can probably change it to Laplacian(axes(S,1))

we should also add laplacian(S) which would match the diff(S) in 1D

[DanielVandH]

PR made. Yeah, I was trying to convert it to the newer version, and got to here before making this issue:

using MultivariateOrthogonalPolynomials, SpecialFunctions
a, b, c = 1, 1, 1
S = WeightedTriangle(a, b, c)
xy = axes(S, 1)
x, y = first.(xy), last.(xy)
∆ = Laplacian(xy)
f = @. 1.0 + erf(5.0(1.0 - 10.0((x - 1/2)^2 + (y-1/2)^2)))
N = 500
M = sparse(∆[Block.(1:N), Block.(1:N)])

StackOverflowError:

Stacktrace:
     [1] unblock(A::QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}, inds::Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}, I::Tuple{BlockArrays.BlockRange{1, Tuple{UnitRange{Int64}}}})
       @ BlockArrays C:\Users\User\.julia\packages\BlockArrays\K0HfQ\src\views.jl:10
     [2] to_indices(A::QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}, inds::Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}, I::Tuple{BlockArrays.BlockRange{1, Tuple{UnitRange{Int64}}}})
       @ ContinuumArrays C:\Users\User\.julia\packages\ContinuumArrays\Uyiyf\src\ContinuumArrays.jl:81
     [3] to_indices(A::QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}, I::Tuple{BlockArrays.BlockRange{1, Tuple{UnitRange{Int64}}}})
       @ ContinuumArrays C:\Users\User\.julia\packages\ContinuumArrays\Uyiyf\src\ContinuumArrays.jl:85
     [4] _getindex(::Type{Tuple{StaticArraysCore.SVector{2, Float64}}}, A::QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}, I::Tuple{BlockArrays.BlockRange{1, Tuple{UnitRange{Int64}}}})
       @ QuasiArrays C:\Users\User\.julia\packages\QuasiArrays\ZCn0F\src\abstractquasiarray.jl:372
     [5] getindex(A::QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}, I::BlockArrays.BlockRange{1, Tuple{UnitRange{Int64}}})
       @ QuasiArrays C:\Users\User\.julia\packages\QuasiArrays\ZCn0F\src\abstractquasiarray.jl:367
--- the last 5 lines are repeated 15995 more times ---
 [79981] unblock(A::QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}, inds::Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}, I::Tuple{BlockArrays.BlockRange{1, Tuple{UnitRange{Int64}}}})
       @ BlockArrays C:\Users\User\.julia\packages\BlockArrays\K0HfQ\src\views.jl:10
 [79982] to_indices(A::QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}, inds::Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}, I::Tuple{BlockArrays.BlockRange{1, Tuple{UnitRange{Int64}}}})
       @ ContinuumArrays C:\Users\User\.julia\packages\ContinuumArrays\Uyiyf\src\ContinuumArrays.jl:81
 [79983] to_indices(A::QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}, I::Tuple{BlockArrays.BlockRange{1, Tuple{UnitRange{Int64}}}})
       @ ContinuumArrays C:\Users\User\.julia\packages\ContinuumArrays\Uyiyf\src\ContinuumArrays.jl:85

What's the correct code for reproducing Example 1 in the paper, i.e. $u_{xx} + u_{yy} = f(x, y)$ for $(x, y) \in T$ with $\left.u\right|_{\partial T} = 0$, where $f(x, y) = 1 + \mathrm{erf}(5(1 - 10((x - 1/2)^2 + (y-1/2)^2)))$?

[dlfivefifty]

Hmm I just checked and I realised that the code for partial derivatives for WeightedTriangle isn't implemented yet. So actually the first task is to get that working. That is, we want to first add a function like:

@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))
    M =  _BandedBlockBandedMatrix((k .+ convert(T, b+c))', axes(k,1), (-1,1), (-1,1)) # TODO: replace with correct formula
    WeightedTriangle(a,b-1,c-1) * M # TODO: deal with the case where b or c are 0
end