JuliaApproximation
Triangle examples in examples/triangleexamples.jl fail
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.
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.
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
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)))$?
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