paraview.simple: get points coordinates from a slice

ParaView Support
1

Hello,

I’m trying to obtain the points coordinates of a slice filter from OpenFOAM data without luck.
Here what I’ve tried so far

Extract data from a plane example (file is also attached) test_mpl.py (2.0 KB)

import os

import numpy as np
import matplotlib.pyplot as plt

from paraview.simple import *
from vtk.util import numpy_support as npvtk

def plane_points(plane_data, verbose=False):
“”"Extract the Cartesian points coordinates (x,y,z) of the triangulated plane

Parameters
----------
plane_data: Paraview slice object

Return
------
x,y,z: arrays of floats
       Cartesian coordinates of the points in the plane
"""

# Get the number of points
nb_points = plane_data.GetNumberOfPoints()

# Put the points coordinates in x, y and z arrays
x_list = []
y_list = []
z_list = []

for i in range(nb_points):
    coord = plane_data.GetPoint(i)
    xx, yy, zz = coord[:3]
    x_list.append(xx)
    y_list.append(yy)
    z_list.append(zz)

x = np.array(x_list)
y = np.array(y_list)
z = np.array(z_list)

# Display some information on the triangulation points
if verbose:
    print("Number of points:", nb_points)
    print("xrange: [{0:5.3f}, {1:5.3f}]".format(x.min(), x.max()))
    print("yrange: [{0:5.3f}, {1:5.3f}]".format(y.min(), y.max()))
    print("zrange: [{0:5.3f}, {1:5.3f}]".format(z.min(), z.max()))

# Return the points coordinates of the triangulated plane
return x, y, z

create a new ‘OpenFOAMReader’

pipeCyclicfoam = OpenFOAMReader(FileName=’./pipeCyclic/pipeCyclic.foam’)
pipeCyclicfoam.MeshRegions = [‘internalMesh’]
pipeCyclicfoam.CellArrays = [‘U’, ‘epsilon’, ‘k’, ‘nut’, ‘p’]

create a new ‘Slice’

slice1 = Slice(Input=pipeCyclicfoam)
slice1.SliceType = ‘Plane’
slice1.HyperTreeGridSlicer = ‘Plane’
slice1.SliceOffsetValues = [0.0]

Properties modified on slice1.SliceType

slice1.SliceType.Origin = [5.0, 0.0, 0.0]
slice1.SliceType.Normal = [1.0, 0.0, 0.0]

plane_data = servermanager.Fetch(slice1)

Extract the cartesian coordinates of the triangulation points

x,y,z = plane_points(plane_data, verbose=False)

I got this error message:

in plane_points(plane_data, verbose)
     31 
     32     for i in range(nb_points):
---> 33         coord = plane_data.GetPoint(i)
     34         xx, yy, zz = coord[:3]
     35         x_list.append(xx)

AttributeError: 'vtkmodules.vtkCommonDataModel.vtkMultiBlockDataSet' object has no attribute 'GetPoint'

Any idea ?
I’m stuck :frowning:

Many thanks !

2

You’re working with a MultiBlockDataSet so you’ll have to either:

  1. Work down to a leaf node using mbds.GetBlock(i) and then use leaf_node.GetPoint(i) or create an iterator to more easily get to the leaf nodes, something like:
itero = plane_data.NewIterator()
itero.SetVisitOnlyLeaves(1)
itero.InitTraversal()
while not itero.IsDoneWithTraversal():
    obj = itero.GetCurrentDataObject()
    nb = obj.GetNumberOfPoints()
    for i in range(nb):
        obj.GetPoint(i)
        ....
    itero.GoToNextItem()
  1. Use a MergeBlocks filter before Fetch to create a single dataset from your MBDS, which you can then directly operate on as you were trying before
3

It works, many many thanks !