JuliaApproximation
Fix plotting
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.
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
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/
Here is a Julia package to save solutions suitable for Paraview. https://github.com/jipolanco/WriteVTK.jl
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?
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.
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.
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/
And https://makie.juliaplots.org/dev/documentation/backends/rprmakie/ which will make our plots Pixar-quality
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.