Connectivity mismatch between WriteVTK and ReadVTK

[Matthew-Whisenant]

I am currently writing a function that can take VTU data read from ReadVTK and use that data to write a new VTU in WriteVTK.

As a test, I am doing the simple write VTU example on WriteVTK here.

Using ReadVTK's get_cells, one can obtain the connectivity list and corresponding offset vector. However, I noticed that the connectivity list is starts at index 0 while WriteVTK starts at 1, resulting in mismatch of connectivity data. Here is a minimal example:

using ReadVTK
using WriteVTK

function convert_cell_structs(cells::VTKCells)

    start_offsets = [0; cells.offsets[1:end-1]] .+ 1
    end_offsets = cells.offsets

    offset_ranges = range.(start_offsets, end_offsets)

    connectivity = getindex.([cells.connectivity], offset_ranges)
    cell_types = VTKCellType.(cells.types)

    return MeshCell.(cell_types, connectivity)

end

points = rand(3, 5)  # 5 points in three dimensions
cells = [MeshCell(VTKCellTypes.VTK_TRIANGLE, [1, 4, 2]),
    MeshCell(VTKCellTypes.VTK_QUAD, [2, 4, 3, 5])]

vtk_grid("filename", points, cells) do vtk
    # add datasets...
end

vtkcase = VTKFile("filename.vtu")
points2 = get_points(vtkcase)
cells2 = convert_cell_structs(get_cells(vtkcase))

The connectivity input can be seen as

2-element Vector{MeshCell{VTKCellType, Vector{Int64}}}:
 MeshCell{VTKCellType, Vector{Int64}}(VTKCellType("VTK_TRIANGLE", 0x05, 3), [1, 4, 2])
 MeshCell{VTKCellType, Vector{Int64}}(VTKCellType("VTK_QUAD", 0x09, 4), [2, 4, 3, 5])

and for ReadVTK:

tempcells = get_cells(vtkcase)
VTKCells{Base.ReinterpretArray{Int64, 1, UInt8, Vector{UInt8}, false}, Base.ReinterpretArray{Int64, 1, UInt8, Vector{UInt8}, false}, Vector{UInt8}}([0, 3, 1, 1, 3, 2, 4], [3, 7], UInt8[0x05, 0x09])

using my function:

2-element Vector{MeshCell{VTKCellType, Vector{Int64}}}:
 MeshCell{VTKCellType, Vector{Int64}}(VTKCellType("VTK_TRIANGLE", 0x05, 3), [0, 3, 1])
 MeshCell{VTKCellType, Vector{Int64}}(VTKCellType("VTK_QUAD", 0x09, 4), [1, 3, 2, 4])

where you can see the connectivity should be +1 for every index. This is seen in the VTKCells structure even before applying my function. Of course, I could and currently am just adding 1 to it so it comes out normal.

I say all this to say, it seems strange since this is suppose to read what WriteVTK outputs, so the connectivity should match too. In fact, it would make sense to me to just output the MeshCell type that WriteVTK is using rather than a new VTKCells type. That way, reading and writing would be relatively low cost as the points from get_points and the cells from get_cells are the same as in the vtk_grid function input.

This would allow for use cases of: external software VTK -> ReadVTK ->Julia processing ->WriteVTK -> Visualize VTKs

Whereas the current use case seems to be: Julia processing ->WriteVTK -> ReadVTK -> more Julia processing?

I do not mind contributing the function that can convert VTKCells -> Vector{MeshCell}.

Sorry for the long post, just thought it may be worth thinking about.

[sloede]

Sorry for the long post, just thought it may be worth thinking about.

Not at all, thanks for contributing to this package!

Generally speaking, I think you are touching upon two (somewhat) orthogonal topics in this issue:

The first topic is the question of zero-based (C/VTK style) and one-based (Julia style) indexing. When storing the point connectivity array, WriteVTK.jl explicitly converts from one-based to zero-based indices: https://github.com/jipolanco/WriteVTK.jl/blob/4b4e406f2cbe5a0b376b7d8c50c56d59323e9772/src/gridtypes/unstructured/unstructured.jl#L118 Thus I think it would be a reasonable change to ReadVTK.jl to also do this conversion automatically. As far as I can tell, it would not be very difficult, as one would only have to modify the get_cells method after this line https://github.com/trixi-framework/ReadVTK.jl/blob/8ceeb4f9b556fee6066aa170074533dad9f50a92/src/ReadVTK.jl#L1036 to add oneunit to the connectivity and offset arrays.

The second topic is the ability to seamlessly process VTK data with WriteVTK.jl/ReadVTK.jl and use both packages together in an intuitive way. I completely agree: One should be able to read a VTK file with ReadVTK.jl and write it out again with WriteVTK.jl and get identical VTK files on disk, which is currently not possible. Thus, you have finally reached the ugly incompatibility that I was hoping nobody would notice too soon 😬 So yes, ReadVTK and WriteVTK should use the same types, but how this can be achieved in a sane manner is not yet clear to me.

I will reach out to the WriteVTK developers regarding the second point. In the meantime, for the first issue (one vs zero based indices), I would be happy to accept a PR where you modify the get_cells method such that it converts all indices to one-based indices. Would you be willing to give it a try @Matthew-Whisenant?

[sloede]

I've brought this up in https://github.com/jipolanco/WriteVTK.jl/issues/23#issuecomment-1490163304

[Matthew-Whisenant]

Yes. I have this update, the previous from yesterday, and another change that allows get_data to read non-appended binary OpenFOAM vtk files (couldn't figure out the error with appended data I was getting).

These along with my function that converts ReadVTK cell data to WriteVTK cell data since I have it anyway and may be useful for anyone that will use my current workflow.

Would these be good for a single PR or should I do multiple separate for clarity? Probably still multiple commits at least?

I just need to finish making a test case and it should be good.

[sloede]

It should definitely be a separate PR for each orthogonal feature. Why don't you go ahead and make a PR for changing the offsets/connectivity from zero-based to one-based indices now?

This is IMHO independent from the PR that converts the cell data to use actual VTKCellType cells, for which we have to wait at least 2 more days for VTKBase.jl to be registered as a package.