Fix plotting

[dlfivefifty]

Plotting on the disk using Plots.jl is basically broken:

I suspect the best thing to do is forget about Plots.jl and focus on Makie.jl which I believe has advanced a lot. The catch is we probably don't want to make Makie.jl a dependency and there isn't a lightweight recipes package like in Plots.jl.

We also need to get it to use inverse transforms to compute the values.

[marcusdavidwebb]

Cool tornado!

[dlfivefifty]

Makie.jl actually has good plotting on a disk via just surface, which accepts matrices for x and y.

Unfortunately it doesn't work for contourf. My solution for triangle plotting was to generate a triangular plotting mesh by hand, the old code below supported both triangles and rectangles this way.

@ioannisPApapadopoulos What did you use for plotting with Firedrake? What would you do to visualise solutions to PDEs in a 3D ball?

using ApproxFun, Makie, MultivariateOrthogonalPolynomials
    import AbstractPlotting: mesh, mesh!
    import ApproxFun: _padua_length

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::Fun{<:Space{<:Triangle}}) =
    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

[ioannisPApapadopoulos]

In Firedrake you can save solutions as pvd files which can be opened in Paraview. Then in paraview there are all sorts of options of how to view the solutions, including rotations, filters, transparency of the solution at certain values etc. https://www.paraview.org/

[ioannisPApapadopoulos]

Here is a Julia package to save solutions suitable for Paraview. https://github.com/jipolanco/WriteVTK.jl

[TSGut]

I've extensively used Makie.jl in the past for 3D plotting of solutions to Schrödinger equations and it was not the best experience, though as you say it may have advanced significantly since then.

My main concern for publication-quality plotting is that Makie only supports vector graphics output for dimensions <= 2 via CairoMakie (the documentation suggests this is still the case), meaning that 3D graphics will be stuck in bitmap format via GLMakie.@ioannisPApapadopoulos I imagine Paraview allows exporting 3D plots as PDF or SVG?

[ioannisPApapadopoulos]

Yes, you can export in PDF, SVG, EPS, etc.

[TSGut]

For what it's worth, contour plotting is easy to get working (with the pyplot Plots.jl backend at least) from the current state.

pyplot()
Z = Zernike(1)
xy = axes(Z,1)
x = first.(xy)
y = last.(xy)
f = Z * (Z \ (y.^2 .+ exp.(x.^2 .+ y)))

grid = MultivariateOrthogonalPolynomials.plotgrid(f)
fl(x,y) = f[(x,y)]

X = getindex.(grid,1,1)
Y = getindex.(grid,2,1)
Z = map(fl,X,Y)

desiredlevels = 50
contour(X, Y, Z, legend=true, fill=true, aspect_ratio=:equal, levels=desiredlevels)

I didn't give this much thought but given that it works I suspect Plots.jl functionality could probably be restored if we want it.

[dlfivefifty]

I actually wrote VTK code 20 years ago.... For VTK do we have to do something similar as above and create the mesh ourselves? That is, provide the grid points but also how they connect together?

I noticed Gridap.jl also uses writevtk. Perhaps also Meshes.jl is useful. It seems like what we really want to do is just generate a mesh which can then be converted to VTK or Makie.jl but not sure the infrastructure is there yet.

But Timon's pyplot suggest looks good for now.

[MikaelSlevinsky]

Makie is the fastest but the PlotlyJS Plots backend looks the best to me. Lossless exporting as HTML https://juliaapproximation.github.io/FastTransforms.jl/dev/generated/disk/

[dlfivefifty]

Note there's now WGLMakie.jl that supports WebGL

[dlfivefifty]

And https://makie.juliaplots.org/dev/documentation/backends/rprmakie/ which will make our plots Pixar-quality

[dlfivefifty]

Actually using pyplot() seems to work as-is both with surface and contourf, here's the tornado:

The only issue is we aren't using the inverse transform so it's really slow af rn.