Cannot get faces from unstructured polyhedral mesh

[pcsanematsu]

I have an XML unstructured grid with polyhedral elements and would like to read it using Python for post-processing analysis of a simulation. I am using Python 3.7.4 and PyVista 0.22.4.

I can successfully read the mesh elements, points, and cell arrays doing the following:

import numpy as np  
import pyvista as pv

# open *.vtu file
mesh = pv.read("cells_0.0000.vtu")

myCells = mesh.cells
print(mesh.n_cells)
print(mesh.cells)

myPoints = mesh.points
print(mesh.n_points)
print(mesh.points)

print(mesh.n_arrays)

# read some arrays
cellEnergy = mesh.get_array('cellEnergy','cell')
print("cellEnergy")
print(cellEnergy)

cellType = mesh.get_array('cellType','cell')
print("cellType")
print(cellType)

I also need to read the faces of polyhedral elements. More specifically, I would like to read the faces and facesoffsets arrays within the Cell category (note: not within CellData). I noticed that myCells stores information from the connectivity and offsets arrays within Cell, but I couldn't find a function to get the faces information. Is there a function to retrieve faces and facesoffsets?

I attached an example file that I have been using to test the Python code above.

cells_0.0000.vtu.zip

category of the XML file)

[banesullivan]

The offsets can be gathered by:

mesh.offset

and the faces can be gathered like the following. There is probably a faster way to do this using the offset array with some numpy magic but this works:

faces = []
i, offset = 0, 0
cc = mesh.cells # fetch up front
while i < mesh.n_cells:
    nn = cc[offset]
    faces.append(cc[offset+1:offset+1+nn])
    offset += nn + 1
    i += 1

and so the faces list contains the point indices that make up every cell.

and just to prove that the above snippet works, we can iterate over every cell in the mesh and make sure the points listed in the faces array match each cell when extracting:

for idx in range(mesh.n_cells):
    pts = mesh.points[faces[idx]]
    cell = mesh.extract_cells(idx)
    assert np.allclose(pts, cell.points)

[banesullivan]

@akaszynski, should we add the snippet above to the UnstructuredGrid class to fetch the faces? Or is there a better way to do this?

[pcsanematsu]

Thanks for your reply. This information was very helpful. I can see the offsets and the point indices of each cell. I'm not sure if I'm missing something in your snippet, but I think there is still some information needed that is missing. Since I am working with polyhedral cells, I not only need the point indices of each cell, but I also need (1) the number of faces of each cell and (2) the number of points in each face. For example, this is the necessary information to describe a polyhedron cell (from https://vtk.org/Wiki/VTK/Polyhedron_Support):

[numberOfCellFaces, (numberOfPointsOfFace0, pointId0, pointId1, … ), (numberOfPointsOfFace1, pointId0, pointId1, …), … ]. 

In the original attachment (cells_0.0000.vtu.zip), the above information is located in the faces and facesoffsets arrays inside Cells, which I can't see that it was retrieved in the snipped.

[banesullivan]

Ah, my bad - I misinterpreted your question. We don't have any convenience methods for accessing the polyhedral's faces like that, but since PyVista is a direct layer on top of VTK, it's all still there via the VTK interface.

So if you want the array

[numberOfCellFaces, (numberOfPointsOfFace0, pointId0, pointId1, … ), (numberOfPointsOfFace1, pointId0, pointId1, …), … ].

then simply call:

>>> pv.convert_array(mesh.GetFaces())
array([   13,     7,     6, ..., 88857, 88856, 88855])

[banesullivan]

And here is a snippet to make a nested list of those faces:

stream = pv.convert_array(mesh.GetFaces())
locs = pv.convert_array(mesh.GetFaceLocations())

cells = []
for i in range(len(locs)):
    if i >= len(locs) - 1:
        stop = -1
    else:
        stop = locs[i+1]
    cells.append(stream[locs[i]:stop])
len(cells)

faces = []
for cell in cells:
    cell_faces = []
    i, offset = 0, 1
    nf = cell[0]
    while i < nf:
        nn = cell[offset]
        cell_faces.append(cell[offset+1:offset+1+nn])
        offset += nn + 1
        i += 1
    faces.append(cell_faces)

such that a pseudo 3D array is returned in faces - first level is for each cell, the second level is for all the faces in that cell, and the third level is the indices of points for that face

[banesullivan]

I'm thinking we may want to put these methods in pyvista's UnstructuredGrid class as helper routines

[pcsanematsu]

I didn't realize that GetFaces and GetFaceLocations get the arrays that I need!

Thank you for creating the snippet in addition to providing the vtk functions. The snipped does exactly what I need and puts the face data in a really nice structure for user-friendly access. Thanks again! Now I can do my data analysis :)

[banesullivan]

happy to help! i’m going to leave this open until we decide if/how we should add those helper methods in pyvista

_{Sent with GitHawk}

[GonzaloSaez]

Hi @banesullivan thank you for that code snippet to better manage faces in polyhedral meshes. Do you know if it is possible to obtain a similar face array for hexahedral or tetrahedra meshes? I now how access node connectivity but I don't know how to extract node to face or face to cell connectivity of non-polyhedral meshes with vtk

Thank you

[pcsanematsu]

stream = pv.convert_array(mesh.GetFaces())
locs = pv.convert_array(mesh.GetFaceLocations())

cells = []
for i in range(len(locs)):
    if i >= len(locs) - 1:
        stop = -1
    else:
        stop = locs[i+1]
    cells.append(stream[locs[i]:stop])
len(cells)

faces = []
for cell in cells:
    cell_faces = []
    i, offset = 0, 1
    nf = cell[0]
    while i < nf:
        nn = cell[offset]
        cell_faces.append(cell[offset+1:offset+1+nn])
        offset += nn + 1
        i += 1
    faces.append(cell_faces)

It's me again. I am working on a different dataset (attached) and noticed that the snippet above may not be storing the very last point of the *.vtu file.

When I print the last face of the last cell, I only get two vertices -- which should not happen since a face must have at least 3 vertices.

print(faces[2114][16][:])
[157259 157258]

Could you please help? I can't find where there might be a bug.

cells_1003.9400.vtu.zip

[ttsesm]

Ah, my bad - I misinterpreted your question. We don't have any convenience methods for accessing the polyhedral's faces like that, but since PyVista is a direct layer on top of VTK, it's all still there via the VTK interface.

So if you want the array

[numberOfCellFaces, (numberOfPointsOfFace0, pointId0, pointId1, … ), (numberOfPointsOfFace1, pointId0, pointId1, …), … ].

then simply call:

>>> pv.convert_array(mesh.GetFaces())
array([   13,     7,     6, ..., 88857, 88856, 88855])

@banesullivan I want to get the faces as well from a PolyData are these helper methods GetFaces() and GetFaceLocations() available somehow because I do not find them anywhere.

[banesullivan]

@ttsesm, can you share a pv.Report()? I suspect you may need to upgrade VTK

@paulasoo, sorry for not responding... I will try to look into that bug

[banesullivan]

@paulasoo, fixed the bug!

    if i >= len(locs) - 1:
        stop = -1

needed to be

    if i >= len(locs) - 1:
        stop = None

so the snippet is:

stream = pv.convert_array(mesh.GetFaces())
locs = pv.convert_array(mesh.GetFaceLocations())

cells = []
for i in range(len(locs)):
    if i >= len(locs) - 1:
        stop = None
    else:
        stop = locs[i+1]
    cells.append(stream[locs[i]:stop])
len(cells)

faces = []
for cell in cells:
    cell_faces = []
    i, offset = 0, 1
    nf = cell[0]
    while i < nf:
        nn = cell[offset]
        cell_faces.append(cell[offset+1:offset+1+nn])
        offset += nn + 1
        i += 1
    faces.append(cell_faces)

[ttsesm]

@banesullivan I am quite sure that have the latest vtk, and here is the report that confirms it:

--------------------------------------------------------------------------------
  Date: Fri Jun 05 10:00:19 2020 W. Europe Daylight Time

                OS : Windows
            CPU(s) : 12
           Machine : AMD64
      Architecture : 64bit
               RAM : 31.7 GB
       Environment : Python
        GPU Vendor : Intel
      GPU Renderer : Intel(R) UHD Graphics 630
       GPU Version : 4.5.0 - Build 26.20.100.7985

  Python 3.8.3 (tags/v3.8.3:6f8c832, May 13 2020, 22:37:02) [MSC v.1924 64 bit
  (AMD64)]

           pyvista : 0.24.2
               vtk : 9.0.0
             numpy : 1.18.4
           imageio : 2.8.0
           appdirs : 1.4.4
            scooby : 0.5.4
        matplotlib : 3.2.1
             PyQt5 : 5.14.2
           IPython : 7.14.0
             scipy : 1.4.1
--------------------------------------------------------------------------------

But if I give the following snippet:

sphere = pv.Sphere(radius=0.85)
faces = sphere.GetFaces()

I am getting an AttributeError: 'PolyData' object has no attribute 'GetFaces'

[banesullivan]

Ah there’s your problem. This is for unstructured meshes only and specifically for polyhedral cells (not all cells). PolyData can’t have polyhedral cells, so this won’t work for your mesh.

However if you want the faces array, that’s already available for your mesh under the faces property

[ttsesm]

I see, yes I know that there is the faces property under the PolyData() class but it returns them in a different format from what I would expect (e.g. nx3). Thus, I need to reshape them and also remove the connectivity indicator from the beginning. So I wrote the following help function:

    def get_faces(polydata_faces, num_of_vertices=3):
        return polydata_faces.reshape(-1,num_of_vertices+1)[:,1:num_of_vertices+1]

Anyways, it is not a big hassle.

[banesullivan]

class but it returns them in a different format from what I would expect (e.g. nx3). Thus, I need to reshape them and also remove the connectivity indicator from the beginning. So I wrote the following help function:

Yep, that's just the way it is. PolyData supports quad faces in addition to triangles so we can't return an nx3 array as quads have four vertices. Mixing face types is supported and a very common thing, so there is no guarantee your helper method will always work - if you have a quad cell in your mesh, you will run into sizing issues.

Here is a safer helper method:

def get_tri_faces(poly):
    """Fetch all triangle faces."""
    stream = poly.faces
    tris = []
    i = 0
    while i < len(stream):
        n = stream[i]
        if n != 3:
            i += n + 1
            continue
        stop = i + n + 1
        tris.append(stream[i+1:stop])
        i = stop
    return np.array(tris)

[pcsanematsu]

@paulasoo, fixed the bug!

    if i >= len(locs) - 1:
        stop = -1

needed to be

    if i >= len(locs) - 1:
        stop = None

so the snippet is:

stream = pv.convert_array(mesh.GetFaces())
locs = pv.convert_array(mesh.GetFaceLocations())

cells = []
for i in range(len(locs)):
    if i >= len(locs) - 1:
        stop = None
    else:
        stop = locs[i+1]
    cells.append(stream[locs[i]:stop])
len(cells)

faces = []
for cell in cells:
    cell_faces = []
    i, offset = 0, 1
    nf = cell[0]
    while i < nf:
        nn = cell[offset]
        cell_faces.append(cell[offset+1:offset+1+nn])
        offset += nn + 1
        i += 1
    faces.append(cell_faces)

It's been a while, but thanks for fixing the bug! I will update my codes :)

[gmhaber]

I didn't read in detail the above discussion but if the purpose is to obtain the cells from an unstructured grid vtu file , check: np.array(pv.read(file_name).cells).reshape((-1, n_vertices+1))[:, 1:]

Much faster than looping through every cell id using a for loop and doesn't rely on vtkCommonDataModelPython.vtkCellArray that some distributions of pyvista or vtk lack.