Getting node IDs from a vtkCell?

Support
1

I’m generating a text list of node connectivity in a dataset, getting the node locations is pretty straightforward, but the connectivity is less so in my particular case.

The dataset is a vtkUnstructuredGrid composed of vtkHexahedrons. The required connectivity is not the hexahedrons themselves, but the faces of the hexahedrons. Right now I am getting this list like this:

# Descriptive constants rather than magic numbers.
NFACES = ds.GetCell(0).GetNumberOfFaces() # I know it is 6 ... always
NNODES = ds.GetCell(0).GetFace(0).GetNumberOfPoints() # I know it is 4 ... always

connectivity = []
for idx in range(ds.GetNumberOfCells()):
    c = ds.GetCell(idx)
    for k in range(NFACES):
        f = c.GetFace(k)
        connectivity.append([f.GetPointId(g) for g in range(NNODES)])

Is there a more elegant way of getting the node list? The connectivity array for the hexahedrons is so much nicer-looking:

vtk_to_numpy(ds.GetCells().GetConnectivityArray())
2

I can’t think of an elegant way to do this in VTK, but there are more efficient ways. My first thought when I read your question was that there must be a simple mapping between the cell connectivity array and the face connectivity array. And since numpy can map arrays to arrays, it should be possible to get numpy to do the work.

So without further ado, here’s the lookup table (expressed as an array) that maps from face point IDs to cell point IDs for a vtkHexahedron. Each row gives the mapping for one of the 6 faces. Later on, I’ll show how this LUT is constructed.

[[0 4 7 3]
 [1 2 6 5]
 [0 1 5 4]
 [3 7 6 2]
 [0 3 2 1]
 [4 5 6 7]]

With this LUT, we only have to loop over the cells. We no longer loop over the faces or the nodes, since numpy takes care of those. This makes the code very efficient.

# lut is initialized to the lookup table shown above
# fca is face connectivity array
# cca is cell connectivity array
for i in range(ds.GetNumberOfCells()):
    fca[i*6:(i+1)*6,:] = cca[lut]
    lut += 8

Here’s a complete example that builds the lookup table and applies it to a synthetic data set:

import vtk
import numpy as np

# build a 3D cell
cell = vtk.vtkHexahedron()

# get 3D cell info
NVERTS = cell.GetNumberOfPoints()
NFACES = cell.GetNumberOfFaces()
NNODES = cell.GetFace(0).GetNumberOfPoints()

# set cell IDs to identity (this is necessary for lookup table)
for i in range(NVERTS):
    cell.GetPointIds().SetId(i, i)

# make a lookup table to convert cell connectivity to face connectivity
lut = np.zeros((NFACES, NNODES), dtype=np.int64)
for i in range(NFACES):
    face = cell.GetFace(i)
    for j in range(NNODES):
        lut[i,j] = face.GetPointId(j)

# make a sample data set
source = vtk.vtkCellTypeSource()
source.SetCellType(vtk.VTK_HEXAHEDRON)
source.SetBlocksDimensions(4,2,1)
source.Update()
ds = source.GetOutput()

# get cell connectivity array as numpy array
cca = np.frombuffer(ds.GetCells().GetConnectivityArray64(), dtype=np.int64)

# use lut to generate face connectivity array (fca) from cell connectivity 
# (this is where all the work is done)
NCELLS = cca.size//NVERTS
fca = np.zeros((NCELLS*NFACES, NNODES), dtype=np.int64)
for i in range(0, NCELLS*NFACES, NFACES):
    # get all faces for this cell
    fca[i:i+NFACES,:] = cca[lut]
    # advance lookup table to next cell
    lut += NVERTS

#print(fca)
3

I sort of follow - but I am unclear on how the lookup table works, it doesn’t appear in the example that it is being used to look things up?

4

The operation cca[lut] is doing the lookup. Consider just the first row of lut:

[[0 4 7 3] ...

Doing cca[lut] creates a new array equal to this

[[ cca[0], cca[4], cca[7], cca[3] ] ...

The new array has the same shape as lut, that is, it has one row for every face of the cell (I’ve only shown the first face, but it works the same for all the other faces, too). The lookup table knows which cell nodes correspond to which face nodes.

The reason that I used vtk.vtkHexahedron() to create a new cell, rather than using an existing cell from the dataset, is that generation of the lookup table requires a cell where all the IDs are set to identity. Generically, I could write this:

cellType = vtk.VTK_HEXAHEDRON

cell = vtk.vtkGenericCell.InstantiateCell(cellType)
for i in range(cell.GetNumberOfPoints()):
    cell.GetPointIds().SetId(i, i)

lut = np.zeros((NFACES, NNODES), dtype=np.int64)
for i in range(NFACES):
    face = cell.GetFace(i)
    for j in range(NNODES):
        lut[i,j] = face.GetPointId(j)

Note that a Python generator expression can be used to create an iterator that will use the lookup table to map the cell connectivity to face connectivity. This is perhaps a bit more elegant than the for() loop that I had at the end of my previous post.

# return an iterator for all faces in dataset
def faces_iterator(ds):
    cca = np.frombuffer(ds.GetCells().GetConnectivityArray64(), dtype=np.int64)
    totalFaces = NFACES * ds.GetNumberOfCells()
    return (cca[lut[i%NFACES] + i//NFACES*NVERTS] for i in range(totalFaces))

# example: iterate through faces
for face in faces_iterator(ds):
    print(face)

# example: generate a numpy array from the iterator (numpy 1.23 or later)
fca = np.fromiter(faces_iterator(ds), dtype=np.dtype((np.int64, NNODES)))
1 Like
5

Educational as always, thanks for your answer.