nschloe / meshio

:spider_web: input/output for many mesh formats
MIT License
1.93k stars 397 forks source link

[BUG] VTK CELL_DATA is not written in correct format #1457

Open reox opened 7 months ago

reox commented 7 months ago

Describe the bug meshio can properly read scalar, vector, and tensorial quantities from cell_data blocks. However, when these are written again as VTK, the structure is different and the files cannot be loaded in ParaView:

ERROR: In C:\glr\builds\paraview\paraview-ci\build\superbuild\paraview\src\VTK\IO\Legacy\vtkDataReader.cxx, line 838
vtkUnstructuredGridReader (000002505AAA51C0): Unsupported cell attribute type: 0.0 for file: S:\tmp\bugreport_meshio_rewritten.vtk

To Reproduce

Here is a small VTK file that contains all three types of Cell Data:

# vtk DataFile Version 2.0
Analysis
ASCII
DATASET UNSTRUCTURED_GRID
POINTS 4 float
-1.0 -1.0 0.0
1.0 -1.0 0.0
0.0 1.0 0.0
0.0 0.0 1.0

CELLS 1 5
4 0 1 2 3

CELL_TYPES 1
10

CELL_DATA 1
SCALARS SomeScalar float 1
LOOKUP_TABLE default
23.0

VECTORS SomeVector float
25.2 24.2 42.1

TENSORS SomeTensor float
2.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 4.0

Write the file again with meshio:

Python 3.9.7 (default, Sep 16 2021, 13:09:58) 
[GCC 7.5.0] :: Anaconda, Inc. on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> import meshio
>>> meshio.__version__
'5.3.5'
>>> mesh = meshio.read('bugreport_meshio.vtk')
>>> mesh.write('bugreport_meshio_rewritten.vtk', binary=False)
Warning: VTK ASCII files are only meant for debugging.

The new mesh looks like this:

# vtk DataFile Version 5.1
written by meshio v5.3.5
ASCII
DATASET UNSTRUCTURED_GRID
POINTS 4 float
-1.0 -1.0 0.0 1.0 -1.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0
CELLS 2 4
OFFSETS vtktypeint64
0
4
CONNECTIVITY vtktypeint64
0
1
2
3
CELL_TYPES 1
10
CELL_DATA 1
FIELD FieldData 3
SomeScalar 1 1 float
23.0
SomeVector 3 1 float
25.200000762939453 24.200000762939453 42.099998474121094
SomeTensor 3 1 float
2.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 4.0

As you can see, the dimensions for the tensor is completely wrong. It seems like, that even though the rest is also different (was that changed in newer VTK versions?) the only line that breaks paraview is actually the wrong dimension of the tensor. If I replace the 3 with 9, it loads.

The problem seems to come from here: https://github.com/nschloe/meshio/blob/b2ee99842e119901349fdeee06b5bf61e01f450a/src/meshio/vtk/_vtk_51.py#L640-L643

As the input data is of shape (N, 3, 3), only the second dimension is checked while the third is never looked at.

I guess the fix here would be to add another elif that checks for tensors and writes the product of the dimensions? However, I'm not familiar with VTK 5 format - all my scripts still write v2 (as you can see above :D)

Diagnose I may ask you to cut and paste the output of the following command.

pip freeze | grep meshio
meshio==5.3.5

Did I help?

If I was able to resolve your problem, consider sponsoring my work on meshio, or buy me a coffee to say thanks.