Question on 3D plotting capabilities

[jlchan]

Hi @j-fu - I'm one of the developers for the high order CFD package Trixi.jl, and we've been looking for a way to visualize 3D flows on unstructured meshes natively in Julia. GridVisualize.jl looks very promising!

I wanted to ask about a few things:

• can GridVisualize.jl visualize 3D scalar data and isosurfaces on unstructured tetrahedral meshes with GLMakie.jl? It seems that way from the docs, but I wanted to double check.
• Would it be computationally feasible to use this to visualize isosurfaces around 10-100 million tetrahedra?
• Does it complicate things to visualize solutions which are discontinuous across element interfaces?

Thank you!

[j-fu]

Hi @j-fu - I'm one of the developers for the high order CFD package Trixi.jl, and we've been looking for a way to visualize 3D flows on unstructured meshes natively in Julia. GridVisualize.jl looks very promising!

I wanted to ask about a few things:

• can GridVisualize.jl visualize 3D scalar data and isosurfaces on unstructured tetrahedral meshes with GLMakie.jl? It seems that way from the docs, but I wanted to double check.

Yes, this is the intention.

• Would it be computationally feasible to use this to visualize isosurfaces around 10-100 million tetrahedra?

Essentially, yes. Given a piecewise linear continuous function, a set of isosurface values, a GeometryBasics.Mesh is created by a marching tetrahedra algorithm, which is then enhances with normals and visualized. The later step is really fast with GLMakie, and GLMakie it is the only package I am aware of which allows to do this. It uses modern OpenGL for this purpose, which is the right thing to do.

It boils down to a scaling test which would be interesting to do anyway.

• Does it complicate things to visualize solutions which are discontinuous across element interfaces?

Just a little bit - this would amount to add something like a "brokenscalarplot" method to the UI. However I would say that for users this may be not very intuitive, and one could think about plotting a continuous reconstruction.

On a different note, AFAIU, Trixi.jl uses octree like grids with adaptive local refinement. In this case it may be worth to think about plotting directly from that data structure. This mainly would amount to adapting marching_tetrahedra (see https://github.com/j-fu/GridVisualize.jl/blob/0bf2c6690bc1d43e661b92b9814dbf81f429debd/src/common.jl#L294 ) .

[j-fu]

Another possibility would be to add an API based on an iterator which produces cell data by demand, some generalization of what I did in my old code https://github.com/j-fu/gltools .

[jlchan]

Essentially, yes. Given a piecewise linear continuous function, a set of isosurface values, a GeometryBasics.Mesh is created by a marching tetrahedra algorithm, which is then enhances with normals and visualized. The later step is really fast with GLMakie, and GLMakie it is the only package I am aware of which allows to do this. It uses modern OpenGL for this purpose, which is the right thing to do.

Thanks - that sounds ideal. I haven't found any other packages which implement marching tetrahedra for unstructured meshes (I believe https://github.com/JuliaGeometry/Meshing.jl only implements it for structured data or signed distance functions).

• Does it complicate things to visualize solutions which are discontinuous across element interfaces?

Just a little bit - this would amount to add something like a "brokenscalarplot" method to the UI. However I would say that for users this may be not very intuitive, and one could think about plotting a continuous reconstruction.

That would certainly be possible on our end too. Visualizing discontinuous functions would mostly be useful for debugging. Let me try out GridVisualize.jl in more detail before exploring this option.

On a different note, AFAIU, Trixi.jl uses octree like grids with adaptive local refinement. In this case it may be worth to think about plotting directly from that data structure.

That's one possibility, but we've been extending our solvers to curved/unstructured high order meshes. We currently use interpolation to a fine mesh for visualizing curved domains and high order solutions in 2D, and were hoping to extend this to 3D using tetrahedralizations.

[j-fu]

Ok, that is a good way to start. For a way to get an ExtendableGrid from TetGen, see the code from https://github.com/j-fu/SimplexGridFactory.jl/blob/master/src/tetgen.jl#L9

I am aware of Meshing.jl but did not find the time yet to lift my marching_tetrahedra to the level of abstraction seemingly appropriate for that package...

[jlchan]

For reference, here are a few timings:

julia> g,f=func3d(n=50); @time scalarplot(g,f; flevel=.25, Plotter=GLMakie)
  1.202893 seconds (4.35 M allocations: 158.493 MiB, 29.25% gc time)

julia> g,f=func3d(n=100); @time scalarplot(g,f; flevel=.25, Plotter=GLMakie)
  3.937516 seconds (6.64 M allocations: 381.042 MiB, 9.93% gc time)

julia> g,f=func3d(n=150); @time scalarplot(g,f; flevel=.25, Plotter=GLMakie)
 12.738594 seconds (10.46 M allocations: 752.109 MiB, 9.72% gc time)

julia> g,f=func3d(n=200); @time scalarplot(g,f; flevel=.25, Plotter=GLMakie)
 26.975054 seconds (15.80 M allocations: 1.221 GiB, 6.73% gc time)

The setup for n=100 has about 1 million nodes, 6 million cells, whilen=200 has about 8 million nodes and 48 million cells.

[j-fu]

Ok for n=100 I see on my laptop:

 @time scalarplot(g,f; Plotter=GLMakie, flevel=0.25,outline=false)
  1.013448 seconds (51 allocations: 24.516 MiB)
  1.069589 seconds (51 allocations: 24.516 MiB)
  2.483476 seconds (2.51 M allocations: 225.351 MiB)

The intermediate timings are for the marching tetrahedra which are called twice as it appears - I would need to figure out why - things are updated by Makie via observables, so it is a bit intricate to figure out. I however think it is possible.

The timings for marching tetrahedra show that this code is nearly allocation free (I cared about this), so it is pure loops/calculations where I am not sure how much can be done. Some loop unrolling probably will give some speedup.

Besides of that it could be multithreaded in a straightforward manner.

[jlchan]

One other possibility would be using lower precision. I admit I haven't looked through the rest of the code in detail so I'm not sure how hard it would be to introduce Float32 or Float16.

[j-fu]

One other possibility would be using lower precision. I admit I haven't looked through the rest of the code in detail so I'm not sure how hard it would be to introduce Float32 or Float16.

As expected, this does not help much (and I already assemble stuff in float32).

[j-fu]

For the lift, see https://github.com/JuliaPlots/Makie.jl/issues/1303 ...

[j-fu]

For the lift, see JuliaPlots/Makie.jl#1303 ...

... tagging a workaround...

[j-fu]

I also have something with vtk, (VTKView.jl in a separate test registry) see https://www.youtube.com/watch?v=DmueA_Lvigs

Got it working only on linux (vtk compilation...), it can run as a backend to GridVisualize as well.

I is not significantly faster for this problem (no wonder as its pipelines also have some overhead).

So multithreading...

[j-fu]

Ah I see that I am always calculating plane intersections though the planes are not visible by default...

[j-fu]

Ok with 0c2cc0fffa86d9b6aa493ec316c6931c28d3c5e7 shifted the blame to GLMakie... marching tets for 1M nodes take now 0.1s (without threading yet). The rest of the time is GLMakie, and a bit too many allocations over there IMHO.

[jlchan]

That's an incredible speed up! That should be more than enough to work for our settings.

[j-fu]

You can speed it up a bit more with outline=false. This prevents the grid boundary elements from being rendered.

[j-fu]

Continued in #8...