Vedo: How to plot streamlines

[Yoorthiziri]

Hello,
Sorry, if this is not the right place to post this.
I followed the different examples provided by vtkplotter to plot streamlines. Unfortunately I couldn't make it work.
Here is an example of a vtk file cavity.zip that contains the velocity vector field U.
Here is my code:

"""
Draw streamlines for the cavity case from OpenFOAM tutorial
"""

from vtkplotter import loadUnStructuredGrid, load, show, interactive, Box, Grid, streamLines, Text

cavity = load('cavity.vtk')

domain = load('cavity.vtk')

comment = Text(__doc__, c='w')
probe = Grid(pos=[0.05,0.1,0.005], normal=[0,1,0], sx=0.01, sy=0.1, resx=30, resy=10)
probe.color('k')

stream = streamLines(domain, probe, direction='forward')
show(domain, stream, probe, comment, axes=8)

when I run the example above, I couldn't see the streamlines
I appreciate your help

[marcomusy]

Hi @Yoorthiziri
this is the right place to ask this type of question :)
please upgrade vtkplotter with the command:
pip install git+https://github.com/marcomusy/vtkplotter.git --upgrade
then you should be able to run:

"""Draw streamlines for the cavity case from OpenFOAM tutorial"""

from vtkplotter import *

# load file as type vtkUnStructuredGrid
ugrid = loadUnStructuredGrid('cavity.vtk')

# make a grid of points to probe as type vtkActor
probe = Grid(pos=(0.05,0.08,0.005), normal=(0,1,0),
             sx=0.01, sy=0.1, resx=4, resy=20, c='k')

# compute stream lines with Runge-Kutta4, return a vtkActor
stream = streamLines(ugrid, probe,
                     activeVectors='U', # name of the active array
                     #tubes={"radius":1e-04, "varyRadius":2},
                     lw=2, # line width
                     )

# make a cloud of points form the ugrid, in order to draw arrows
domain = pointCloudFrom(ugrid)
coords = domain.getPoints()
vects  = domain.getPointArray('U')/200
arrows = Arrows(coords-vects, coords+vects, c='jet_r') # use colormap
box    = domain.box().c('k') # build a box frame of the domain
comment= Text(__doc__)

show(stream, arrows, box, probe, comment, axes=5, bg='white')

Let me know if you encounter any problems running it.

PS; may I add this example to the gallery of examples? Thanks for your feedback.

[Yoorthiziri]

Thank you very much for the detailed answer.
Of course you could add it to the examples. You could find many others in the tutorials directory of OpenFOAM, just use foamToVTK to convert the results to vtk files.

I have just another question about how to choose the field to plot if you have more than one. For example if I convert the original OpenFOAM cavity results, the final VTK will contain the field p for pressure and U for velocity. In this case how can I choose which field to plot in vtkplotter?

Thanks

[marcomusy]

Cool! I guess that setting activeVectors='U', # name of the active array in the streamLines method should work?

[Yoorthiziri]

Yes you are right about that for streamlines.
I don't think that this question belongs to this issue but just to clarify: I mean if I have a vtk file, example cavity.vtk that has multiple fields (p, U, ...etc) when I run it like this:

vtkplotter cavity.vtk

By default it show the pressure controus, hw can I switch to another field from within vtkplotter window? is there any shortcut, command line argument for that?

Thank you

[marcomusy]

Yes, press 4 to cycle through existing arrays.
If you press h, you will get the full list of key options.

[Yoorthiziri]

Thank you very much.

If you press h, you will get the full list of key options.

Unfortunately the description in the help is not clear, here is the description:

4 use scalars as colors, if present
If I am not mistaken, this is not precise because vector fields are not scalar.
Thanks

[marcomusy]

You are right, although what is visualized should be the magnitude of the vector.

[Yoorthiziri]

Hello again!
if the mesh is very dense, how can we plot the velocity vectors at only some locations (that is, by skipping some points)?

Another question, how could we plot the velocity vectors at thr cell centers? instead of using cell vertices?

Thank you

[marcomusy]

Hi, you could do something like this:

• clean() imposes a minimum tolerance between close points (in percent of the bounding box)
• cutWithPlane() removes one plane of points (but any numpy method here would also work fine), and z() shifts the whole set of arrows.

from vtkplotter import *

# load file as type vtkUnStructuredGrid
ugrid = loadUnStructuredGrid(datadir+'cavity.vtk')

# make a grid of points to probe as type Mesh(vtkActor)
probe = Grid(pos=(0.05,0.08,0.005), normal=(0,1,0),
             sx=0.1, sy=0.01, resx=20, resy=4, c='k')

# compute stream lines with Runge-Kutta4, return a Mesh(vtkActor)
stream = streamLines(ugrid, probe,
                     activeVectors='U', # name of the active array
                     #tubes={"radius":1e-04, "varyRadius":2},
                     lw=2, # line width
                     )

# make a cloud of points form the ugrid, in order to draw arrows
domain = pointCloudFrom(ugrid)
box = domain.box().c('k') # build a box frame of the domain

domain.clean(tol=0.07).cutWithPlane([0,0,0.0001], normal=(0,0,1)).z(-0.005)

coords = domain.points()
vects  = domain.getPointArray('U')/50
arrows = Arrows(coords-vects, coords+vects, c='jet_r') # use colormap

show(stream, arrows, box, probe, axes=5)

[Yoorthiziri]

Thank you very much for your answer.

Regarding plotting vectors at cell centers, I am trying to use the following:

 # make a cloud of points form the ugrid, in order to draw arrows
domain = pointCloudFrom(ugrid).mapPointsToCells()
coords = domain.cellCenters()
vects  = domain.getCellArray('U')/200
arrows = Arrows(coords-vects, coords+vects, c='jet_r') # use colormap
box    = domain.box().c('k') # build a box frame of the domain
show(stream, arrows, box, probe, axes=5)

Is that the correct way to plot the vector at the cell centers? I am not able to show the mesh cells to see if the vectors are plotted at the centers because when I press 'w' to show the wireframe, it toggles the wireframes/surface for the vectors instead of the domain.

Thank you

[Yoorthiziri]

Sorry if my question was not clear
I mean the velocity vectors at the cell centers. For example something like this:

[marcomusy]

..then you need to do something more complicated.. because the vector field is defined at the mesh vertices so you need to interpolate it to a different set of points (the cell centers).
This can be controlled in various ways, one possibility is to use RBF (radial basis functions) from scypi:

from vtkplotter import *

# load file as type vtkUnStructuredGrid
ugrid = loadUnStructuredGrid(datadir+'cavity.vtk')

#################################################
# use native vtk commands to get cell centers
from vtk import vtkCellCenters
from vtk.util.numpy_support import vtk_to_numpy
vcen = vtkCellCenters()
vcen.SetInputData(ugrid)
vcen.Update()
cns = vtk_to_numpy(vcen.GetOutput().GetPoints().GetData())

#################################################
# make a grid of points to probe as type Mesh(vtkActor)
probe = Grid(pos=(0.05,0.08,0.005), normal=(0,1,0),
             sx=0.1, sy=0.01, resx=20, resy=4, c='k')

# compute stream lines with Runge-Kutta4, return a Mesh(vtkActor)
stream = streamLines(ugrid, probe,
                     activeVectors='U', # name of the active array
                     #tubes={"radius":1e-04, "varyRadius":2},
                     lw=2, # line width
                     )

# make a cloud of points form the ugrid, in order to draw arrows
domain = pointCloudFrom(ugrid)
box = domain.box().c('k') # build a box frame of the domain

coords = domain.points()
vects  = domain.getPointArray('U')/50


################################################# scipy RBF
from scipy.interpolate import Rbf
import numpy as np

cns = Points(cns).clean(tol=0.07).points() # reduce the density

xr, yr, zr = cns[:, 0], cns[:, 1], cns[:, 2]
x,  y,  z  = coords[:, 0], coords[:, 1], coords[:, 2]
dx, dy, dz = vects[:, 0], vects[:, 1], vects[:, 2]

itrx = Rbf(x, y, z, dx)  # Radial Basis Function interpolator:
itry = Rbf(x, y, z, dy)  #  interoplate the deltas in each separate
itrz = Rbf(x, y, z, dz)  #  cartesian dimension

positions_x = itrx(xr, yr, zr) + xr
positions_y = itry(xr, yr, zr) + yr
positions_z = itrz(xr, yr, zr) + zr
positions_rbf = np.c_[positions_x, positions_y, positions_z]

arrows = Arrows(cns, positions_rbf, c='jet_r')

show(stream, arrows, box, probe, axes=5)

PS: the script you show cannot work because the cell positions corresponding to a point cloud is still the same point cloud..

[Yoorthiziri]

Thank you very much. I appreciate your help.