Hi there!
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 coordinate array:
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 coordinate:
print(connectivity.max())
139903
print(len(connectivity[connectivity>coordinate.shape[0]]))
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?