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:
