3D volume plot for spherical coordinate grid points

[raman_kumar]

How to get volume plot for the following spherical coordinate positions?
I have a lot of blocks each having such r, θ and ϕ spherical coordinate position vectors and i want to create 3D volume plot for each block vectors and finally combine them to create a single 3D volume plot.

r = Float32[1.677044, 1.8385518, 2.0156136, 2.2097273, 2.4225352, 2.6558375, 2.911608, 3.1920104, 3.499417, 3.8364286, 4.205896, 4.6109447, 5.055002, 5.5418243, 6.0755296, 6.660634, 7.3020864, 8.005314]
θ = Float32[0.049087387, 0.14726216, 0.24543692, 0.3436117]
ϕ = Float32[0.19634955, 0.5890486, 0.9817477, 1.3744467, 1.7671459, 2.1598449, 2.552544, 2.9452431, 3.3379421, 3.7306414, 4.12334, 4.5160394, 4.9087386, 5.3014374, 5.6941366, 6.086836]

[JM_Beckers]

What do you mean by “volume plot”? Your vectors do not have the same length, so what do they represent ?

@raman_kumar

x = Float32[1.677044, 1.8385518, 2.0156136, 2.2097273, 2.4225352, 2.6558375, 2.911608, 3.1920104, 3.499417, 3.8364286, 4.205896, 4.6109447, 5.055002, 5.5418243, 6.0755296, 6.660634, 7.3020864, 8.005314]
y = Float32[0.049087387, 0.14726216, 0.24543692, 0.3436117]
z = Float32[0.19634955, 0.5890486, 0.9817477, 1.3744467, 1.7671459, 2.1598449, 2.552544, 2.9452431, 3.3379421, 3.7306414, 4.12334, 4.5160394, 4.9087386, 5.3014374, 5.6941366, 6.086836]

[juliohm]

Here is one way to define a RectilinearGrid with Spherical coordinates:

using Meshes
using CoordRefSystems

import GLMakie as Mke

r = Float32[1.677044, 1.8385518, 2.0156136, 2.2097273, 2.4225352, 2.6558375, 2.911608, 3.1920104, 3.499417, 3.8364286, 4.205896, 4.6109447, 5.055002, 5.5418243, 6.0755296, 6.660634, 7.3020864, 8.005314]
θ = Float32[0.049087387, 0.14726216, 0.24543692, 0.3436117]
ϕ = Float32[0.19634955, 0.5890486, 0.9817477, 1.3744467, 1.7671459, 2.1598449, 2.552544, 2.9452431, 3.3379421, 3.7306414, 4.12334, 4.5160394, 4.9087386, 5.3014374, 5.6941366, 6.086836]

g = RectilinearGrid{𝔼,typeof(Spherical(0,0,0))}(r, θ, ϕ)

viz(g)

The type parameters in the constructor are a bit low-level, but I think they match what you have in mind. The first type parameter 𝔼 refers to the manifold where the geometries live. In this case this is the Euclidean manifold. The second type parameter is the coordinate reference system type, which I am producing with a typeof(Spherical(0,0,0)) call.

You can also define the periodicity of each coordinate. Below you can find an alternative code that considers one of the angular coordinates to be periodic. First, it creates a GridTopology with the number of elements in the grid and the periodicity information, and then it feeds the topology as the second argument of the RectilinearGrid:

# grid topology with given number of elements per dimension
# the last coordinate is periodic so we set (false, false, true)
t = GridTopology((length(r)-1, length(θ)-1, length(ϕ)-1), (false, false, true))

g = RectilinearGrid{𝔼,typeof(Spherical(0,0,0))}((r, θ, ϕ), t)

viz(g)

Notice how the cone is now full-filled compared with the previous one.

[juliohm]

I’ve updated the answer to include a variation that you might be interested in.

[raman_kumar]

@juliohm Thank you very much. I will also map Density on these 3D points.

@juliohm

Notice how the cone is now full-filled compared with the previous one.

For me first plot would be fine.

[juliohm]

You can map colors to every element of the grid with

viz(g, color=1:nelements(g))

Notice the length of the vector. It must match the number of elements in the grid, not the number of vertices. Some viz methods also support colors for every vertex, but I don’t remember if this is the case with RectilinearGrid. You can try and see.

[raman_kumar]

I got this error below. Please update you reply and replace first line with

using Meshes, CoordRefSystems, GLMakie

so that people don’t have misunderstanding while reading your solution later.

julia> g = RectilinearGrid{𝔼,typeof(Spherical(0,0,0))}(r, θ, ϕ)
ERROR: UndefVarError: Spherical not defined in Main
Suggestion: check for spelling errors or missing imports.
Hint: a global variable of this name may be made accessible by importing CoordRefSystems in the current active module Main

[JoshuaLampert]

@raman_kumar

Can many such block plots be combined to get final single plot?

Sure, just create another RectilinearGrid in the same way with the other coordinates and call viz! instead of viz to plot it in the same figure.

[raman_kumar]

Thanks to you all. I got this plot of Athena++ Black Hole simulation but color of density plot looks not so smooth and real.

and for alpha = 0.3 it looks

Is there any way to interpolate color to smooth them?

[sdanisch]

First thing should be to use transparency=true for the second image

[raman_kumar]

On using

viz!(g, color = nrho, transparency=true, alpha = 0.3)

i got this plot. But still it don’t look smooth.

and cone segment looks as shown below.

[juliohm]

We currently don’t forward the transparency option to the underlying Makie.jl function. It should be easy to submit a PR with a fix.

[sdanisch]

As a work around, if you can access the plot object from what viz! returns, you could do plot.transparency=true.
If it’s a recipe, you’ll find the mesh plot in plot.plots.

[raman_kumar]

scene = viz!(g, color = nrho, transparency=true, showsegments = false, alpha = 0.3)
scene.plots[1].transparency = true

On using above code i still get patchy rough plot.

[juliohm]

@raman_kumar by “patchy rough plot” you mean discrete color transitions? The color option is assigned to each element of the grid/mesh so this is an expected result. If the colors were assigned to vertices, the graphics card would perform an additional interpolation to make the final color map look “smooth”.

[raman_kumar]

Thanks, I should plot to vertices and not elements.

@juliohm

If the colors were assigned to vertices, the graphics card would perform an additional interpolation to make the final color map look “smooth”.

I have nelements(g) = 765 and nvertices(g) = 1152 but i have size(nrho) = (1152,).

• How to assign color to vertices and not elements?

[juliohm]

In theory all you need to do is pass the vector of colors with the correct length. Some viz recipes will have specializations for vertices, but I don’t think the recipe with RectilinearGrid has this specialization at the moment.

You can try with a CartesianGrid(10,10) to see the difference:

viz(CartesianGrid(10,10), color=1:100) # color elements
viz(Cartesiangrid(10,10), color=1:121) # color vertices

[sdanisch]

@raman_kumar

scene.plots[1].transparency = true

Hm, that still doesn’t look like propper transparency…
Is this GLMakie?

[raman_kumar]

Yes, I am using GLMakie.jl.

@raman_kumar

r = Float32[1.677044, 1.8385518, 2.0156136, 2.2097273, 2.4225352, 2.6558375, 2.911608, 3.1920104, 3.499417, 3.8364286, 4.205896, 4.6109447, 5.055002, 5.5418243, 6.0755296, 6.660634, 7.3020864, 8.005314]
θ = Float32[0.049087387, 0.14726216, 0.24543692, 0.3436117]
ϕ = Float32[0.19634955, 0.5890486, 0.9817477, 1.3744467, 1.7671459, 2.1598449, 2.552544, 2.9452431, 3.3379421, 3.7306414, 4.12334, 4.5160394, 4.9087386, 5.3014374, 5.6941366, 6.086836]

@juliohm Above are the r, θ and ϕ values of cell centers along x1/x2/x3-direction in mesh.
But when i supply in code different r, θ and ϕ denoting values of interface locations along x1/x2/x3-direction

r =  Float32[1.6, 1.754088, 1.9230156, 2.1082115, 2.311243, 2.5338273, 2.7778475, 3.0453684, 3.3386526, 3.6601818, 4.012676, 4.3991165, 4.8227735, 5.287231, 5.7964177, 6.354642, 6.9666257, 7.637547, 8.373081]
θ =  Float32[0.0, 0.09817477, 0.19634955, 0.2945243, 0.3926991]
ϕ =  Float32[0.0, 0.3926991, 0.7853982, 1.1780972, 1.5707964, 1.9634954, 2.3561945, 2.7488935, 3.1415927, 3.5342917, 3.9269907, 4.3196898, 4.712389, 5.105088, 5.497787, 5.8904862, 6.2831855]

then i get grid without cracks as in second plot of your post .

@juliohm

Notice how the cone is now full-filled compared with the previous one.

[juliohm]

If these are grid centers you will need some additional processing to convert into vertices for the RectilinearGrid