JuliaVTK / WriteVTK.jl

Julia package for writing VTK XML files
Other
151 stars 32 forks source link

Missing XML attribute for `RationalWeights` and `HigherOrderDegrees` in Bézier cells #131

Closed m1ka05 closed 8 months ago

m1ka05 commented 8 months ago

Hi, I've been experimenting with Bézier cells and found that the rational weights and higher order degrees are not correctly exported. At least Paraview v5.11.2 using VTK v9.2.20220823 does not pick them up.

I attach two screenshots of the tetra quartic solid sphere octant test from bezier.jl with and without correct attribute definition. The following fixes the test:

    outfiles = vtk_grid(VTK_BASENAME, points, cells; vtkversion = v"1.0") do vtk
        vtk["RationalWeights"] = rational_weights

        # fix
        grid = WriteVTK.find_element(WriteVTK.root(vtk.xdoc), vtk.grid_type)
        piece = WriteVTK.find_element(grid, "Piece")
        pointdata = WriteVTK.find_element(piece, "PointData")
        WriteVTK.set_attribute(pointdata, "RationalWeights", "RationalWeights")
    end

Here is also a single cell quarter annulus PoC with anisotropic degrees:

using WriteVTK

cell_type = VTKCellTypes.VTK_BEZIER_QUADRILATERAL
points = [1.0 0.0 0.0 2.0 1.0 2.0; 0.0 1.0 2.0 0.0 1.0 2.0]
connectivity = 1:size(points, 2)
rational_weights = [1.0, 1.0, 1.0, 1.0, sqrt(0.5), sqrt(0.5)]
cells = [MeshCell(cell_type, connectivity)]

outfiles = vtk_grid("bezier_test", points, cells; vtkversion=v"1.0", compress=false, append=false, ascii=true) do vtk
    vtk["HigherOrderDegrees", VTKCellData()] = [2.0 1.0 1.0]
    vtk["RationalWeights", VTKPointData()] = rational_weights

    grid = WriteVTK.find_element(WriteVTK.root(vtk.xdoc), vtk.grid_type)
    piece = WriteVTK.find_element(grid, "Piece")
    celldata = WriteVTK.find_element(piece, "CellData")
    pointdata = WriteVTK.find_element(piece, "PointData")
    WriteVTK.set_attribute(celldata, "HigherOrderDegrees", "HigherOrderDegrees")
    WriteVTK.set_attribute(pointdata, "RationalWeights", "RationalWeights")
end

Fixed tetra quartic solid sphere octant

fixed

Broken tetra quartic solid sphere octant

broken

Quarter annulus

qa

jipolanco commented 8 months ago

Hi, thank you for catching this! I will fix the test, and maybe add your second example as an additional test (unless you want to do it yourself in a PR!). It may also be worth it to provide a simpler way of setting these kinds of attributes.

There's a slightly simpler way of setting an attribute to "PointData" and "CellData" elements. It is not very well documented, but it's mentioned at the very end of this section in the docs (and also here). Basically, your fix can be written as:

    outfiles = vtk_grid(VTK_BASENAME, points, cells; vtkversion = v"1.0") do vtk
        vtk["RationalWeights"] = rational_weights

        # fix
        vtk[VTKPointData()] = "RationalWeights" => "RationalWeights"
    end
m1ka05 commented 8 months ago

Thanks, I've missed that line in the documentation, great! I'll open a PR with the tests later today so you can link it with this issue.