diff --git a/Manifest.toml b/Manifest.toml
index b40d3fd7..bb44ef8e 100644
--- a/Manifest.toml
+++ b/Manifest.toml
@@ -148,7 +148,7 @@ version = "0.10.20"
 
 [[Gridap]]
 deps = ["AbstractTrees", "BSON", "BlockArrays", "Combinatorics", "DocStringExtensions", "FastGaussQuadrature", "FileIO", "FillArrays", "ForwardDiff", "JLD2", "JSON", "LineSearches", "LinearAlgebra", "NLsolve", "NearestNeighbors", "QuadGK", "Random", "SparseArrays", "SparseMatricesCSR", "StaticArrays", "Test", "WriteVTK"]
-git-tree-sha1 = "192cb4f55c2eeefcd7e39dd208a58e46fc9840a5"
+git-tree-sha1 = "d3bfeb769259b08f7385bdad91681f38c50fc519"
 repo-rev = "gridap_distributed"
 repo-url = "https://github.com/gridap/Gridap.jl.git"
 uuid = "56d4f2e9-7ea1-5844-9cf6-b9c51ca7ce8e"
diff --git a/Project.toml b/Project.toml
index f4981849..27dddced 100644
--- a/Project.toml
+++ b/Project.toml
@@ -10,6 +10,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195"
 PartitionedArrays = "5a9dfac6-5c52-46f7-8278-5e2210713be9"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
 
 [compat]
 Gridap = "0.17"
diff --git a/src/GridapDistributed.jl b/src/GridapDistributed.jl
index 2cf7e16f..65eff54c 100644
--- a/src/GridapDistributed.jl
+++ b/src/GridapDistributed.jl
@@ -15,6 +15,7 @@ using PartitionedArrays
 const PArrays = PartitionedArrays
 
 using SparseArrays
+using WriteVTK
 
 import Gridap.TensorValues: inner, outer, double_contraction, symmetric_part
 import LinearAlgebra: det, tr, cross, dot, ⋅
diff --git a/src/Visualization.jl b/src/Visualization.jl
index bf2b2e62..ca5dc3a0 100644
--- a/src/Visualization.jl
+++ b/src/Visualization.jl
@@ -6,7 +6,7 @@ function Base.getproperty(x::DistributedVisualizationData, sym::Symbol)
   if sym == :grid
     map_parts(i->i.grid,x.visdata)
   elseif sym == :filebase
-    map_parts(i->i.filebase,x.visdata)
+    get_part(x.visdata).filebase
   elseif sym == :celldata
     map_parts(i->i.celldata,x.visdata)
   elseif sym == :nodaldata
@@ -32,8 +32,7 @@ function Visualization.visualization_data(
   vd = map_parts(
     parts,model.models,model.gids.partition,labels.labels) do part,model,gids,labels
 
-    n = lpad(part,ceil(Int,log10(nparts)),'0')
-    vd = visualization_data(model,"$(filebase)_$(n)";labels=labels)
+    vd = visualization_data(model,filebase;labels=labels)
     vd_cells = vd[end]
     push!(vd_cells.celldata, "gid" => gids.lid_to_gid)
     push!(vd_cells.celldata, "part" => gids.lid_to_part)
@@ -63,11 +62,16 @@ function Visualization.visualization_data(
 
   vd = map_parts(
     parts,trians,cdat,fdat) do part,trian,celldata,cellfields
-    n = lpad(part,ceil(Int,log10(nparts)),'0')
+    _celldata = Dict{Any,Any}(celldata)
+    # we do not use "part" since it is likely to be used by the user
+    if haskey(_celldata,"piece")
+      @unreachable "piece is a reserved cell data name"
+    end
+    _celldata["piece"] = fill(part,num_cells(trian))
     vd = visualization_data(
-      trian,"$(filebase)_$(n)";
+      trian,filebase;
       order=order,nsubcells=nsubcells,
-      celldata=celldata,cellfields=cellfields)
+      celldata=_celldata,cellfields=cellfields)
     @assert length(vd) == 1
     vd[1]
   end
@@ -131,19 +135,21 @@ end
 
 # Vtk related
 
-# TODO use pvd format
-
 function Visualization.write_vtk_file(
   grid::AbstractPData{<:Grid}, filebase; celldata, nodaldata)
-  map_parts(grid,filebase,celldata,nodaldata) do g,f,c,n
-    write_vtk_file(g,f;celldata=c,nodaldata=n)
-  end
+  pvtk = Visualization.create_vtk_file(grid,filebase;celldata=celldata,nodaldata=nodaldata)
+  map_parts(vtk_save,pvtk)
 end
 
 function Visualization.create_vtk_file(
   grid::AbstractPData{<:Grid}, filebase; celldata, nodaldata)
-  map_parts(grid,filebase,celldata,nodaldata) do g,f,c,n
-    create_vtk_file(g,f;celldata=c,nodaldata=n)
+  parts = get_part_ids(grid)
+  nparts = length(parts)
+  map_parts(parts,grid,celldata,nodaldata) do part,g,c,n
+    Visualization.create_pvtk_file(
+      g,filebase;
+      pvtkargs=[:part=>part,:nparts=>nparts],
+      celldata=c,nodaldata=n)
   end
 end