From 13c550b33cc5071b03d5d6e8ed4bbdb46e86e6d6 Mon Sep 17 00:00:00 2001
From: Juan Ignacio Polanco <jipolanc@gmail.com>
Date: Tue, 18 Jan 2022 08:52:54 +0100
Subject: [PATCH 01/18] Rework extents and image data writers

---
 src/gridtypes/common.jl                 | 47 +++++++++++++-------
 src/gridtypes/structured/imagedata.jl   | 57 +++++++++++--------------
 src/gridtypes/structured/rectilinear.jl |  4 +-
 src/gridtypes/structured/structured.jl  |  4 +-
 test/extent.jl                          | 27 ++++++++++++
 test/imagedata.jl                       |  8 ++--
 test/runtests.jl                        |  2 +
 7 files changed, 93 insertions(+), 56 deletions(-)
 create mode 100644 test/extent.jl

diff --git a/src/gridtypes/common.jl b/src/gridtypes/common.jl
index ad02e315..49c1adab 100644
--- a/src/gridtypes/common.jl
+++ b/src/gridtypes/common.jl
@@ -14,31 +14,48 @@ function vtk_xml_write_header(vtk::DatasetFile)
     xroot::XMLElement
 end
 
+# Specifies the "extent" of a structured grid, e.g. (1:4, 0:3, 1:42)
+const Extent{N} = Tuple{Vararg{AbstractUnitRange{<:Integer}, N}} where {N}
 
 # Returns the "extent" attribute required for structured (including rectilinear)
 # grids.
-extent_attribute(Ni, Nj, Nk, ::Nothing=nothing) =
-    "0 $(Ni - 1) 0 $(Nj - 1) 0 $(Nk - 1)"
+function extent_attribute(Ns::Dims{3}, extent::Extent{3})
+    all(pair -> pair[1] == length(pair[2]), zip(Ns, extent)) ||
+        throw(DimensionMismatch("extent is not consistent with dataset dimensions."))
+    iter = Iterators.map(e -> string(first(e), " ", last(e)), extent)
+    join(iter, " ")
+end
+
+function extent_attribute(Ns::Dims{N}, extent::Extent{N}) where {N}
+    if N > 3
+        throw(ArgumentError("dimensionalities N > 3 not supported (got N = $N)"))
+    elseif N < 3
+        M = 3 - N
+        extent_attribute(
+            (Ns..., ntuple(d -> 1, Val(M))...),
+            (extent..., ntuple(d -> 0:0, Val(M))...),
+        )
+    end
+end
 
-function extent_attribute(Ni, Nj, Nk, extent::Array{T}) where T <: Integer
-    length(extent) == 6 || throw(ArgumentError("extent must have length 6."))
-    (extent[2] - extent[1] + 1 == Ni) &&
-    (extent[4] - extent[3] + 1 == Nj) &&
-    (extent[6] - extent[5] + 1 == Nk) ||
-    throw(ArgumentError("extent is not consistent with dataset dimensions."))
-    join(extent, " ")
+function extent_attribute(Ns::Dims, ::Nothing = nothing)
+    ext = map(N -> 0:(N - 1), Ns)
+    extent_attribute(Ns, ext)
+end
+
+# This is left for compatibility...
+function extent_attribute(Ns::Dims{N}, ext::Array{T}) where {N, T <: Integer}
+    length(ext) == 2N || throw(ArgumentError("extent must have length $(2N)."))
+    rs = ntuple(i -> (ext[2i - 1]):(ext[2i]), Val(N))
+    extent_attribute(Ns, rs)
 end
 
 # Number of cells in structured grids.
 # In 3D, all cells are hexahedrons (i.e. VTK_HEXAHEDRON), and the number of
 # cells is (Ni-1)*(Nj-1)*(Nk-1). In 2D, they are quadrilaterals (VTK_QUAD), and in
 # 1D they are line segments (VTK_LINE).
-function num_cells_structured(Ni, Nj, Nk)
-    Ncls = one(Ni)
-    for N in (Ni, Nj, Nk)
-        Ncls *= max(1, N - 1)
-    end
-    Ncls
+function num_cells_structured(Ns::Dims)
+    prod(N -> max(1, N - 1), Ns)
 end
 
 # Accept passing coordinates as a tuple.
diff --git a/src/gridtypes/structured/imagedata.jl b/src/gridtypes/structured/imagedata.jl
index b853662a..afe6fdd8 100644
--- a/src/gridtypes/structured/imagedata.jl
+++ b/src/gridtypes/structured/imagedata.jl
@@ -1,11 +1,12 @@
-function vtk_grid(dtype::VTKImageData, filename::AbstractString,
-                  Nx::Integer, Ny::Integer, Nz::Integer=1;
-                  origin = (0.0, 0.0, 0.0),
-                  spacing = (1.0, 1.0, 1.0),
-                  extent = nothing, kwargs...)
-    Npts = Nx*Ny*Nz
-    Ncls = num_cells_structured(Nx, Ny, Nz)
-    ext = extent_attribute(Nx, Ny, Nz, extent)
+function vtk_grid(
+        dtype::VTKImageData, filename::AbstractString, Ns::Dims{N};
+        origin::NTuple{N} = ntuple(d -> 0.0, Val(N)),
+        spacing::NTuple{N} = ntuple(d -> 1.0, Val(N)),
+        extent = nothing, kwargs...,
+    ) where {N}
+    Npts = prod(Ns)
+    Ncls = num_cells_structured(Ns)
+    ext = extent_attribute(Ns, extent)
 
     xvtk = XMLDocument()
     vtk = DatasetFile(dtype, xvtk, filename, Npts, Ncls; kwargs...)
@@ -17,28 +18,8 @@ function vtk_grid(dtype::VTKImageData, filename::AbstractString,
     xGrid = new_child(xroot, vtk.grid_type)
     set_attribute(xGrid, "WholeExtent", ext)
 
-    No = length(origin)
-    Ns = length(spacing)
-
-    if No > 3
-        throw(ArgumentError("origin array must have length <= 3"))
-    elseif Ns > 3
-        throw(ArgumentError("spacing array must have length <= 3"))
-    end
-
-    origin_str = join(origin, " ")
-    spacing_str = join(spacing, " ")
-
-    while No != 3
-        # Fill additional dimensions (e.g. the z dimension if 2D grid)
-        origin_str *= (" 0.0")
-        No += 1
-    end
-
-    while Ns != 3
-        spacing_str *= (" 1.0")
-        Ns += 1
-    end
+    origin_str = _tuple_to_str3(origin, zero(eltype(origin)))
+    spacing_str = _tuple_to_str3(spacing, one(eltype(origin)))
 
     set_attribute(xGrid, "Origin", origin_str)
     set_attribute(xGrid, "Spacing", spacing_str)
@@ -50,8 +31,18 @@ function vtk_grid(dtype::VTKImageData, filename::AbstractString,
     vtk
 end
 
-vtk_grid(filename::AbstractString, xyz::Vararg{Integer}; kwargs...) =
-    vtk_grid(VTKImageData(), filename, xyz...; kwargs...)
+_tuple_to_str3(ts::NTuple{3, T}, default::T) where {T} = join(ts, " ")
+
+function _tuple_to_str3(ts::NTuple{N, T}, default::T) where {N, T}
+    @assert N < 3
+    M = 3 - N
+    tt = (ts..., ntuple(d -> default, Val(M))...)
+    _tuple_to_str3(tt, default)
+end
+
+function vtk_grid(filename::AbstractString, Ns::Vararg{Integer}; kwargs...)
+    vtk_grid(filename, map(N -> 1:N, Ns)...; kwargs...)
+end
 
 """
     vtk_grid(filename, x::AbstractRange{T}, y::AbstractRange{T}, [z::AbstractRange{T}];
@@ -86,6 +77,6 @@ function vtk_grid(filename::AbstractString, xyz::Vararg{AbstractRange}; kwargs..
     Nxyz = length.(xyz)
     origin = first.(xyz)
     spacing = step.(xyz)
-    vtk_grid(VTKImageData(), filename, Nxyz...;
+    vtk_grid(VTKImageData(), filename, Nxyz;
              origin=origin, spacing=spacing, kwargs...)
 end
diff --git a/src/gridtypes/structured/rectilinear.jl b/src/gridtypes/structured/rectilinear.jl
index e26b818c..44fc5e2f 100644
--- a/src/gridtypes/structured/rectilinear.jl
+++ b/src/gridtypes/structured/rectilinear.jl
@@ -3,8 +3,8 @@ function vtk_grid(dtype::VTKRectilinearGrid, filename::AbstractString,
                   extent=nothing, kwargs...)
     Ni, Nj, Nk = length(x), length(y), length(z)
     Npts = Ni*Nj*Nk
-    Ncls = num_cells_structured(Ni, Nj, Nk)
-    ext = extent_attribute(Ni, Nj, Nk, extent)
+    Ncls = num_cells_structured((Ni, Nj, Nk))
+    ext = extent_attribute((Ni, Nj, Nk), extent)
 
     xvtk = XMLDocument()
     vtk = DatasetFile(dtype, xvtk, filename, Npts, Ncls; kwargs...)
diff --git a/src/gridtypes/structured/structured.jl b/src/gridtypes/structured/structured.jl
index d34d4bfa..1f69bef6 100644
--- a/src/gridtypes/structured/structured.jl
+++ b/src/gridtypes/structured/structured.jl
@@ -15,8 +15,8 @@ function vtk_grid(dtype::VTKStructuredGrid, filename::AbstractString,
     Ni, Nj, Nk = structured_dims(xyz)
     Npts = Ni * Nj * Nk
     Ncomp = num_components(xyz, Npts)
-    Ncls = num_cells_structured(Ni, Nj, Nk)
-    ext = extent_attribute(Ni, Nj, Nk, extent)
+    Ncls = num_cells_structured((Ni, Nj, Nk))
+    ext = extent_attribute((Ni, Nj, Nk), extent)
 
     if Ncomp != 3  # three components (x, y, z)
         msg = "coordinate array `xyz` has incorrect dimensions.\n" *
diff --git a/test/extent.jl b/test/extent.jl
new file mode 100644
index 00000000..033c3d22
--- /dev/null
+++ b/test/extent.jl
@@ -0,0 +1,27 @@
+using Test
+using WriteVTK
+
+@testset "Extent" begin
+    @testset "N = 2" begin
+        Ns = (6, 8)
+        ext = (1:6, 3:10)
+        @test_throws DimensionMismatch WriteVTK.extent_attribute(Ns .- 1, ext)
+        @test WriteVTK.extent_attribute(Ns, ext) == "1 6 3 10 0 0"
+        @test WriteVTK.extent_attribute(Ns) == "0 5 0 7 0 0"
+    end
+
+    @testset "N = 3" begin
+        Ns = (6, 8, 42)
+        ext = (0:5, 0:7, 0:2)
+        @test_throws DimensionMismatch WriteVTK.extent_attribute(Ns, ext)
+        ext = (0:5, 0:7, 1:42)
+        @test WriteVTK.extent_attribute(Ns, ext) == "0 5 0 7 1 42"
+        @test WriteVTK.extent_attribute(Ns) == "0 5 0 7 0 41"
+    end
+
+    @testset "N = 4" begin
+        Ns = (3, 4, 5, 6)
+        ext = Base.OneTo.(Ns)
+        @test_throws ArgumentError WriteVTK.extent_attribute(Ns, ext)
+    end
+end
diff --git a/test/imagedata.jl b/test/imagedata.jl
index 17d709d2..f8120973 100644
--- a/test/imagedata.jl
+++ b/test/imagedata.jl
@@ -33,9 +33,9 @@ function main()
     end
 
     # These are all optional:
-    extent = [1, Ni, 1, Nj, 1, Nk] .+ 42
-    origin = [1.2, 4.3, -3.1]
-    spacing = [0.1, 0.5, 1.2]
+    extent = map(r -> r .+ 42, (1:Ni, 1:Nj, 1:Nk))
+    origin = (1.2, 4.3, -3.1)
+    spacing = (0.1, 0.5, 1.2)
 
     # Initialise new vti file (image data).
     @time outfiles = vtk_grid(vtk_filename_noext, Ni, Nj, Nk, extent=extent,
@@ -59,7 +59,7 @@ function main()
 
     # Test 2D dataset
     @time outfiles_2D = vtk_grid(vtk_filename_noext * "_2D", Ni, Nj,
-                                 # extent=extent[1:2],  # doesn't work for now (TODO)
+                                 # extent=extent[1:2],
                                  origin=origin[1:2], spacing=spacing[1:2]) do vtk
         vtk["p_values"] = p[:, :, 1]
         vtk["myVector"] = vec[:, :, :, 1]
diff --git a/test/runtests.jl b/test/runtests.jl
index 10314dc6..ed0c6fe9 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -4,6 +4,8 @@ using WriteVTK
 using Test: @test
 using SHA: sha1
 
+include("extent.jl")
+
 const tests = [
     "polyhedron_cube.jl",
     "multiblock.jl",

From f7c2bbfd41e86f07ae77f6a647224acfbba8d14c Mon Sep 17 00:00:00 2001
From: Juan Ignacio Polanco <jipolanc@gmail.com>
Date: Tue, 18 Jan 2022 09:05:42 +0100
Subject: [PATCH 02/18] Fix tests for vtk_write_array

---
 src/gridtypes/structured/imagedata.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/gridtypes/structured/imagedata.jl b/src/gridtypes/structured/imagedata.jl
index afe6fdd8..a322cbaa 100644
--- a/src/gridtypes/structured/imagedata.jl
+++ b/src/gridtypes/structured/imagedata.jl
@@ -41,7 +41,8 @@ function _tuple_to_str3(ts::NTuple{N, T}, default::T) where {N, T}
 end
 
 function vtk_grid(filename::AbstractString, Ns::Vararg{Integer}; kwargs...)
-    vtk_grid(filename, map(N -> 1:N, Ns)...; kwargs...)
+    # We put the origin at (0.0, 0.0, 0.0)
+    vtk_grid(filename, map(N -> range(0.0; length = N), Ns)...; kwargs...)
 end
 
 """

From a2fbafadc937f266afb2a1c5070db495ab2cf9ad Mon Sep 17 00:00:00 2001
From: Juan Ignacio Polanco <jipolanc@gmail.com>
Date: Tue, 18 Jan 2022 09:14:21 +0100
Subject: [PATCH 03/18] Refactor rectilinear grid writer

---
 src/gridtypes/structured/rectilinear.jl | 33 ++++++++++++++-----------
 1 file changed, 18 insertions(+), 15 deletions(-)

diff --git a/src/gridtypes/structured/rectilinear.jl b/src/gridtypes/structured/rectilinear.jl
index 44fc5e2f..05b66640 100644
--- a/src/gridtypes/structured/rectilinear.jl
+++ b/src/gridtypes/structured/rectilinear.jl
@@ -1,10 +1,12 @@
-function vtk_grid(dtype::VTKRectilinearGrid, filename::AbstractString,
-                  x::AbstractVector, y::AbstractVector, z::AbstractVector;
-                  extent=nothing, kwargs...)
-    Ni, Nj, Nk = length(x), length(y), length(z)
-    Npts = Ni*Nj*Nk
-    Ncls = num_cells_structured((Ni, Nj, Nk))
-    ext = extent_attribute((Ni, Nj, Nk), extent)
+function vtk_grid(
+        dtype::VTKRectilinearGrid, filename::AbstractString,
+        xs::Tuple{Vararg{AbstractVector, 3}};
+        extent = nothing, kwargs...,
+    )
+    Ns = length.(xs)
+    Npts = prod(Ns)
+    Ncls = num_cells_structured(Ns)
+    ext = extent_attribute(Ns, extent)
 
     xvtk = XMLDocument()
     vtk = DatasetFile(dtype, xvtk, filename, Npts, Ncls; kwargs...)
@@ -24,6 +26,7 @@ function vtk_grid(dtype::VTKRectilinearGrid, filename::AbstractString,
     xPoints = new_child(xPiece, "Coordinates")
 
     # DataArray node
+    x, y, z = xs
     data_to_xml(vtk, xPoints, x, "x")
     data_to_xml(vtk, xPoints, y, "y")
     data_to_xml(vtk, xPoints, z, "z")
@@ -48,11 +51,11 @@ VTK file 'abc.vtr' (RectilinearGrid file, open)
 ```
 
 """
-vtk_grid(filename::AbstractString, x::AbstractVector{T},
-         y::AbstractVector{T}, z::AbstractVector{T}; kwargs...) where T =
-    vtk_grid(VTKRectilinearGrid(), filename, x, y, z; kwargs...)
-
-# 2D variant
-vtk_grid(filename::AbstractString, x::AbstractVector{T},
-         y::AbstractVector{T}; kwargs...) where T =
-    vtk_grid(VTKRectilinearGrid(), filename, x, y, Zeros{T}(1); kwargs...)
+function vtk_grid(
+        filename::AbstractString,
+        x::AbstractVector{T}, y::AbstractVector{T},
+        z::AbstractVector{T} = Zeros{T}(1);
+        kwargs...,
+    ) where {T}
+    vtk_grid(VTKRectilinearGrid(), filename, (x, y, z); kwargs...)
+end

From 61483b70de0035e7f2668654d9edb4be0f25669a Mon Sep 17 00:00:00 2001
From: Juan Ignacio Polanco <jipolanc@gmail.com>
Date: Tue, 18 Jan 2022 11:29:39 +0100
Subject: [PATCH 04/18] Support whole_extent keyword argument

---
 src/gridtypes/structured/imagedata.jl   | 11 ++++++----
 src/gridtypes/structured/rectilinear.jl | 11 ++++++----
 src/gridtypes/structured/structured.jl  | 27 ++++++++++++++++---------
 3 files changed, 31 insertions(+), 18 deletions(-)

diff --git a/src/gridtypes/structured/imagedata.jl b/src/gridtypes/structured/imagedata.jl
index a322cbaa..2efce518 100644
--- a/src/gridtypes/structured/imagedata.jl
+++ b/src/gridtypes/structured/imagedata.jl
@@ -2,11 +2,10 @@ function vtk_grid(
         dtype::VTKImageData, filename::AbstractString, Ns::Dims{N};
         origin::NTuple{N} = ntuple(d -> 0.0, Val(N)),
         spacing::NTuple{N} = ntuple(d -> 1.0, Val(N)),
-        extent = nothing, kwargs...,
+        extent = nothing, whole_extent = extent, kwargs...,
     ) where {N}
     Npts = prod(Ns)
     Ncls = num_cells_structured(Ns)
-    ext = extent_attribute(Ns, extent)
 
     xvtk = XMLDocument()
     vtk = DatasetFile(dtype, xvtk, filename, Npts, Ncls; kwargs...)
@@ -16,7 +15,9 @@ function vtk_grid(
 
     # ImageData node
     xGrid = new_child(xroot, vtk.grid_type)
-    set_attribute(xGrid, "WholeExtent", ext)
+    let ext = extent_attribute(Ns, whole_extent)
+        set_attribute(xGrid, "WholeExtent", ext)
+    end
 
     origin_str = _tuple_to_str3(origin, zero(eltype(origin)))
     spacing_str = _tuple_to_str3(spacing, one(eltype(origin)))
@@ -26,7 +27,9 @@ function vtk_grid(
 
     # Piece node
     xPiece = new_child(xGrid, "Piece")
-    set_attribute(xPiece, "Extent", ext)
+    let ext = extent_attribute(Ns, extent)
+        set_attribute(xPiece, "Extent", ext)
+    end
 
     vtk
 end
diff --git a/src/gridtypes/structured/rectilinear.jl b/src/gridtypes/structured/rectilinear.jl
index 05b66640..792dc07f 100644
--- a/src/gridtypes/structured/rectilinear.jl
+++ b/src/gridtypes/structured/rectilinear.jl
@@ -1,12 +1,11 @@
 function vtk_grid(
         dtype::VTKRectilinearGrid, filename::AbstractString,
         xs::Tuple{Vararg{AbstractVector, 3}};
-        extent = nothing, kwargs...,
+        extent = nothing, whole_extent = extent, kwargs...,
     )
     Ns = length.(xs)
     Npts = prod(Ns)
     Ncls = num_cells_structured(Ns)
-    ext = extent_attribute(Ns, extent)
 
     xvtk = XMLDocument()
     vtk = DatasetFile(dtype, xvtk, filename, Npts, Ncls; kwargs...)
@@ -16,11 +15,15 @@ function vtk_grid(
 
     # RectilinearGrid node
     xGrid = new_child(xroot, vtk.grid_type)
-    set_attribute(xGrid, "WholeExtent", ext)
+    let ext = extent_attribute(Ns, whole_extent)
+        set_attribute(xGrid, "WholeExtent", ext)
+    end
 
     # Piece node
     xPiece = new_child(xGrid, "Piece")
-    set_attribute(xPiece, "Extent", ext)
+    let ext = extent_attribute(Ns, extent)
+        set_attribute(xPiece, "Extent", ext)
+    end
 
     # Coordinates node
     xPoints = new_child(xPiece, "Coordinates")
diff --git a/src/gridtypes/structured/structured.jl b/src/gridtypes/structured/structured.jl
index 1f69bef6..f75a3289 100644
--- a/src/gridtypes/structured/structured.jl
+++ b/src/gridtypes/structured/structured.jl
@@ -10,18 +10,21 @@ structured_dims(xyz::Array4) = ntuple(d -> size(xyz, d + 1), 3)
 structured_dims(xyz::Array3Tuple3) = size(first(xyz))
 structured_dims(xyz::Array3ofVec3) = size(xyz)
 
-function vtk_grid(dtype::VTKStructuredGrid, filename::AbstractString,
-                  xyz::StructuredCoords; extent=nothing, kwargs...)
-    Ni, Nj, Nk = structured_dims(xyz)
-    Npts = Ni * Nj * Nk
+function vtk_grid(
+        dtype::VTKStructuredGrid, filename::AbstractString,
+        xyz::StructuredCoords;
+        extent = nothing, whole_extent = extent, kwargs...,
+    )
+    Ns = structured_dims(xyz)
+    Npts = prod(Ns)
     Ncomp = num_components(xyz, Npts)
-    Ncls = num_cells_structured((Ni, Nj, Nk))
-    ext = extent_attribute((Ni, Nj, Nk), extent)
+    Ncls = num_cells_structured(Ns)
 
     if Ncomp != 3  # three components (x, y, z)
+        # TODO the error message is only accurate when xyz isa Array4
         msg = "coordinate array `xyz` has incorrect dimensions.\n" *
-              "Expected dimensions: $((3, Ni, Nj, Nk)).\n" *
-              "Actual dimensions: $((Ncomp, Ni, Nj, Nk))"
+              "Expected dimensions: $((3, Ns...)).\n" *
+              "Actual dimensions: $((Ncomp, Ns...))"
         throw(DimensionMismatch(msg))
     end
 
@@ -33,11 +36,15 @@ function vtk_grid(dtype::VTKStructuredGrid, filename::AbstractString,
 
     # StructuredGrid node
     xGrid = new_child(xroot, vtk.grid_type)
-    set_attribute(xGrid, "WholeExtent", ext)
+    let ext = extent_attribute(Ns, whole_extent)
+        set_attribute(xGrid, "WholeExtent", ext)
+    end
 
     # Piece node
     xPiece = new_child(xGrid, "Piece")
-    set_attribute(xPiece, "Extent", ext)
+    let ext = extent_attribute(Ns, extent)
+        set_attribute(xPiece, "Extent", ext)
+    end
 
     # Points node
     xPoints = new_child(xPiece, "Points")

From 04af520341c68c80ffd0f7722f90ae89e95c83cf Mon Sep 17 00:00:00 2001
From: Juan Ignacio Polanco <jipolanc@gmail.com>
Date: Tue, 18 Jan 2022 11:35:37 +0100
Subject: [PATCH 05/18] Fix doctests

---
 src/gridtypes/structured/imagedata.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/gridtypes/structured/imagedata.jl b/src/gridtypes/structured/imagedata.jl
index 2efce518..fe726454 100644
--- a/src/gridtypes/structured/imagedata.jl
+++ b/src/gridtypes/structured/imagedata.jl
@@ -78,9 +78,9 @@ VTK file 'def.vti' (ImageData file, open)
 
 """
 function vtk_grid(filename::AbstractString, xyz::Vararg{AbstractRange}; kwargs...)
-    Nxyz = length.(xyz)
-    origin = first.(xyz)
-    spacing = step.(xyz)
+    Nxyz = promote(length.(xyz)...)
+    origin = promote(first.(xyz)...)
+    spacing = promote(step.(xyz)...)
     vtk_grid(VTKImageData(), filename, Nxyz;
              origin=origin, spacing=spacing, kwargs...)
 end

From 1462d99f3b1518f616deedd23c51889785087f71 Mon Sep 17 00:00:00 2001
From: Juan Ignacio Polanco <jipolanc@gmail.com>
Date: Tue, 18 Jan 2022 14:21:23 +0100
Subject: [PATCH 06/18] Disable extent checks for WholeExtent

---
 src/gridtypes/common.jl                 | 38 +++++++++++++++----------
 src/gridtypes/structured/imagedata.jl   |  4 +--
 src/gridtypes/structured/rectilinear.jl |  4 +--
 src/gridtypes/structured/structured.jl  |  4 +--
 test/extent.jl                          |  1 +
 5 files changed, 30 insertions(+), 21 deletions(-)

diff --git a/src/gridtypes/common.jl b/src/gridtypes/common.jl
index 49c1adab..ee4e7a36 100644
--- a/src/gridtypes/common.jl
+++ b/src/gridtypes/common.jl
@@ -17,37 +17,45 @@ end
 # Specifies the "extent" of a structured grid, e.g. (1:4, 0:3, 1:42)
 const Extent{N} = Tuple{Vararg{AbstractUnitRange{<:Integer}, N}} where {N}
 
-# Returns the "extent" attribute required for structured (including rectilinear)
-# grids.
-function extent_attribute(Ns::Dims{3}, extent::Extent{3})
-    all(pair -> pair[1] == length(pair[2]), zip(Ns, extent)) ||
+function maybe_check_extent(Ns::Dims{N}, ext::Extent{N}; check = true) where {N}
+    check || return true
+    all(pair -> pair[1] == length(pair[2]), zip(Ns, ext)) ||
         throw(DimensionMismatch("extent is not consistent with dataset dimensions."))
-    iter = Iterators.map(e -> string(first(e), " ", last(e)), extent)
-    join(iter, " ")
 end
 
-function extent_attribute(Ns::Dims{N}, extent::Extent{N}) where {N}
+function to_extent(Ns::Dims{3}, extent::Extent{3}; kws...)
+    maybe_check_extent(Ns, extent; kws...)
+    extent
+end
+
+function to_extent(Ns::Dims{N}, extent::Extent{N}; kws...) where {N}
     if N > 3
         throw(ArgumentError("dimensionalities N > 3 not supported (got N = $N)"))
     elseif N < 3
+        maybe_check_extent(Ns, extent; kws...)
         M = 3 - N
-        extent_attribute(
-            (Ns..., ntuple(d -> 1, Val(M))...),
-            (extent..., ntuple(d -> 0:0, Val(M))...),
-        )
+        (extent..., ntuple(d -> 0:0, Val(M))...)
     end
 end
 
-function extent_attribute(Ns::Dims, ::Nothing = nothing)
+function to_extent(Ns::Dims, ::Nothing = nothing; kws...)
     ext = map(N -> 0:(N - 1), Ns)
-    extent_attribute(Ns, ext)
+    to_extent(Ns, ext; kws...)
 end
 
 # This is left for compatibility...
-function extent_attribute(Ns::Dims{N}, ext::Array{T}) where {N, T <: Integer}
+function to_extent(Ns::Dims{N}, ext::Array{T}; kws...) where {N, T <: Integer}
     length(ext) == 2N || throw(ArgumentError("extent must have length $(2N)."))
     rs = ntuple(i -> (ext[2i - 1]):(ext[2i]), Val(N))
-    extent_attribute(Ns, rs)
+    to_extent(Ns, rs; kws...)
+end
+
+# Returns the "extent" attribute required for structured (including rectilinear)
+# grids.
+function extent_attribute(Ns::Dims, extent = nothing; kws...)
+    ext = to_extent(Ns, extent; kws...)
+    iter = Iterators.map(e -> string(first(e), " ", last(e)), ext)
+    join(iter, " ")
 end
 
 # Number of cells in structured grids.
diff --git a/src/gridtypes/structured/imagedata.jl b/src/gridtypes/structured/imagedata.jl
index fe726454..a77f9035 100644
--- a/src/gridtypes/structured/imagedata.jl
+++ b/src/gridtypes/structured/imagedata.jl
@@ -15,7 +15,7 @@ function vtk_grid(
 
     # ImageData node
     xGrid = new_child(xroot, vtk.grid_type)
-    let ext = extent_attribute(Ns, whole_extent)
+    let ext = extent_attribute(Ns, whole_extent; check = false)
         set_attribute(xGrid, "WholeExtent", ext)
     end
 
@@ -27,7 +27,7 @@ function vtk_grid(
 
     # Piece node
     xPiece = new_child(xGrid, "Piece")
-    let ext = extent_attribute(Ns, extent)
+    let ext = extent_attribute(Ns, extent; check = true)
         set_attribute(xPiece, "Extent", ext)
     end
 
diff --git a/src/gridtypes/structured/rectilinear.jl b/src/gridtypes/structured/rectilinear.jl
index 792dc07f..3e6ff068 100644
--- a/src/gridtypes/structured/rectilinear.jl
+++ b/src/gridtypes/structured/rectilinear.jl
@@ -15,13 +15,13 @@ function vtk_grid(
 
     # RectilinearGrid node
     xGrid = new_child(xroot, vtk.grid_type)
-    let ext = extent_attribute(Ns, whole_extent)
+    let ext = extent_attribute(Ns, whole_extent; check = false)
         set_attribute(xGrid, "WholeExtent", ext)
     end
 
     # Piece node
     xPiece = new_child(xGrid, "Piece")
-    let ext = extent_attribute(Ns, extent)
+    let ext = extent_attribute(Ns, extent; check = true)
         set_attribute(xPiece, "Extent", ext)
     end
 
diff --git a/src/gridtypes/structured/structured.jl b/src/gridtypes/structured/structured.jl
index f75a3289..dc35aeba 100644
--- a/src/gridtypes/structured/structured.jl
+++ b/src/gridtypes/structured/structured.jl
@@ -36,13 +36,13 @@ function vtk_grid(
 
     # StructuredGrid node
     xGrid = new_child(xroot, vtk.grid_type)
-    let ext = extent_attribute(Ns, whole_extent)
+    let ext = extent_attribute(Ns, whole_extent; check = false)
         set_attribute(xGrid, "WholeExtent", ext)
     end
 
     # Piece node
     xPiece = new_child(xGrid, "Piece")
-    let ext = extent_attribute(Ns, extent)
+    let ext = extent_attribute(Ns, extent; check = true)
         set_attribute(xPiece, "Extent", ext)
     end
 
diff --git a/test/extent.jl b/test/extent.jl
index 033c3d22..3ea891c6 100644
--- a/test/extent.jl
+++ b/test/extent.jl
@@ -6,6 +6,7 @@ using WriteVTK
         Ns = (6, 8)
         ext = (1:6, 3:10)
         @test_throws DimensionMismatch WriteVTK.extent_attribute(Ns .- 1, ext)
+        @test WriteVTK.extent_attribute(Ns .- 1, ext; check = false) == "1 6 3 10 0 0"  # disabling checks
         @test WriteVTK.extent_attribute(Ns, ext) == "1 6 3 10 0 0"
         @test WriteVTK.extent_attribute(Ns) == "0 5 0 7 0 0"
     end

From 771d7088e43b9cbb2379805c317f82aec04517a6 Mon Sep 17 00:00:00 2001
From: Juan Ignacio Polanco <jipolanc@gmail.com>
Date: Thu, 20 Jan 2022 09:28:11 +0100
Subject: [PATCH 07/18] Simplify VTK extents

- remove checks and simplify signatures
- they must now be passed as one-based ranges, and converted to
  zero-based for writing to VTK files
---
 src/gridtypes/common.jl                 | 39 +++++++++----------------
 src/gridtypes/structured/imagedata.jl   |  6 ++--
 src/gridtypes/structured/rectilinear.jl |  8 ++---
 src/gridtypes/structured/structured.jl  |  6 ++--
 test/extent.jl                          | 14 ++++-----
 5 files changed, 29 insertions(+), 44 deletions(-)

diff --git a/src/gridtypes/common.jl b/src/gridtypes/common.jl
index ee4e7a36..f29af329 100644
--- a/src/gridtypes/common.jl
+++ b/src/gridtypes/common.jl
@@ -14,46 +14,35 @@ function vtk_xml_write_header(vtk::DatasetFile)
     xroot::XMLElement
 end
 
-# Specifies the "extent" of a structured grid, e.g. (1:4, 0:3, 1:42)
+# Specifies the "extent" of a structured grid, e.g. (1:4, 1:3, 1:42)
 const Extent{N} = Tuple{Vararg{AbstractUnitRange{<:Integer}, N}} where {N}
 
-function maybe_check_extent(Ns::Dims{N}, ext::Extent{N}; check = true) where {N}
-    check || return true
-    all(pair -> pair[1] == length(pair[2]), zip(Ns, ext)) ||
-        throw(DimensionMismatch("extent is not consistent with dataset dimensions."))
-end
+# Switch to zero-based extent used by VTK.
+to_vtk_extent(r::AbstractUnitRange) = oftype(r, r .- 1)
 
-function to_extent(Ns::Dims{3}, extent::Extent{3}; kws...)
-    maybe_check_extent(Ns, extent; kws...)
-    extent
+function to_vtk_extent(extent::Union{Tuple, AbstractArray})
+    ext = to_extent3(extent)
+    map(to_vtk_extent, ext)
 end
 
-function to_extent(Ns::Dims{N}, extent::Extent{N}; kws...) where {N}
+# Convert different extent specifications to Extent{3}.
+to_extent3(extent::Extent{3}) = extent
+
+function to_extent3(extent::Extent{N}) where {N}
     if N > 3
         throw(ArgumentError("dimensionalities N > 3 not supported (got N = $N)"))
     elseif N < 3
-        maybe_check_extent(Ns, extent; kws...)
         M = 3 - N
-        (extent..., ntuple(d -> 0:0, Val(M))...)
+        (extent..., ntuple(d -> 1:1, Val(M))...) :: Extent{3}
     end
 end
 
-function to_extent(Ns::Dims, ::Nothing = nothing; kws...)
-    ext = map(N -> 0:(N - 1), Ns)
-    to_extent(Ns, ext; kws...)
-end
-
-# This is left for compatibility...
-function to_extent(Ns::Dims{N}, ext::Array{T}; kws...) where {N, T <: Integer}
-    length(ext) == 2N || throw(ArgumentError("extent must have length $(2N)."))
-    rs = ntuple(i -> (ext[2i - 1]):(ext[2i]), Val(N))
-    to_extent(Ns, rs; kws...)
-end
+to_extent3(Ns::Dims) = to_extent3(map(N -> 1:N, Ns))
 
 # Returns the "extent" attribute required for structured (including rectilinear)
 # grids.
-function extent_attribute(Ns::Dims, extent = nothing; kws...)
-    ext = to_extent(Ns, extent; kws...)
+function extent_attribute(extent_in)
+    ext = to_vtk_extent(extent_in)
     iter = Iterators.map(e -> string(first(e), " ", last(e)), ext)
     join(iter, " ")
 end
diff --git a/src/gridtypes/structured/imagedata.jl b/src/gridtypes/structured/imagedata.jl
index a77f9035..3fdde994 100644
--- a/src/gridtypes/structured/imagedata.jl
+++ b/src/gridtypes/structured/imagedata.jl
@@ -2,7 +2,7 @@ function vtk_grid(
         dtype::VTKImageData, filename::AbstractString, Ns::Dims{N};
         origin::NTuple{N} = ntuple(d -> 0.0, Val(N)),
         spacing::NTuple{N} = ntuple(d -> 1.0, Val(N)),
-        extent = nothing, whole_extent = extent, kwargs...,
+        extent = Ns, whole_extent = extent, kwargs...,
     ) where {N}
     Npts = prod(Ns)
     Ncls = num_cells_structured(Ns)
@@ -15,7 +15,7 @@ function vtk_grid(
 
     # ImageData node
     xGrid = new_child(xroot, vtk.grid_type)
-    let ext = extent_attribute(Ns, whole_extent; check = false)
+    let ext = extent_attribute(whole_extent)
         set_attribute(xGrid, "WholeExtent", ext)
     end
 
@@ -27,7 +27,7 @@ function vtk_grid(
 
     # Piece node
     xPiece = new_child(xGrid, "Piece")
-    let ext = extent_attribute(Ns, extent; check = true)
+    let ext = extent_attribute(extent)
         set_attribute(xPiece, "Extent", ext)
     end
 
diff --git a/src/gridtypes/structured/rectilinear.jl b/src/gridtypes/structured/rectilinear.jl
index 3e6ff068..b837a619 100644
--- a/src/gridtypes/structured/rectilinear.jl
+++ b/src/gridtypes/structured/rectilinear.jl
@@ -1,9 +1,9 @@
 function vtk_grid(
         dtype::VTKRectilinearGrid, filename::AbstractString,
         xs::Tuple{Vararg{AbstractVector, 3}};
-        extent = nothing, whole_extent = extent, kwargs...,
+        extent = map(length, xs), whole_extent = extent, kwargs...,
     )
-    Ns = length.(xs)
+    Ns = map(length, xs)
     Npts = prod(Ns)
     Ncls = num_cells_structured(Ns)
 
@@ -15,13 +15,13 @@ function vtk_grid(
 
     # RectilinearGrid node
     xGrid = new_child(xroot, vtk.grid_type)
-    let ext = extent_attribute(Ns, whole_extent; check = false)
+    let ext = extent_attribute(whole_extent)
         set_attribute(xGrid, "WholeExtent", ext)
     end
 
     # Piece node
     xPiece = new_child(xGrid, "Piece")
-    let ext = extent_attribute(Ns, extent; check = true)
+    let ext = extent_attribute(extent)
         set_attribute(xPiece, "Extent", ext)
     end
 
diff --git a/src/gridtypes/structured/structured.jl b/src/gridtypes/structured/structured.jl
index dc35aeba..49fcbaec 100644
--- a/src/gridtypes/structured/structured.jl
+++ b/src/gridtypes/structured/structured.jl
@@ -13,7 +13,7 @@ structured_dims(xyz::Array3ofVec3) = size(xyz)
 function vtk_grid(
         dtype::VTKStructuredGrid, filename::AbstractString,
         xyz::StructuredCoords;
-        extent = nothing, whole_extent = extent, kwargs...,
+        extent = structured_dims(xyz), whole_extent = extent, kwargs...,
     )
     Ns = structured_dims(xyz)
     Npts = prod(Ns)
@@ -36,13 +36,13 @@ function vtk_grid(
 
     # StructuredGrid node
     xGrid = new_child(xroot, vtk.grid_type)
-    let ext = extent_attribute(Ns, whole_extent; check = false)
+    let ext = extent_attribute(whole_extent)
         set_attribute(xGrid, "WholeExtent", ext)
     end
 
     # Piece node
     xPiece = new_child(xGrid, "Piece")
-    let ext = extent_attribute(Ns, extent; check = true)
+    let ext = extent_attribute(extent)
         set_attribute(xPiece, "Extent", ext)
     end
 
diff --git a/test/extent.jl b/test/extent.jl
index 3ea891c6..50adfaf1 100644
--- a/test/extent.jl
+++ b/test/extent.jl
@@ -5,24 +5,20 @@ using WriteVTK
     @testset "N = 2" begin
         Ns = (6, 8)
         ext = (1:6, 3:10)
-        @test_throws DimensionMismatch WriteVTK.extent_attribute(Ns .- 1, ext)
-        @test WriteVTK.extent_attribute(Ns .- 1, ext; check = false) == "1 6 3 10 0 0"  # disabling checks
-        @test WriteVTK.extent_attribute(Ns, ext) == "1 6 3 10 0 0"
         @test WriteVTK.extent_attribute(Ns) == "0 5 0 7 0 0"
+        @test WriteVTK.extent_attribute(ext) == "0 5 2 9 0 0"
     end
 
     @testset "N = 3" begin
         Ns = (6, 8, 42)
-        ext = (0:5, 0:7, 0:2)
-        @test_throws DimensionMismatch WriteVTK.extent_attribute(Ns, ext)
-        ext = (0:5, 0:7, 1:42)
-        @test WriteVTK.extent_attribute(Ns, ext) == "0 5 0 7 1 42"
+        ext = (1:6, 1:8, 1:3)
+        ext = (1:6, 1:8, 2:43)
+        @test WriteVTK.extent_attribute(ext) == "0 5 0 7 1 42"
         @test WriteVTK.extent_attribute(Ns) == "0 5 0 7 0 41"
     end
 
     @testset "N = 4" begin
         Ns = (3, 4, 5, 6)
-        ext = Base.OneTo.(Ns)
-        @test_throws ArgumentError WriteVTK.extent_attribute(Ns, ext)
+        @test_throws ArgumentError WriteVTK.extent_attribute(Ns)
     end
 end

From 6c3262be25686b3e901c7c96b29077ed3c114c9f Mon Sep 17 00:00:00 2001
From: Juan Ignacio Polanco <jipolanc@gmail.com>
Date: Thu, 20 Jan 2022 09:31:03 +0100
Subject: [PATCH 08/18] Support extents using Base.OneTo

---
 src/gridtypes/common.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/gridtypes/common.jl b/src/gridtypes/common.jl
index f29af329..2af2119e 100644
--- a/src/gridtypes/common.jl
+++ b/src/gridtypes/common.jl
@@ -18,7 +18,7 @@ end
 const Extent{N} = Tuple{Vararg{AbstractUnitRange{<:Integer}, N}} where {N}
 
 # Switch to zero-based extent used by VTK.
-to_vtk_extent(r::AbstractUnitRange) = oftype(r, r .- 1)
+to_vtk_extent(r::AbstractUnitRange) = r .- 1
 
 function to_vtk_extent(extent::Union{Tuple, AbstractArray})
     ext = to_extent3(extent)
@@ -33,11 +33,11 @@ function to_extent3(extent::Extent{N}) where {N}
         throw(ArgumentError("dimensionalities N > 3 not supported (got N = $N)"))
     elseif N < 3
         M = 3 - N
-        (extent..., ntuple(d -> 1:1, Val(M))...) :: Extent{3}
+        (extent..., ntuple(d -> Base.OneTo(1), Val(M))...) :: Extent{3}
     end
 end
 
-to_extent3(Ns::Dims) = to_extent3(map(N -> 1:N, Ns))
+to_extent3(Ns::Dims) = to_extent3(map(Base.OneTo, Ns))
 
 # Returns the "extent" attribute required for structured (including rectilinear)
 # grids.

From 84ffd536280255abe52e73593c0b073e007c347e Mon Sep 17 00:00:00 2001
From: Juan Ignacio Polanco <jipolanc@gmail.com>
Date: Thu, 20 Jan 2022 09:31:55 +0100
Subject: [PATCH 09/18] Image data: shift origin according to extent

---
 src/gridtypes/structured/imagedata.jl | 11 ++++++++++-
 1 file changed, 10 insertions(+), 1 deletion(-)

diff --git a/src/gridtypes/structured/imagedata.jl b/src/gridtypes/structured/imagedata.jl
index 3fdde994..d0322226 100644
--- a/src/gridtypes/structured/imagedata.jl
+++ b/src/gridtypes/structured/imagedata.jl
@@ -79,8 +79,17 @@ VTK file 'def.vti' (ImageData file, open)
 """
 function vtk_grid(filename::AbstractString, xyz::Vararg{AbstractRange}; kwargs...)
     Nxyz = promote(length.(xyz)...)
-    origin = promote(first.(xyz)...)
     spacing = promote(step.(xyz)...)
+    origin = promote(first.(xyz)...)
+
+    if (extent = get(kwargs, :extent, nothing)) !== nothing
+        # Shift origin accordingly
+        length.(extent) == Nxyz ||
+            throw(DimensionMismatch("`extent` argument doesn't match grid dimensions"))
+        origin_new = @. origin + (1 - first(extent)) * spacing
+        origin = oftype(origin, origin_new)
+    end
+
     vtk_grid(VTKImageData(), filename, Nxyz;
              origin=origin, spacing=spacing, kwargs...)
 end

From 18f06786dab091065bd68a9abe6d395fc8f09693 Mon Sep 17 00:00:00 2001
From: Juan Ignacio Polanco <jipolanc@gmail.com>
Date: Thu, 20 Jan 2022 09:32:26 +0100
Subject: [PATCH 10/18] Update tests for extent changes

---
 test/imagedata.jl     | 2 +-
 test/pvdCollection.jl | 2 +-
 test/rectilinear.jl   | 2 +-
 test/structured.jl    | 5 +----
 4 files changed, 4 insertions(+), 7 deletions(-)

diff --git a/test/imagedata.jl b/test/imagedata.jl
index f8120973..54a5db82 100644
--- a/test/imagedata.jl
+++ b/test/imagedata.jl
@@ -33,7 +33,7 @@ function main()
     end
 
     # These are all optional:
-    extent = map(r -> r .+ 42, (1:Ni, 1:Nj, 1:Nk))
+    extent = map(r -> r .+ 43, (1:Ni, 1:Nj, 1:Nk))
     origin = (1.2, 4.3, -3.1)
     spacing = (0.1, 0.5, 1.2)
 
diff --git a/test/pvdCollection.jl b/test/pvdCollection.jl
index 43a975ef..aaf5c333 100644
--- a/test/pvdCollection.jl
+++ b/test/pvdCollection.jl
@@ -54,7 +54,7 @@ function main()
     cdata = zeros(FloatType, Ni - 1, Nj - 1, Nk - 1)
 
     # Test extents (this is optional!!)
-    ext = [0, Ni-1, 0, Nj-1, 0, Nk-1] .+ 42
+    ext = map(N -> (1:N) .+ 42, (Ni, Nj, Nk))
 
     # Initialise pvd container file
     @time outfiles = paraview_collection(vtk_filename_noext) do pvd
diff --git a/test/rectilinear.jl b/test/rectilinear.jl
index d21a9821..f31d20fa 100644
--- a/test/rectilinear.jl
+++ b/test/rectilinear.jl
@@ -72,7 +72,7 @@ function main()
         end
 
         # Test extents (this is optional!!)
-        ext = [0, Ni-1, 0, Nj-1, 0, Nk-1] .+ 42
+        ext = map(N -> (1:N) .+ 42, (Ni, Nj, Nk))
 
         # Initialise new vtr file (rectilinear grid).
         @time begin
diff --git a/test/structured.jl b/test/structured.jl
index d9dc3c78..413599f0 100644
--- a/test/structured.jl
+++ b/test/structured.jl
@@ -69,9 +69,6 @@ function generate_structured(grid_format, ::Val{dim}) where {dim}
         end
     end
 
-    # Test extents (this is optional!!)
-    ext = [0, Ni-1, 0, Nj-1, 0, Nk-1]
-
     # We test defining grid points using different formats.
     points = if grid_format === SingleArray()
         xyz
@@ -96,7 +93,7 @@ function generate_structured(grid_format, ::Val{dim}) where {dim}
     @time begin
         # Initialise new vts file (structured grid).
         fname = "$(vtk_filename_noext)_$(dim)D_$(name(grid_format))"
-        vtk = vtk_grid(fname, points; extent=ext)
+        vtk = vtk_grid(fname, points)
 
         # Add data.
         vtk["p_values"] = psub

From 10d0135ecd171e921441a83c4d3e4b2ecec3b14d Mon Sep 17 00:00:00 2001
From: Juan Ignacio Polanco <jipolanc@gmail.com>
Date: Thu, 20 Jan 2022 09:36:59 +0100
Subject: [PATCH 11/18] Extents: remove unused signature

---
 src/gridtypes/common.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/gridtypes/common.jl b/src/gridtypes/common.jl
index 2af2119e..49a388b5 100644
--- a/src/gridtypes/common.jl
+++ b/src/gridtypes/common.jl
@@ -20,7 +20,7 @@ const Extent{N} = Tuple{Vararg{AbstractUnitRange{<:Integer}, N}} where {N}
 # Switch to zero-based extent used by VTK.
 to_vtk_extent(r::AbstractUnitRange) = r .- 1
 
-function to_vtk_extent(extent::Union{Tuple, AbstractArray})
+function to_vtk_extent(extent::Tuple)
     ext = to_extent3(extent)
     map(to_vtk_extent, ext)
 end

From e57731a05cd71357bd1303a862c02fec608c47be Mon Sep 17 00:00:00 2001
From: Corentin Lothode <corentin.lothode@univ-rouen.fr>
Date: Thu, 13 Jan 2022 18:09:53 +0100
Subject: [PATCH 12/18] Adding Coordinates, WholeExtent and Extent that were
 missing for Structured Parallel files

---
 src/gridtypes/pvtk_grid.jl | 47 ++++++++++++++++++++++++++++++--------
 1 file changed, 37 insertions(+), 10 deletions(-)

diff --git a/src/gridtypes/pvtk_grid.jl b/src/gridtypes/pvtk_grid.jl
index 5ebc1bd7..c36d97ea 100644
--- a/src/gridtypes/pvtk_grid.jl
+++ b/src/gridtypes/pvtk_grid.jl
@@ -116,18 +116,45 @@ function _init_pvtk!(pvtk::PVTKFile)
     # Recover point type and number of components
     vtk_root = root(vtk.xdoc)
     vtk_grid = find_element(vtk_root, vtk.grid_type)
+    # adding whole extent if necessary
+    whole_extent = attribute(vtk_grid, "WholeExtent")
+    if ~isnothing(whole_extent)
+        set_attribute(pvtk_grid, "WholeExtent", whole_extent)
+    end
     vtk_piece = find_element(vtk_grid, "Piece")
+    # adding extent if necessary
+    extent = attribute(vtk_piece, "Extent")
+    if ~isnothing(extent)
+        for c in child_elements(pvtk_grid)
+            set_attribute(c, "Extent", extent)
+        end
+    end
+
+    # If serial VTK has points
     vtk_points = find_element(vtk_piece, "Points")
-    vtk_data_array = find_element(vtk_points, "DataArray")
-    point_type = attribute(vtk_data_array, "type")
-    Nc = attribute(vtk_data_array, "NumberOfComponents")
-
-    ## PPoints node
-    pvtk_ppoints = new_child(pvtk_grid, "PPoints")
-    pvtk_pdata_array = new_child(pvtk_ppoints, "PDataArray")
-    set_attribute(pvtk_pdata_array, "type", point_type)
-    set_attribute(pvtk_pdata_array, "Name", "Points")
-    set_attribute(pvtk_pdata_array, "NumberOfComponents", Nc)
+    if ~isnothing(vtk_points)
+        vtk_data_array = find_element(vtk_points, "DataArray")
+        point_type = attribute(vtk_data_array, "type")
+        Nc = attribute(vtk_data_array, "NumberOfComponents")
+        ## PPoints node
+        pvtk_ppoints = new_child(pvtk_grid, "PPoints")
+        pvtk_pdata_array = new_child(pvtk_ppoints, "PDataArray")
+        set_attribute(pvtk_pdata_array, "type", point_type)
+        set_attribute(pvtk_pdata_array, "Name", "Points")
+        set_attribute(pvtk_pdata_array, "NumberOfComponents", Nc)
+    end
+
+    # If serial VTK has coordinates
+    vtk_coordinates = find_element(vtk_piece, "Coordinates")
+    if ~isnothing(vtk_coordinates)
+        pvtk_pcoordinates = new_child(pvtk_grid, "PCoordinates")
+        for c in child_elements(vtk_coordinates)
+            pvtk_pdata_array = new_child(pvtk_pcoordinates, "PDataArray")
+            set_attribute(pvtk_pdata_array, "type", attribute(c, "type"))
+            set_attribute(pvtk_pdata_array, "Name", attribute(c, "Name"))
+            set_attribute(pvtk_pdata_array, "NumberOfComponents", attribute(c, "NumberOfComponents"))
+        end
+    end
 
     pvtk
 end

From 55b5407cda600a0313b105cd77bf964840116ec6 Mon Sep 17 00:00:00 2001
From: Corentin Lothode <corentin.lothode@univ-rouen.fr>
Date: Thu, 13 Jan 2022 22:06:08 +0100
Subject: [PATCH 13/18] PImageData support

---
 src/gridtypes/pvtk_grid.jl | 13 +++++++++++++
 1 file changed, 13 insertions(+)

diff --git a/src/gridtypes/pvtk_grid.jl b/src/gridtypes/pvtk_grid.jl
index c36d97ea..fd8a0b5c 100644
--- a/src/gridtypes/pvtk_grid.jl
+++ b/src/gridtypes/pvtk_grid.jl
@@ -115,12 +115,25 @@ function _init_pvtk!(pvtk::PVTKFile)
 
     # Recover point type and number of components
     vtk_root = root(vtk.xdoc)
+
+    # Getting original grid informations
     vtk_grid = find_element(vtk_root, vtk.grid_type)
     # adding whole extent if necessary
     whole_extent = attribute(vtk_grid, "WholeExtent")
     if ~isnothing(whole_extent)
         set_attribute(pvtk_grid, "WholeExtent", whole_extent)
     end
+    # adding origin and spacing if necessary
+    origin = attribute(vtk_grid, "Origin")
+    if ~isnothing(origin)
+        set_attribute(pvtk_grid, "Origin", origin)
+    end
+    spacing = attribute(vtk_grid, "Spacing")
+    if ~isnothing(origin)
+        set_attribute(pvtk_grid, "Spacing", spacing)
+    end
+
+    # Getting original piece informations
     vtk_piece = find_element(vtk_grid, "Piece")
     # adding extent if necessary
     extent = attribute(vtk_piece, "Extent")

From 52fcd6970bbe0428e67b022bbc6da8bc5e31e8c8 Mon Sep 17 00:00:00 2001
From: Juan Ignacio Polanco <jipolanc@gmail.com>
Date: Fri, 21 Jan 2022 10:50:32 +0100
Subject: [PATCH 14/18] Support pvti (parallel image data) files

---
 src/gridtypes/pvtk_grid.jl | 123 ++++++++++++++++++++++++++++---------
 1 file changed, 95 insertions(+), 28 deletions(-)

diff --git a/src/gridtypes/pvtk_grid.jl b/src/gridtypes/pvtk_grid.jl
index fd8a0b5c..fdaa325e 100644
--- a/src/gridtypes/pvtk_grid.jl
+++ b/src/gridtypes/pvtk_grid.jl
@@ -14,8 +14,33 @@ struct PVTKFile <: VTKFile
     path::String
 end
 
+# Returns true if the arguments do *not* contain any cell vectors.
+_pvtk_is_structured(x::AbstractVector{<:AbstractMeshCell}, args...) = Val(false)
+_pvtk_is_structured(x, args...) = _pvtk_is_structured(args...)
+_pvtk_is_structured() = Val(true)
+
+_pvtk_nparts(structured::Val{false}; nparts::Integer, etc...) = nparts
+_pvtk_nparts(structured::Val{true}; extents::AbstractArray, etc...) = length(extents)
+
+_pvtk_extents(structured::Val{false}; etc...) = nothing
+_pvtk_extents(structured::Val{true}; extents::AbstractArray, etc...) = extents
+
+# Filter keyword arguments to be passed to vtk_grid.
+_remove_parallel_kwargs(; nparts = nothing, extents = nothing, kws...) =
+    NamedTuple(kws)
+
+# Determine whole extent from the local extents associated to each process.
+function compute_whole_extent(extents::AbstractArray{<:Extent})
+    # Compute the minimum and maximum across each dimension
+    begins = mapreduce(ext -> map(first, ext), min, extents)
+    ends = mapreduce(ext -> map(last, ext), max, extents)
+    map((b, e) -> b:e, begins, ends)
+end
+
+compute_whole_extent(::Nothing) = nothing
+
 """
-    pvtk_grid(args...; part, nparts, ismain = (part == 1), ghost_level = 0, kwargs...)
+    pvtk_grid(args...; part, nparts, extents, ismain = (part == 1), ghost_level = 0, kwargs...)
 
 Returns a handler representing a parallel vtk file, which can be
 eventually written to file with `vtk_save`.
@@ -28,29 +53,71 @@ The extra keyword arguments only apply to parallel VTK file formats.
 
 Mandatory ones are:
 
-- `part` current (1-based) part id,
-- `nparts` total number of parts,
+- `part`: current (1-based) part id,
+- `nparts`: total number of parts (only **unstructured** grids),
+- `extents`: array specifying the partitioning of a **structured** grid across
+  different processes (see below for details).
 
 Optional ones are:
 
 - `ismain` true if the current part id `part` is the main (the only one that will write the `.pvtk` file),
 - `ghost_level` ghost level.
+
+## Specifying extents for a structured grid
+
+For structured grids, the partitioning of the dataset across different processes
+must be specified via the `extents` argument.
+This is an array where each element represents the data extent associated to a
+given process.
+
+For instance, for a dataset of global dimensions `15×12×4` distributed across 4
+processes, this array may look like the following:
+
+```julia
+extents = [
+    ( 1:10,  1:5, 1:4),  # process 1
+    (10:15,  1:5, 1:4),  # process 2
+    ( 1:10, 5:12, 1:4),  # process 3
+    (10:15, 5:12, 1:4),  # process 4
+]
+```
+
+Some important notes:
+
+- the `extents` argument must be the same for all processes;
+- the extents **must overlap**, or VTK / ParaView will complain when trying to
+  open the files.
+
 """
 function pvtk_grid(
         filename::AbstractString, args...;
-        part, nparts, ismain = (part == 1), ghost_level = 0, kwargs...,
+        part, ismain = (part == 1), ghost_level = 0, kwargs...,
     )
-    bname = basename(filename)
+    is_structured = _pvtk_is_structured(args...)
+    nparts = _pvtk_nparts(is_structured; kwargs...)
+    extents = _pvtk_extents(is_structured; kwargs...)
+
     mkpath(filename)
+    bname = basename(filename)
     prefix = joinpath(filename, bname)
     fn = _serial_filename(part, nparts, prefix, "")
+
+    vtk = let kws_vtk = _remove_parallel_kwargs(; kwargs...)
+        kws = if extents === nothing
+            kws_vtk
+        else
+            (; kws_vtk..., extent = extents[part])
+        end
+        vtk_grid(fn, args...; kws...)
+    end
+
     pvtkargs = PVTKArgs(part, nparts, ismain, ghost_level)
-    xdoc  = XMLDocument()
-    vtk = vtk_grid(fn, args...; kwargs...)
+    xdoc = XMLDocument()
     _, ext = splitext(vtk.path)
     path = filename * ".p" * ext[2:end]
     pvtk = PVTKFile(pvtkargs, xdoc, vtk, path)
-    _init_pvtk!(pvtk)
+    _init_pvtk!(pvtk, extents)
+
     pvtk
 end
 
@@ -82,7 +149,7 @@ function _serial_filename(part, nparts, prefix, extension)
     prefix * "_$p" * extension
 end
 
-function _init_pvtk!(pvtk::PVTKFile)
+function _init_pvtk!(pvtk::PVTKFile, extents)
     # Recover some data
     vtk = pvtk.vtk
     pvtkargs = pvtk.pvtkargs
@@ -107,45 +174,45 @@ function _init_pvtk!(pvtk::PVTKFile)
     set_attribute(pvtk_grid, "GhostLevel", string(pvtkargs.ghost_level))
 
     # Pieces (i.e. Pointers to serial files)
-    for piece in 1:npieces
+    for piece ∈ 1:npieces
         pvtk_piece = new_child(pvtk_grid, "Piece")
         fn = _serial_filename(piece, npieces, prefix, ext)
         set_attribute(pvtk_piece, "Source", fn)
+
+        # Add local extent if necessary
+        if extents !== nothing
+            set_attribute(pvtk_piece, "Extent", extent_attribute(extents[piece]))
+        end
     end
 
-    # Recover point type and number of components
-    vtk_root = root(vtk.xdoc)
+    # Add whole extent if necessary
+    whole_extent = compute_whole_extent(extents)
+    if whole_extent !== nothing
+        set_attribute(pvtk_grid, "WholeExtent", extent_attribute(whole_extent))
+    end
 
     # Getting original grid informations
+    # Recover point type and number of components
+    vtk_root = root(vtk.xdoc)
     vtk_grid = find_element(vtk_root, vtk.grid_type)
-    # adding whole extent if necessary
-    whole_extent = attribute(vtk_grid, "WholeExtent")
-    if ~isnothing(whole_extent)
-        set_attribute(pvtk_grid, "WholeExtent", whole_extent)
-    end
+
     # adding origin and spacing if necessary
     origin = attribute(vtk_grid, "Origin")
-    if ~isnothing(origin)
+    if origin !== nothing
         set_attribute(pvtk_grid, "Origin", origin)
     end
+
     spacing = attribute(vtk_grid, "Spacing")
-    if ~isnothing(origin)
+    if spacing !== nothing
         set_attribute(pvtk_grid, "Spacing", spacing)
     end
 
     # Getting original piece informations
     vtk_piece = find_element(vtk_grid, "Piece")
-    # adding extent if necessary
-    extent = attribute(vtk_piece, "Extent")
-    if ~isnothing(extent)
-        for c in child_elements(pvtk_grid)
-            set_attribute(c, "Extent", extent)
-        end
-    end
 
     # If serial VTK has points
     vtk_points = find_element(vtk_piece, "Points")
-    if ~isnothing(vtk_points)
+    if vtk_points !== nothing
         vtk_data_array = find_element(vtk_points, "DataArray")
         point_type = attribute(vtk_data_array, "type")
         Nc = attribute(vtk_data_array, "NumberOfComponents")
@@ -159,7 +226,7 @@ function _init_pvtk!(pvtk::PVTKFile)
 
     # If serial VTK has coordinates
     vtk_coordinates = find_element(vtk_piece, "Coordinates")
-    if ~isnothing(vtk_coordinates)
+    if vtk_coordinates !== nothing
         pvtk_pcoordinates = new_child(pvtk_grid, "PCoordinates")
         for c in child_elements(vtk_coordinates)
             pvtk_pdata_array = new_child(pvtk_pcoordinates, "PDataArray")

From 70557390f4fa170f352d882760fb8a0834335e29 Mon Sep 17 00:00:00 2001
From: Juan Ignacio Polanco <jipolanc@gmail.com>
Date: Fri, 21 Jan 2022 10:57:32 +0100
Subject: [PATCH 15/18] Add test for pvti files

---
 test/checksums.sha1 |  7 +++++++
 test/pvtk_grid.jl   | 42 +++++++++++++++++++++++++++++++++++++++++-
 test/runtests.jl    |  3 +--
 3 files changed, 49 insertions(+), 3 deletions(-)

diff --git a/test/checksums.sha1 b/test/checksums.sha1
index 1daeec1c..2652f83e 100644
--- a/test/checksums.sha1
+++ b/test/checksums.sha1
@@ -52,3 +52,10 @@ e0180e5cf18bb7c6f00df0bb9ded769f90b2574e  arraydata_3D.vti
 9e21f61ebf8641fdcc325836d4bd707e52932441  arrays.vti
 2626bef24152fcfe28c810ca68adf473c8fad90c  simulation.pvtu
 9a0c9ff25c52c807adac85d7fdd1001d18422952  simulation/simulation_1.vtu
+e23c8f0db331e1db635aa5d98b9f1375132e7ffd  pimage.pvti
+e581d4824d67655584f573743f42788fb26e162c  pimage/pimage_1.vti
+150ce17cb8ddbb974344954a3193a3373ee1c3fc  pimage/pimage_2.vti
+e389b7948af0bd9cd9d178dcc6d0fb86493738f3  pimage/pimage_3.vti
+669fdb06ab07c289cb5c01dd511da380fb06bdfd  pimage/pimage_4.vti
+1ea223999c053e38ab84726c0a4bfdd0468435ef  pimage/pimage_5.vti
+21dd9d15e51720f9d5822bf663c718974094d779  pimage/pimage_6.vti
diff --git a/test/pvtk_grid.jl b/test/pvtk_grid.jl
index c56a1524..248af02a 100644
--- a/test/pvtk_grid.jl
+++ b/test/pvtk_grid.jl
@@ -1,7 +1,8 @@
+using Base.Threads: @spawn
 using WriteVTK
 using Test
 
-function main()
+function pvtk_unstructured()
     # Suppose that the mesh is made of 5 points:
     cells = [
         MeshCell(VTKCellTypes.VTK_TRIANGLE, [1, 4, 2]),
@@ -24,4 +25,43 @@ function main()
     outfiles
 end
 
+function make_structured_partition()
+    Ns = (10, 15, 3)
+    xp = (1:5, 5:10)  # NOTE: overlapping is required by VTK
+    yp = (1:3, 3:12, 12:15)
+    zp = (1:3, )
+    extents = collect(Iterators.product(xp, yp, zp))
+    Ns, extents
+end
+
+function pvtk_imagedata()
+    Ns, extents = make_structured_partition()
+    xs_whole = map(N -> range(-1; step = 0.2, length = N), Ns)  # full grid
+    nparts = length(extents)  # number of "processes"
+    filenames = Vector{Vector{String}}(undef, nparts)
+    @sync for (n, extent) ∈ enumerate(extents)
+        @spawn begin
+            xs = getindex.(xs_whole, extent)  # local grid
+            point_data = map(x -> +(x...), Iterators.product(xs...))
+            processid = fill(n, length.(xs) .- 1)  # cell data
+            filenames[n] = pvtk_grid(
+                "pimage", xs...;
+                part = n, extents = extents,
+                append = false, compress = false,
+            ) do vtk
+                vtk["point_data"] = point_data
+                vtk["process_id"] = processid
+            end
+        end
+    end
+    collect(Iterators.flatten(filenames))
+end
+
+function main()
+    vcat(
+        pvtk_unstructured(),
+        pvtk_imagedata(),
+    )
+end
+
 main()
diff --git a/test/runtests.jl b/test/runtests.jl
index ed0c6fe9..916d0da1 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,5 +1,3 @@
-#!/usr/bin/env julia
-
 using WriteVTK
 using Test: @test
 using SHA: sha1
@@ -45,6 +43,7 @@ for test in tests
             write(csio, sha_str)
         else
             # Returns `nothing` if string is not found.
+            @info "Verifying $file"
             @test findfirst(sha_str, checksum_list) !== nothing
         end
     end

From 3f596d9549327bdb81e83a9b146db126c58ee8cb Mon Sep 17 00:00:00 2001
From: Juan Ignacio Polanco <jipolanc@gmail.com>
Date: Fri, 21 Jan 2022 11:57:47 +0100
Subject: [PATCH 16/18] Add tests for .pvtr and .pvts

---
 test/checksums.sha1 | 14 +++++++++++
 test/pvtk_grid.jl   | 57 +++++++++++++++++++++++++++++++++++++++++++--
 2 files changed, 69 insertions(+), 2 deletions(-)

diff --git a/test/checksums.sha1 b/test/checksums.sha1
index 2652f83e..260c6b4a 100644
--- a/test/checksums.sha1
+++ b/test/checksums.sha1
@@ -59,3 +59,17 @@ e389b7948af0bd9cd9d178dcc6d0fb86493738f3  pimage/pimage_3.vti
 669fdb06ab07c289cb5c01dd511da380fb06bdfd  pimage/pimage_4.vti
 1ea223999c053e38ab84726c0a4bfdd0468435ef  pimage/pimage_5.vti
 21dd9d15e51720f9d5822bf663c718974094d779  pimage/pimage_6.vti
+d35ceff83c6c43b43eefe07007cbba8752cead15  prectilinear.pvtr
+b7edbc7bb8907b472551d456f60cbc7f0c915b53  prectilinear/prectilinear_1.vtr
+d26dd7cf0570b1fbad8714e16e784f09bb6affed  prectilinear/prectilinear_2.vtr
+e36a533248051994e7da804bff1b168d2d866bba  prectilinear/prectilinear_3.vtr
+4b73039f56d1a92fe958a98c86dde64b3fb86b63  prectilinear/prectilinear_4.vtr
+87a07f8911949b002501da5bdd5daf497a65d09a  prectilinear/prectilinear_5.vtr
+4fd1d4de3867230a0636bef0260417ed17b7c17f  prectilinear/prectilinear_6.vtr
+fc3a26b061585c093520bbb720714d180ca22335  pstructured.pvts
+584520780f94c58a9919cb216d116733f8466e0f  pstructured/pstructured_1.vts
+0ff33138288a96f8f6f4377f2e061748aa702d97  pstructured/pstructured_2.vts
+e365dc765a17a288640e04f726cf8a7421559f39  pstructured/pstructured_3.vts
+0daf5cccfaec20cdcfae94294513ef211954c9e1  pstructured/pstructured_4.vts
+a2b4e0e63dda63eee19466b08d5d4ae8d34dcfd5  pstructured/pstructured_5.vts
+0adbea4b9bebd2b95569398666028c19f3d8cb4a  pstructured/pstructured_6.vts
diff --git a/test/pvtk_grid.jl b/test/pvtk_grid.jl
index 248af02a..a2cdd014 100644
--- a/test/pvtk_grid.jl
+++ b/test/pvtk_grid.jl
@@ -1,4 +1,5 @@
 using Base.Threads: @spawn
+using StaticArrays
 using WriteVTK
 using Test
 
@@ -20,7 +21,6 @@ function pvtk_unstructured()
     end
     @test isfile(vtufile)
     @test vtufile ∈ outfiles
-    println("Saved:  ", join(outfiles, "  "))
     @test WriteVTK._serial_filename(3, 100, "prefix", ".ext") == "prefix_003.ext"
     outfiles
 end
@@ -42,7 +42,7 @@ function pvtk_imagedata()
     @sync for (n, extent) ∈ enumerate(extents)
         @spawn begin
             xs = getindex.(xs_whole, extent)  # local grid
-            point_data = map(x -> +(x...), Iterators.product(xs...))
+            point_data = map(sum, Iterators.product(xs...))
             processid = fill(n, length.(xs) .- 1)  # cell data
             filenames[n] = pvtk_grid(
                 "pimage", xs...;
@@ -57,10 +57,63 @@ function pvtk_imagedata()
     collect(Iterators.flatten(filenames))
 end
 
+function pvtk_rectilinear()
+    Ns, extents = make_structured_partition()
+    nparts = length(extents)  # number of "processes"
+    filenames = Vector{Vector{String}}(undef, nparts)
+    @sync for (n, extent) ∈ enumerate(extents)
+        @spawn begin
+            xs = map(is -> [sqrt(i) for i ∈ is], extent)  # local grid
+            point_data = map(sum, Iterators.product(xs...))
+            processid = fill(n, length.(xs) .- 1)  # cell data
+            filenames[n] = pvtk_grid(
+                "prectilinear", xs...;
+                part = n, extents = extents,
+                append = false, compress = false,
+            ) do vtk
+                vtk["point_data"] = point_data
+                vtk["process_id"] = processid
+            end
+        end
+    end
+    collect(Iterators.flatten(filenames))
+end
+
+function pvtk_structured()
+    Ns, extents = make_structured_partition()
+    nparts = length(extents)  # number of "processes"
+    filenames = Vector{Vector{String}}(undef, nparts)
+    @sync for (n, extent) ∈ enumerate(extents)
+        @spawn begin
+            points = [
+                SVector(
+                    I[1] + sqrt(I[2]),
+                    I[2] + sqrt(I[1]),
+                    I[3] + sqrt(I[1] + I[2]),
+                )
+                for I ∈ CartesianIndices(extent)
+            ]
+            point_data = map(sum, points)
+            processid = fill(n, length.(extent) .- 1)  # cell data
+            filenames[n] = pvtk_grid(
+                "pstructured", points;
+                part = n, extents = extents,
+                append = false, compress = false,
+            ) do vtk
+                vtk["point_data"] = point_data
+                vtk["process_id"] = processid
+            end
+        end
+    end
+    collect(Iterators.flatten(filenames))
+end
+
 function main()
     vcat(
         pvtk_unstructured(),
         pvtk_imagedata(),
+        pvtk_rectilinear(),
+        pvtk_structured(),
     )
 end
 

From 1974f60a1dc016a45141a24f426a23380c71c919 Mon Sep 17 00:00:00 2001
From: Juan Ignacio Polanco <jipolanc@gmail.com>
Date: Fri, 21 Jan 2022 12:35:21 +0100
Subject: [PATCH 17/18] Bump version

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 5a090862..6923c7c8 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "WriteVTK"
 uuid = "64499a7a-5c06-52f2-abe2-ccb03c286192"
 authors = ["Juan Ignacio Polanco <jipolanc@gmail.com>"]
-version = "1.13.0"
+version = "1.14.0"
 
 [deps]
 Base64 = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

From 17fe736e67cd8e8b1a1ae34dbbe6c9d4a573c328 Mon Sep 17 00:00:00 2001
From: Juan Ignacio Polanco <jipolanc@gmail.com>
Date: Fri, 21 Jan 2022 14:48:13 +0100
Subject: [PATCH 18/18] Add docs for parallel structured grids

---
 README.md                     |   2 +-
 docs/src/index.md             |   2 +-
 docs/src/metadata/parallel.md | 164 ++++++++++++++++++++++++++--------
 src/gridtypes/pvtk_grid.jl    |  29 +++---
 4 files changed, 146 insertions(+), 51 deletions(-)

diff --git a/README.md b/README.md
index fe0a74d1..aa186949 100644
--- a/README.md
+++ b/README.md
@@ -73,7 +73,7 @@ Supported dataset formats include:
 Moreover, the following metadata formats are supported:
 - [multiblock files](https://jipolanco.github.io/WriteVTK.jl/dev/metadata/multiblock/) (`.vtm`),
 - [ParaView collections](https://jipolanco.github.io/WriteVTK.jl/dev/metadata/paraview_collections/) (`.pvd`, typically used for time series),
-- [parallel files](https://jipolanco.github.io/WriteVTK.jl/dev/metadata/parallel/) (`.pvt*`, only partial support for now).
+- [parallel files](https://jipolanco.github.io/WriteVTK.jl/dev/metadata/parallel/) (`.pvt*`).
 
 ## Further reading
 
diff --git a/docs/src/index.md b/docs/src/index.md
index 93c29427..2773f331 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -66,7 +66,7 @@ Supported dataset formats include:
 Moreover, the following metadata formats are supported:
 - [multiblock files](@ref Multiblock-files) (`.vtm`),
 - [ParaView collections](@ref ParaView-collections) (`.pvd`, typically used for time series),
-- [parallel files](@ref Parallel-files) (`.pvt*`, only partial support for now).
+- [parallel files](@ref Parallel-files) (`.pvt*`).
 
 ## Authors
 
diff --git a/docs/src/metadata/parallel.md b/docs/src/metadata/parallel.md
index 98f6dd69..413531a1 100644
--- a/docs/src/metadata/parallel.md
+++ b/docs/src/metadata/parallel.md
@@ -3,85 +3,173 @@
 The parallel file formats do not actually store any data in the file.
 Instead, the data is broken into pieces, each of which is stored in a serial file,
 and an extra header file is created containing pointers to the corresponding serial files.
-The header file extension is the serial extension pre-appended with a `p`.
+The header file extension is the serial extension prepended with a `p`.
 For instance, for serial `vtu` files, the corresponding header file extension is `pvtu`.
 
-!!! note "Supported dataset types"
-
-    For now, parallel file formats have only been tested for [unstructured grid](@ref Unstructured-grid) formats (`.vtu`), and are currently unlikely to work for other kinds of datasets.
-    Feel free to open an issue or to submit a PR to add support for other dataset types.
-
 ## Generating a parallel data file
 
 The parallel header file and the corresponding serial files are generated
-using the [`pvtk_grid`](@ref) function.
+by a single call to [`pvtk_grid`](@ref).
 Its signature is
 
 ```julia
 pvtk_grid(
-    args...;
-    part,
-    nparts,
-    ismain = (part == 1),
-    ghost_level = 0,
-    kwargs...,
+    filename, args...;
+    part, nparts, ismain = (part == 1), ghost_level = 0, kwargs...,
 )
 ```
+
 which returns a handler representing a parallel VTK file that can be
 appended with cell and point data and eventually written to disk with
 [`vtk_save`](@ref) as usual.
 In an MPI job, `vtk_save` will cause each rank to write a serial file and just
 a single rank (e.g., rank 0) will write the header file.
 
-Positional and keyword arguments in `args` and `kwargs`
-are passed to `vtk_grid` verbatim in order to generate the serial files
-(with the exception of file names that are augmented with the
-corresponding part id).
+This signature is valid for **unstructured grids**.
+For the case of structured grids, there are small differences detailed [further below](#Parallel-structured-files).
+
+Positional and keyword arguments in `args` and `kwargs` are passed to `vtk_grid`
+verbatim.
+Note that serial filenames are automatically generated from `filename` and from
+the process id `part`.
 
-The extra keyword arguments only apply to parallel VTK file formats.
+The following keyword arguments only apply to parallel VTK file formats.
 
 Mandatory ones are:
 
-- `part` current (1-based) part id (typically MPI rank + 1),
-- `nparts` total number of parts (typically the MPI communicator size).
+- `part`: current (1-based) part id (typically MPI rank + 1),
+- `nparts`: total number of parts (typically the MPI communicator size).
 
 Optional ones are:
 
-- `ismain` true if the current part id `part` is the main (the only one that will write the header file),
-- `ghost_level` ghost level.
+- `ismain`: `true` if the current part id `part` is the main (the only one that will write the header file),
+- `ghost_level`: ghost level.
 
-## Example
+## Parallel structured files
 
-This generates the header file and a single serial file:
+For structured grids, one needs to specify the portion of the grid associated to each process.
+This is done via the `extents` keyword argument, which must be an array containing the data ranges along each dimension associated to each process.
+For example, for a dataset of global dimensions ``15×12×4`` distributed across 4
+processes, this array may look like the following:
 
 ```julia
-cells = [
-    MeshCell(VTKCellTypes.VTK_TRIANGLE, [1, 4, 2]),
-    MeshCell(VTKCellTypes.VTK_QUAD, [2, 4, 3, 5]),
+extents = [
+    ( 1:10,  1:5, 1:4),  # process 1
+    (10:15,  1:5, 1:4),  # process 2
+    ( 1:10, 5:12, 1:4),  # process 3
+    (10:15, 5:12, 1:4),  # process 4
 ]
+```
+
+In practice, in parallel applications, all processes need to have this
+information, and the `extents` argument must be the same for all processes.
+Also note that the length of the `extents` array gives the number of processes,
+and therefore the `nparts` argument is redundant and not needed for structured
+grids.
+
+Finally, note that in the example above the extents for different processes
+**overlap** (for instance, the ranges `1:10` and `10:15` overlap at the index
+`i = 10`).
+This is a requirement of VTK, and without it the full data cannot be visualised
+in ParaView.
+For MPI applications, this typically means that ghost data need to be exchanged
+before writing VTK files.
+
+## Example 1: Unstructured files
+
+This generates two serial files (typically held by two different processes) and
+a header file combining them:
+
+```julia
+all_data = [
+    # Process 1
+    (
+        points = rand(3, 5),  # 5 points on process 1
+        cells = [             # 2 cells  on process 1
+            MeshCell(VTKCellTypes.VTK_TRIANGLE, [1, 4, 2]),
+            MeshCell(VTKCellTypes.VTK_QUAD, [2, 4, 3, 5]),
+        ],
+    ),
+
+    # Process 2
+    (
+        points = rand(3, 4),  # 4 points on process 2
+        cells = [             # 1 cell   on process 2
+            MeshCell(VTKCellTypes.VTK_QUAD, [1, 2, 3, 4]),
+        ]
+    ),
+]
+
+saved_files = Vector{Vector{String}}(undef, 2)  # files saved by each "process"
 
-x = rand(5)
-y = rand(5)
+for part = 1:2
+    data = all_data[part]
+    saved_files[part] = pvtk_grid(
+            "simulation", data.points, data.cells;
+            part = part, nparts = 2,
+        ) do pvtk
+        pvtk["Pressure"] = sum(data.points; dims = 1)
+    end
+end
+```
+
+In this example, `saved_files` lists the files saved by each "process":
+
+```julia
+julia> saved_files
+2-element Vector{Vector{String}}:
+ ["simulation.pvtu", "simulation/simulation_1.vtu"]
+ ["simulation/simulation_2.vtu"]
+```
+
+Note that the files containing the actual data (in this case `simulation_*.vtu`) are stored in a separate `simulation` directory.
+
+## Example 2: Structured files
+
+This generates 4 serial [image data](@ref Image-data) files (`.vti`) and
+a header file (`.pvti`) combining them:
+
+```julia
+# Global grid
+xs_global = range(0, 2; length = 15)
+ys_global = range(-1, 1; length = 12)
+zs_global = range(0, 1; length = 4)
+
+extents = [
+    ( 1:10,  1:5, 1:4),  # process 1
+    (10:15,  1:5, 1:4),  # process 2
+    ( 1:10, 5:12, 1:4),  # process 3
+    (10:15, 5:12, 1:4),  # process 4
+]
 
-saved_files = pvtk_grid("simulation", x, y, cells; part = 1, nparts = 1) do pvtk
-    pvtk["Pressure"] = x
-    pvtk["Processor"] = rand(2)
+saved_files = Vector{Vector{String}}(undef, 4)  # files saved by each "process"
+
+for part = 1:4
+    is, js, ks = extents[part]  # local indices
+    xs, ys, zs = xs_global[is], ys_global[js], zs_global[ks]  # local grid
+    saved_files[part] = pvtk_grid(
+            "fields", xs, ys, zs;
+            part = part, extents = extents,
+        ) do pvtk
+        pvtk["Temperature"] = [x + 2y + 3z for x ∈ xs, y ∈ ys, z ∈ zs]
+    end
 end
 ```
 
-In this example, `saved_files` lists the two files saved at the end of the do-block:
+As in the previous example, `saved_files` lists the files saved by each "process":
 
 ```julia
 julia> saved_files
-2-element Vector{String}:
- "simulation.pvtu"
- "simulation/simulation_1.vtu"
+4-element Vector{Vector{String}}:
+ ["fields.pvti", "fields/fields_1.vti"]
+ ["fields/fields_2.vti"]
+ ["fields/fields_3.vti"]
+ ["fields/fields_4.vti"]
 ```
 
-Note that the files containing the actual data (in this case `simulation_1.vtu`) are stored in a separate `simulation` directory.
 
 ## Acknowledgements
 
 Thanks to [Francesc Verdugo](https://www.francescverdugo.com/) and [Alberto
 F. Martin](https://research.monash.edu/en/persons/alberto-f-martin) for
-the initial parallel file format implementation.
+the initial parallel file format implementation, and to [Corentin Lothode](https://lmrs.univ-rouen.fr/fr/persopage/corentin-lothode) for the initial work on structured grids.
diff --git a/src/gridtypes/pvtk_grid.jl b/src/gridtypes/pvtk_grid.jl
index fdaa325e..9c747e78 100644
--- a/src/gridtypes/pvtk_grid.jl
+++ b/src/gridtypes/pvtk_grid.jl
@@ -40,28 +40,33 @@ end
 compute_whole_extent(::Nothing) = nothing
 
 """
-    pvtk_grid(args...; part, nparts, extents, ismain = (part == 1), ghost_level = 0, kwargs...)
+    pvtk_grid(
+        filename, args...;
+        part, nparts, extents, ismain = (part == 1), ghost_level = 0,
+        kwargs...,
+    )
 
-Returns a handler representing a parallel vtk file, which can be
+Returns a handler representing a parallel VTK file, which can be
 eventually written to file with `vtk_save`.
 
-Positional and keyword arguments in `args` and `kwargs`
-are passed to `vtk_grid` verbatim (except file names that are augmented with the
-corresponding part id).
+Positional and keyword arguments in `args` and `kwargs` are passed to `vtk_grid`
+verbatim.
+Note that serial filenames are automatically generated from `filename` and from
+the process id `part`.
 
-The extra keyword arguments only apply to parallel VTK file formats.
+The following keyword arguments only apply to parallel VTK file formats.
 
 Mandatory ones are:
 
 - `part`: current (1-based) part id,
-- `nparts`: total number of parts (only **unstructured** grids),
+- `nparts`: total number of parts (only required for **unstructured** grids),
 - `extents`: array specifying the partitioning of a **structured** grid across
   different processes (see below for details).
 
 Optional ones are:
 
-- `ismain` true if the current part id `part` is the main (the only one that will write the `.pvtk` file),
-- `ghost_level` ghost level.
+- `ismain`: `true` if the current part id `part` is the main (the only one that will write the `.pvtk` file),
+- `ghost_level`: ghost level.
 
 ## Specifying extents for a structured grid
 
@@ -70,7 +75,7 @@ must be specified via the `extents` argument.
 This is an array where each element represents the data extent associated to a
 given process.
 
-For instance, for a dataset of global dimensions `15×12×4` distributed across 4
+For example, for a dataset of global dimensions ``15×12×4`` distributed across 4
 processes, this array may look like the following:
 
 ```julia
@@ -86,7 +91,9 @@ Some important notes:
 
 - the `extents` argument must be the same for all processes;
 - the extents **must overlap**, or VTK / ParaView will complain when trying to
-  open the files.
+  open the files;
+- the length of the `extents` array gives the number of processes.
+  Therefore, the `nparts` argument is redundant and does not need to be passed.
 
 """
 function pvtk_grid(