I’m working on a Python project that requires reading parallel output from a finite element model on an unstructured hexahedral grid. The output files are .pvtu (one per timestep) and .vtu (one per timestep per processor) with different segments of the unstructured mesh split across .vtu files. The goal of this project is to add VTK support to an existing 3D plotting library which has a general unstructured mesh renderer, so I just need to load the files and pass along the connectivity, coordinate and data arrays appropriately. I started off using the python
meshio package to load individual .vtu files listed in the .pvtu, which I then collect and pass to my mesh renderer. There are lots of places to mess up here, but I’m struggling with something that I think is a basic question:
I’ve found that the connectivity array within one of my .vtu files contains elements that are not within that .vtu file : is this possible? Is it an error? How do I deal with it?
I think this will be clearer with an example:
After parsing the “problematic” .vtu file, the
connectivity array has 17480 elements, 8 vertices per element:
print(connectivity.shape) (17480, 8) print(connectivity.size) 139840
The values of
connectivity should point to indices of the
print(coordinate.shape) (139840, 3)
The first dimension of the
coordinate array matches the length of the
connectivity – makes sense so far!
But there are values in the
connectivity array that exceed the dimensions of
print(connectivity.max()) 139903 print(len(connectivity[connectivity>coordinate.shape])) 63
So that evaluating the coordinate values, for example
coordinate[connectivity.ravel(),0], will throw an error.
So what are these indices pointing to?
GhostLevel=0 for these files, so I didn’t expect overlap between
.vtu files, but even if these indices refer to coordinates in other .vtu files, how do I know which ones? The connectivity values and coordinate arrays reset between .vtu files.
I’m fairly certain the .pvtu and .vtu files are correct – they load and plot as expected in paraview. So I think I’m missing some piece of understanding here. Any ideas?