How to generate a correct VTS grid using WriteVTK.jl?

Ahmed_Salih · October 30, 2024, 10:07pm

Hello!

I have the following code:

using WriteVTK
# Cell dimensions
dx = 1.0
dy = 1.0

# Initialize lists for storing points and connectivity
points = Vector{SVector{3, Float64}}()    # List to store unique SVector points
cells = Int32[]                           # List for connectivity of each cell
offsets = Int32[]                         # List for offset positions
cell_types = Int8[]                       # VTK cell type array

# Example set of CartesianIndex cells (5x5 grid)
UniqueCells = [CartesianIndex(i, j) for i in 1:5, j in 1:2]

# Loop through each CartesianIndex cell
iter = 
for cell in UniqueCells
    iter += 1
    # Get x and y from the CartesianIndex and calculate cell center
    xi, yi = cell.I
    x_center = (xi - 0.5) * dx
    y_center = (yi - 0.5) * dy

    # Define corners in counterclockwise order
    x0, y0 = x_center - dx / 2, y_center - dy / 2
    x1, y1 = x_center + dx / 2, y_center - dy / 2
    x2, y2 = x_center + dx / 2, y_center + dy / 2
    x3, y3 = x_center - dx / 2, y_center + dy / 2

    # Add corners as individual points
    push!(points, SVector(x0, y0, 0.0)); push!(cells, length(points))
    push!(points, SVector(x1, y1, 0.0)); push!(cells, length(points))
    push!(points, SVector(x2, y2, 0.0)); push!(cells, length(points))
    push!(points, SVector(x3, y3, 0.0)); push!(cells, length(points))
    
    # Ensure correct connectivity for each cell
    push!(offsets, length(cells))       # Track the end of the cell in the cells array
    push!(cell_types, 9)                  # VTK_QUAD type for square cells
end

# Write to VTK as unstructured grid
vtk_grid("ExtractedCells", points) do vtk
    vtk["cells"] = cells
    vtk["offsets"] = offsets
    vtk["celltypes"] = cell_types
end

Which generates:

This is of course wrong, since I want a grid of squares.

I want to use this and no other way, since my real use case is where UniqueCells are a lot of disjointed indices and this seems to work well, I am just struggling with the connectivity… maybe someone else has solved this before?

Kind regards

jipolanco · October 31, 2024, 8:41am

Hi, if I understand correctly what you’re trying to do, then the VTS file format is not the right format, as it represents structured grids whose cell types are fixed by the format (as described here).

If you want to choose the cell type, then you need to look at the unstructured formats, and in particular to the more general VTU format.

I’ve modified your example in order to create a VTU file with square cells:

using WriteVTK
using StaticArrays

# Cell dimensions
dx = 1.0
dy = 1.0

# Initialise lists for storing points and cells
points = Vector{SVector{3, Float64}}()    # List to store unique SVector points
cells = empty!([MeshCell(VTKCellTypes.VTK_QUAD, 1:4)])  # empty vector of cells (with the right element type)
cells_id = Int[]  # not needed, but useful for visualisation purposes

# Example set of CartesianIndex cells (5x5 grid)
UniqueCells = [CartesianIndex(i, j) for i in 1:5, j in 1:2]

# Loop through each CartesianIndex cell
for (id, cell) in enumerate(UniqueCells)
    # Get x and y from the CartesianIndex and calculate cell center
    xi, yi = cell.I
    x_center = (xi - 0.5) * dx
    y_center = (yi - 0.5) * dy

    # Define corners in counterclockwise order
    x0, y0 = x_center - dx / 2, y_center - dy / 2
    x1, y1 = x_center + dx / 2, y_center - dy / 2
    x2, y2 = x_center + dx / 2, y_center + dy / 2
    x3, y3 = x_center - dx / 2, y_center + dy / 2

    # Create quad cell connecting the 4 new points
    n = length(points)
    cell = MeshCell(VTKCellTypes.VTK_QUAD, (n + 1):(n + 4))
    push!(cells, cell)
    push!(cells_id, id)

    # Add corners as individual points
    push!(points, SVector(x0, y0, 0.0))
    push!(points, SVector(x1, y1, 0.0))
    push!(points, SVector(x2, y2, 0.0))
    push!(points, SVector(x3, y3, 0.0))
end

# Write to VTK as unstructured grid
vtk_grid("ExtractedCells", points, cells) do vtk
    vtk["cell_id"] = cells_id
end

Result:

Ahmed_Salih · October 31, 2024, 9:11pm

Thank you so much! This was exactly what I needed to be able to visualize the grid I use for my simulation tool (SPHExample.jl) and now I can also do it every time step to do motion!

I need to do it using VTKHDF in the future for performance, but I think what you gave me is a very good head start, so thanks a lot!

Kind regards

Ahmed_Salih · November 1, 2024, 10:58pm

Here is how to do it through vtkhdf:

    function SaveCellGridVTKHDF(FilePath, SimConstants, UniqueCells)
        # Cell dimensions
        dx = SimConstants.H
        dy = SimConstants.H
    
        # Initialize lists for storing points and cells
        points = Vector{SVector{3, Float64}}()  # List to store unique SVector points
        connectivity = Int[]                    # Connectivity for each cell
        offsets = Int[]                         # Offsets for each cell
        cell_types = Int[]                      # Cell types (for VTK_QUAD)
        cell_data  = Int[]
    
        push!(offsets, 0)
        # Loop through each CartesianIndex cell
        for (id, cell) in enumerate(UniqueCells)
            # Get x and y from the CartesianIndex and calculate cell center
            xi, yi = cell.I
            x_center = (xi - 0.5) * dx
            y_center = (yi - 0.5) * dy
    
            # Define corners individually
            x0, y0 = x_center - dx / 2, y_center - dy / 2
            x1, y1 = x_center + dx / 2, y_center - dy / 2
            x2, y2 = x_center + dx / 2, y_center + dy / 2
            x3, y3 = x_center - dx / 2, y_center + dy / 2
    
            # Add each corner point individually and update connectivity
            n = length(points)
            push!(points, SVector(x0, y0, 0.0))
            push!(connectivity, n)
            n += 1
    
            push!(points, SVector(x1, y1, 0.0))
            push!(connectivity, n)
            n += 1
    
            push!(points, SVector(x2, y2, 0.0))
            push!(connectivity, n)
            n += 1
    
            push!(points, SVector(x3, y3, 0.0))
            push!(connectivity, n)
    
            # Define cell type and offsets
            push!(offsets, length(connectivity))
            push!(cell_types, VTKCellTypes.VTK_QUAD.vtk_id)  # Convert to Int

            push!(cell_data, id)
        end
    
        # Open HDF5 file for writing
        io = h5open(FilePath, "w")

        # Create top-level group "VTKHDF"
        gtop = HDF5.create_group(io, "VTKHDF")

        # Set the Version attribute
        HDF5.attrs(gtop)["Version"] = [2, 1]

        # Write Type attribute as ASCII string
        let s = "UnstructuredGrid"
            dtype = HDF5.datatype(s)
            HDF5.API.h5t_set_cset(dtype.id, HDF5.API.H5T_CSET_ASCII)
            dspace = HDF5.dataspace(s)
            attr = HDF5.create_attribute(gtop, "Type", dtype, dspace)
            HDF5.write_attribute(attr, dtype, s)
        end

        # Write Number of Points, Number of Cells, and Number of Connectivity IDs
        gtop["NumberOfPoints"]          = [length(points)]
        gtop["NumberOfCells"]           = [length(cell_types)]
        gtop["NumberOfConnectivityIds"] = [length(connectivity)]

        # Write Points
        gtop["Points"] = reinterpret(reshape, eltype(eltype(points)), points)

        # Write Connectivity, Offsets, and Types
        gtop["Connectivity"] = connectivity
        gtop["Offsets"] = offsets
        gtop["Types"] = [VTKCellTypes.VTK_QUAD.vtk_id for _ in cell_types]  # Use VTKCellTypes Quad vtk_id for all types

        # Write CellData (cell-level variables)
        let cell_group = HDF5.create_group(gtop, "CellData")
            cell_group["CellData"] = cell_data
            close(cell_group)
        end

        # Write PointData (point-level variables)
        # let point_group = HDF5.create_group(gtop, "PointData")
        #     for (name, data) in point_data
        #         point_group[name] = data
        #     end
        #     close(point_group)
        # end

        # Write an empty FieldData group (placeholder for additional data)
        create_group(gtop, "FieldData")
        
        # Close file
        close(io)
    end

Performance difference between WriteVTK and using hdf is:

12B-WriteVTK Output CellGrid         200    7.23s    6.2%  36.1ms    451MiB   32.4%  2.25MiB
12B-VTKHDF   Output CellGrid         200    4.77s    4.1%  23.9ms    449MiB   32.2%  2.24MiB

Code to test:

github.com

AhmedSalih3d/SPHExample/blob/1eee2732f1b1fe617549cf9c83eb93a603c88dc4/example/Dambreak2d.jl

using SPHExample

let
    Dimensions = 2
    FloatType  = Float64

    SimConstantsDambreak = SimulationConstants{FloatType}(dx=0.02,c₀=88.14487860902641, δᵩ = 0.1, CFL=0.2, α = 0.02)

    # Create Geometry instances
    FixedBoundary = Geometry{Dimensions, FloatType}(
        CSVFile     = "./input/dam_break_2d/DamBreak2d_Dp0.02_Bound_ThreeLayers.csv",
        GroupMarker = 1,
        Type        = Fixed,   # Using the enum value Fixed
        Motion      = nothing
    )

    Water = Geometry{Dimensions, FloatType}(
        CSVFile     = "./input/dam_break_2d/DamBreak2d_Dp0.02_Fluid_OneLayer.csv",
        GroupMarker = 2,
        Type        = Fluid,   # Using the enum value Fluid

This file has been truncated. show original

(Requires installing the package from this specific branch and running this specific file)

jipolanco · November 2, 2024, 12:48pm

Just curious, did you disable compression when trying out WriteVTK.jl? By default the datasets are written in compressed form^[1] which can significantly slow down writing files. You can disable compression by passing compress = false to vtk_grid.

I admit this is not necessarily a good default. The decision was made many years ago, but it can be changed. ↩︎

Ahmed_Salih · November 3, 2024, 10:00am

No, I unfortunately did not get to test that

I think tbh that it is an okay default, since in reality most of the time you want the smallest file as possible.

In my case the reason I don’t use WriteVTK is that I need to close my files at the end only for performance reasons (i.e. all threads stop working if I start saving files and I have not found a good fix).

I also think that VTKHDF will be the future of the format, so I am trying to leave behind the xml versions, since being able to load in VTKHDF directly using a HDFJ5.jl library is so strong, with no special wrappers

Kind regards