Colour Isosurfaces by axial position

[ElHouas]

Dear Pyvista Team,

I wanted to ask you how can I colour a isosurface with the values of the axial position.

I have been already able to extract the desired surface as you can see the in the image attached.

The problem is that I want to colour a surface with a certain value as the example in this link: https://www.researchgate.net/figure/Isosurfaces-of-vorticity-magnitude-coloured-by-axial-position-z-where-the-vorticity_fig1_278646580

I tried using the scalars parameters in several ways or the plot_curvature option with clim, but I am unsuccessful...

Here is my code:

# Create the spatial reference
grid = pv.UniformGrid()

grid.dimensions = np.array(Enstrophy.shape) + 1

grid.origin = (0, 0, 0)  # The bottom left corner of the data set
grid.spacing = (1, 1, 1)  # These are the cell sizes along each axis

grid.cell_arrays["values"] = Enstrophy.flatten(order="F")  # Flatten the array!

nodal = grid.cell_data_to_point_data()

geom = nodal.extract_geometry()

print(geom.number_of_points)

contours = nodal.contour(isosurfaces=5, scalars="axial")
pv.set_plot_theme("document")
contours.plot(cmap='Spectral')

cpos = contours.plot_curvature(curv_type='mean', clim=[0, 10] )

plotter = pv.Plotter()  # instantiate the plotter
plotter.add_mesh(contours)    # add a dataset to the scene
cpos = plotter.show()     # show the rendering window

plotter = pv.Plotter(off_screen=True)
plotter.add_mesh(contours, cmap='Spectral')
plotter.camera_position = cpos
plotter.show(auto_close=False)
plotter.screenshot('vorticity.png')
plotter.close()

Can you point me in the direction I need to work on to achieve this figure?

I really appreciate your time and attention.

Ismael

[banesullivan]

Hi @ElHouas! Thanks for posting your question! I should be able to help with this one - could you share with me a direct download link to the data file in this example? I can't seem to find the data file on research gate

[ElHouas]

Dear @banesullivan,

Really appreciate your quick response :)

I will attach you the data and the beginning of the code which does the preprocessing part:

import pyvista as pv
import numpy as np
import pickle
import matplotlib as mpl
import matplotlib.pyplot as plt

D = 1   # <!--- diameter of the wind turbine-->
Lx = 4 * D # <!---Length of the domain on x-->
Ly = 1.5 * D
Lz = 1.5 * D
nx = 415
ny = 145
nz = 145
dx = Lx/nx
dy = Ly/ny
dz = Lz/nz

uvec_file = open("U_0_415_snap600_steady.dms", "rb")
uvec_0_415_time = np.array(pickle.load(uvec_file))

vvec_file = open("V_0_415_snap600_steady.dms", "rb")
vvec_0_415_time = np.array(pickle.load(vvec_file))

wvec_file = open("W_0_415_snap600_steady.dms", "rb")
wvec_0_415_time = np.array(pickle.load(wvec_file))

x0 = np.linspace(0, 4, 415)     #Lx
x1 = np.linspace(0, 1.5, 145)   #Ly
x2 = np.linspace(0, 1.5, 145)   #Lz

x0grid, x1grid = np.meshgrid(x0, x2)

# <!---Compute vorticity and save in file-->

dw_dy = np.gradient(wvec_0_415_time, 2*dy, axis=1)
dv_dz = np.gradient(vvec_0_415_time, 2*dz, axis=0)
du_dz = np.gradient(uvec_0_415_time, 2*dz, axis=0)
dw_dx = np.gradient(wvec_0_415_time, 2*dx, axis=2)
dv_dx = np.gradient(vvec_0_415_time, 2*dx, axis=2)
du_dy = np.gradient(uvec_0_415_time, 2*dy, axis=1)

wx = dw_dy - dv_dz
wy = du_dz - dw_dx
wz = dv_dx - du_dy

Enstrophy = np.sqrt(wx**2 + wy**2 + wz**2)

The data files are big so I cannot attach them. Do you prefer an url to my github or one drive ?

Thanks :) :)

[banesullivan]

A URL to either works for me!

[ElHouas]

Here are the three links for the pickle files. Let me know if everything is okay

File U_0_415_snap600_steady: https://imperiallondon-my.sharepoint.com/:u:/g/personal/ie3518_ic_ac_uk/EV1489B-gvdOjiWXTyUeAbABCIsWDIcGNPsql_gFxTdiYw?e=8dgOmF

File V_0_415_snap600_steady: https://imperiallondon-my.sharepoint.com/:u:/g/personal/ie3518_ic_ac_uk/EaaBlu4YP_ZLugw3k9cTBaEBYrX_Y-gIXYOGZwIvdhp5Ew?e=hMeVUv

File W_0_415_snap600_steady: https://imperiallondon-my.sharepoint.com/:u:/g/personal/ie3518_ic_ac_uk/Ee9UhS-7lLxDoniA72W40AYBWwBopciYG9DK-X_0y0lUyw?e=eqvItK

Thanks a lot :)

[banesullivan]

Okay so normally this is quite easy with all mesh types - however, I just realized there's a pretty major bug with pyvista.UniformGrid.points where the points aren't ordered properly...

Normally, you would simply perform:

grid['axial'] = grid.points[:,2]
grid.plot(scalars='axial')

to have the grid colored by the z component of its spatial reference. Unfortunately, the UniformGrid class returns the points in the improper order for this to work with UniformGrid. This leaves you with two options until the next release (I fixed this in https://github.com/pyvista/pyvista/commit/13ae67df03ed055fce258f2d16a2f7f56bdf01c5):

1) casting to pyvista.UnstructuredGrid 2) using the .elevation() filter

Casting as an unstructured grid is costly and should be avoided. Try using the .elevation() filter which by default will add an array for the z component of the mesh which is what I believe you are attempting:

... # the code to create the dataset

# PyVista can handle F-ordered 3D NumPy arrays!
grid = pv.wrap(Enstrophy)
axial = grid.elevation()

# A nice camera position I found manually
cpos = [(446.9417095818895, -223.12196318769293, -327.0090049982658),
 (66.75124640066495, 95.17559002456603, 188.9825771285702),
 (-0.23712139635742802, -0.8951474486847606, 0.377471175724276)]

contours = axial.contour(isosurfaces=5, scalars='values')
contours.plot(scalars='Elevation', cmap='rainbow', cpos=cpos)

[banesullivan]

Also, this is a really cool dataset and is an excellent example for volume rendering:

grid.active_scalar_name = 'values'
plotter = pv.Plotter(notebook=False)
plotter.add_volume(grid, scalars='values', cmap='viridis', opacity='sigmoid_6')
plotter.show(cpos=cpos)

[banesullivan]

And here's a neat way to get a profile view (I took log of the values here):

grid['log(v)'] = np.log10(grid['values'])
slc = grid.slice()

plotter = pv.Plotter(shape=(2,1), notebook=False)
plotter.add_mesh(slc, 
                 scalars='log(v)', 
                 cmap='viridis', 
                 lighting=False)
plotter.view_vector([1,0,0], viewup=[0, 1, 0])

plotter.subplot(0,1)
plotter.add_mesh(slc, 
                 scalars='values', 
                 cmap='viridis', 
                 lighting=False)

plotter.link_views()
plotter.show()

[banesullivan]

@ElHouas - this dataset is quite cool and useful for examples - would it be okay for us to include this or another one similar to it (that you've made) in our examples module?

[ElHouas]

Hi @banesullivan !!

Yesterday I was so tired from work and I went sleep directly, sorry for the late answer.

Once, I come from work I will continue working on my thesis and try what you explained me. I will also ask my supervisor first, in order to avoid misunderstandings. I will be very happy to share these examples with the community and help that faced the same problems.

I like the features that Pyvista offers, so I will continue exploring its capabilities. I would let you know if I have any new issue.

I was struggling with it, so I really appreciate your help :) :)

Have a nice day!