diff --git a/CHANGELOG.md b/CHANGELOG.md
index 69acecbf52..dfd8bacd81 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -190,6 +190,11 @@ more discussion).
   + add!(dh, :u, Lagrange{RefTriangle, 1}())
   ```
 
+- **VTK export**: Ferrite no longer extends methods from `WriteVTK.jl`, instead the new types 
+  `VTKFile` and `VTKFileCollection` should be used instead. New methods exists for writing to 
+  a `VTKFile`, e.g. `write_solution`, `write_cell_data`, `write_node_data`, and `write_projection`.
+  See [#692][github-692].
+
 ### Added
 
 - `InterfaceValues` for computing jumps and averages over interfaces. ([#743][github-743])
@@ -858,6 +863,7 @@ poking into Ferrite internals:
 [github-684]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/684
 [github-687]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/687
 [github-688]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/688
+[github-692]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/692
 [github-694]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/694
 [github-695]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/695
 [github-697]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/697
@@ -896,4 +902,3 @@ poking into Ferrite internals:
 [github-835]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/835
 [github-855]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/855
 [github-880]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/880
-
diff --git a/docs/src/literate-gallery/helmholtz.jl b/docs/src/literate-gallery/helmholtz.jl
index 89795ed995..52c97c34b7 100644
--- a/docs/src/literate-gallery/helmholtz.jl
+++ b/docs/src/literate-gallery/helmholtz.jl
@@ -162,9 +162,9 @@ K, f = doassemble(cellvalues, facetvalues, K, dh);
 apply!(K, f, dbcs)
 u = Symmetric(K) \ f;
 
-vtkfile = vtk_grid("helmholtz", dh)
-vtk_point_data(vtkfile, dh, u)
-vtk_save(vtkfile)
+vtk = VTKFile("helmholtz", dh)
+write_solution(vtk, dh, u)
+close(vtk)
 using Test #src
 #src this test catches unexpected changes in the result over time.
 #src the true maximum is slightly bigger then 1.0
diff --git a/docs/src/literate-gallery/landau.jl b/docs/src/literate-gallery/landau.jl
index 7eb3fb9108..b84283c869 100644
--- a/docs/src/literate-gallery/landau.jl
+++ b/docs/src/literate-gallery/landau.jl
@@ -113,10 +113,10 @@ function LandauModel(α, G, gridsize, left::Vec{DIM, T}, right::Vec{DIM, T}, elp
 end
 
 # utility to quickly save a model
-function Ferrite.vtk_save(path, model, dofs=model.dofs)
-    vtkfile = vtk_grid(path, model.dofhandler)
-    vtk_point_data(vtkfile, model.dofhandler, dofs)
-    vtk_save(vtkfile)
+function save_landau(path, model, dofs=model.dofs)
+    VTKFile(path, model.dofhandler) do vtk
+        write_solution(vtk, model.dofhandler, dofs)
+    end
 end
 
 # ## Assembly
@@ -253,9 +253,9 @@ left = Vec{3}((-75.,-25.,-2.))
 right = Vec{3}((75.,25.,2.))
 model = LandauModel(α, G, (50, 50, 2), left, right, element_potential)
 
-vtk_save("landauorig", model)
+save_landau("landauorig", model)
 @time minimize!(model)
-vtk_save("landaufinal", model)
+save_landau("landaufinal", model)
 
 # as we can see this runs very quickly even for relatively large gridsizes.
 # The key to get high performance like this is to minimize the allocations inside the threaded loops,
diff --git a/docs/src/literate-gallery/quasi_incompressible_hyperelasticity.jl b/docs/src/literate-gallery/quasi_incompressible_hyperelasticity.jl
index 9bdcf043fd..40e699a3cf 100644
--- a/docs/src/literate-gallery/quasi_incompressible_hyperelasticity.jl
+++ b/docs/src/literate-gallery/quasi_incompressible_hyperelasticity.jl
@@ -305,7 +305,7 @@ function solve(interpolation_u, interpolation_p)
     Δt = 0.1;
     NEWTON_TOL = 1e-8
 
-    pvd = paraview_collection("hyperelasticity_incomp_mixed.pvd");
+    pvd = VTKFileCollection("hyperelasticity_incomp_mixed", grid);
     for t ∈ 0.0:Δt:Tf
         ## Perform Newton iterations
         Ferrite.update!(dbc, t)
@@ -334,13 +334,11 @@ function solve(interpolation_u, interpolation_p)
         end;
 
         ## Save the solution fields
-        vtk_grid("hyperelasticity_incomp_mixed_$t.vtu", dh) do vtkfile
-            vtk_point_data(vtkfile, dh, w)
-            vtk_save(vtkfile)
-            pvd[t] = vtkfile
+        addstep!(pvd, t) do io 
+            write_solution(io, dh, w)
         end
     end;
-    vtk_save(pvd);
+    close(pvd);
     vol_def = calculate_volume_deformed_mesh(w, dh, cellvalues_u);
     print("Deformed volume is $vol_def")
     return vol_def;
diff --git a/docs/src/literate-gallery/topology_optimization.jl b/docs/src/literate-gallery/topology_optimization.jl
index e56309afa5..e3b2c154f2 100644
--- a/docs/src/literate-gallery/topology_optimization.jl
+++ b/docs/src/literate-gallery/topology_optimization.jl
@@ -492,16 +492,16 @@ function topopt(ra,ρ,n,filename; output=:false)
             i = @sprintf("%3.3i", it)
             filename_it = string(filename, "_", i)
 
-            vtk_grid(filename_it, grid) do vtk
-                vtk_cell_data(vtk, χ, "density")
+            VTKFile(filename_it, grid) do vtk
+                write_cell_data(vtk, χ, "density")
             end
         end
     end
 
     ## export converged results
     if(!output)
-        vtk_grid(filename, grid) do vtk
-            vtk_cell_data(vtk, χ, "density")
+        VTKFile(filename, grid) do vtk
+            write_cell_data(vtk, χ, "density")
         end
     end
     @printf "Rel. stiffness: %.4f \n" compliance^(-1)/compliance_0^(-1)
diff --git a/docs/src/literate-howto/postprocessing.jl b/docs/src/literate-howto/postprocessing.jl
index 79586f0640..ef2e38818a 100644
--- a/docs/src/literate-howto/postprocessing.jl
+++ b/docs/src/literate-howto/postprocessing.jl
@@ -87,8 +87,8 @@ q_projected = project(projector, q_gp, qr);
 # To visualize the heat flux, we export the projected field `q_projected`
 # to a VTK-file, which can be viewed in e.g. [ParaView](https://www.paraview.org/).
 # The result is also visualized in *Figure 1*.
-vtk_grid("heat_equation_flux", grid) do vtk
-    vtk_point_data(vtk, projector, q_projected, "q")
+VTKFile("heat_equation_flux", grid) do vtk
+    write_projection(vtk, projector, q_projected, "q")
 end;
 
 # ## Point Evaluation
diff --git a/docs/src/literate-howto/threaded_assembly.jl b/docs/src/literate-howto/threaded_assembly.jl
index 9c3466ec2e..2f21fb33b1 100644
--- a/docs/src/literate-howto/threaded_assembly.jl
+++ b/docs/src/literate-howto/threaded_assembly.jl
@@ -25,9 +25,9 @@ function create_example_2d_grid()
     grid = generate_grid(Quadrilateral, (10, 10), Vec{2}((0.0, 0.0)), Vec{2}((10.0, 10.0)))
     colors_workstream = create_coloring(grid; alg=ColoringAlgorithm.WorkStream)
     colors_greedy = create_coloring(grid; alg=ColoringAlgorithm.Greedy)
-    vtk_grid("colored", grid) do vtk
-        vtk_cell_data_colors(vtk, colors_workstream, "workstream-coloring")
-        vtk_cell_data_colors(vtk, colors_greedy, "greedy-coloring")
+    VTKFile("colored", grid) do vtk
+        Ferrite.write_cell_colors(vtk, grid, colors_workstream, "workstream-coloring")
+        Ferrite.write_cell_colors(vtk, grid, colors_greedy, "greedy-coloring")
     end
 end
 
diff --git a/docs/src/literate-tutorials/computational_homogenization.jl b/docs/src/literate-tutorials/computational_homogenization.jl
index 8e341de172..11a5da477c 100644
--- a/docs/src/literate-tutorials/computational_homogenization.jl
+++ b/docs/src/literate-tutorials/computational_homogenization.jl
@@ -533,16 +533,16 @@ round.(ev; digits=-8)
 
 uM = zeros(ndofs(dh))
 
-vtk_grid("homogenization", dh) do vtk
+VTKFile("homogenization", dh) do vtk
     for i in 1:3
         ## Compute macroscopic solution
         apply_analytical!(uM, dh, :u, x -> εᴹ[i] ⋅ x)
         ## Dirichlet
-        vtk_point_data(vtk, dh, uM + u.dirichlet[i], "_dirichlet_$i")
-        vtk_point_data(vtk, projector, σ.dirichlet[i], "σvM_dirichlet_$i")
+        write_solution(vtk, dh, uM + u.dirichlet[i], "_dirichlet_$i")
+        write_projection(vtk, projector, σ.dirichlet[i], "σvM_dirichlet_$i")
         ## Periodic
-        vtk_point_data(vtk, dh, uM + u.periodic[i], "_periodic_$i")
-        vtk_point_data(vtk, projector, σ.periodic[i], "σvM_periodic_$i")
+        write_solution(vtk, dh, uM + u.periodic[i], "_periodic_$i")
+        write_projection(vtk, projector, σ.periodic[i], "σvM_periodic_$i")
     end
 end;
 
diff --git a/docs/src/literate-tutorials/dg_heat_equation.jl b/docs/src/literate-tutorials/dg_heat_equation.jl
index a44a9f02f0..2e9d39f692 100644
--- a/docs/src/literate-tutorials/dg_heat_equation.jl
+++ b/docs/src/literate-tutorials/dg_heat_equation.jl
@@ -334,8 +334,8 @@ K, f = assemble_global(cellvalues, facetvalues, interfacevalues, K, dh, order, d
 
 apply!(K, f, ch)
 u = K \ f;
-vtk_grid("dg_heat_equation", dh) do vtk
-    vtk_point_data(vtk, dh, u)
+VTKFile("dg_heat_equation", dh) do vtk
+    write_solution(vtk, dh, u)
 end;
 
 ## test the result                #src
diff --git a/docs/src/literate-tutorials/heat_equation.jl b/docs/src/literate-tutorials/heat_equation.jl
index fea8386f79..b924ce653a 100644
--- a/docs/src/literate-tutorials/heat_equation.jl
+++ b/docs/src/literate-tutorials/heat_equation.jl
@@ -213,8 +213,8 @@ u = K \ f;
 # ### Exporting to VTK
 # To visualize the result we export the grid and our field `u`
 # to a VTK-file, which can be viewed in e.g. [ParaView](https://www.paraview.org/).
-vtk_grid("heat_equation", dh) do vtk
-    vtk_point_data(vtk, dh, u)
+VTKFile("heat_equation", dh) do vtk
+    write_solution(vtk, dh, u)
 end
 
 ## test the result                #src
diff --git a/docs/src/literate-tutorials/hyperelasticity.jl b/docs/src/literate-tutorials/hyperelasticity.jl
index 7007f91cb9..9db1554cb7 100644
--- a/docs/src/literate-tutorials/hyperelasticity.jl
+++ b/docs/src/literate-tutorials/hyperelasticity.jl
@@ -408,8 +408,8 @@ function solve()
 
     ## Save the solution
     @timeit "export" begin
-        vtk_grid("hyperelasticity", dh) do vtkfile
-            vtk_point_data(vtkfile, dh, u)
+        VTKFile("hyperelasticity", dh) do vtk
+            write_solution(vtk, dh, u)
         end
     end
 
diff --git a/docs/src/literate-tutorials/incompressible_elasticity.jl b/docs/src/literate-tutorials/incompressible_elasticity.jl
index 36cfcc661f..14a1e0c69a 100644
--- a/docs/src/literate-tutorials/incompressible_elasticity.jl
+++ b/docs/src/literate-tutorials/incompressible_elasticity.jl
@@ -270,13 +270,14 @@ function solve(ν, interpolation_u, interpolation_p)
     ## Export the solution and the stress
     filename = "cook_" * (interpolation_u == Lagrange{RefTriangle, 1}()^2 ? "linear" : "quadratic") *
                          "_linear"
-    vtk_grid(filename, dh) do vtkfile
-        vtk_point_data(vtkfile, dh, u)
+
+    VTKFile(filename, grid) do vtk
+        write_solution(vtk, dh, u)
         for i in 1:3, j in 1:3
             σij = [x[i, j] for x in σ]
-            vtk_cell_data(vtkfile, σij, "sigma_$(i)$(j)")
+            write_cell_data(vtk, σij, "sigma_$(i)$(j)")
         end
-        vtk_cell_data(vtkfile, σvM, "sigma von Mise")
+        write_cell_data(vtk, σvM, "sigma von Mises")
     end
     return u
 end
diff --git a/docs/src/literate-tutorials/linear_shell.jl b/docs/src/literate-tutorials/linear_shell.jl
index 07890e18e6..54ef96bbb9 100644
--- a/docs/src/literate-tutorials/linear_shell.jl
+++ b/docs/src/literate-tutorials/linear_shell.jl
@@ -113,8 +113,8 @@ a = K\f
 
 # Output results.
 #+
-vtk_grid("linear_shell", dh) do vtk
-    vtk_point_data(vtk, dh, a)
+VTKFile("linear_shell", dh) do vtk
+    write_solution(vtk, dh, a)
 end
 
 end; #end main functions
diff --git a/docs/src/literate-tutorials/ns_vs_diffeq.jl b/docs/src/literate-tutorials/ns_vs_diffeq.jl
index 7139096447..352af3ffcc 100644
--- a/docs/src/literate-tutorials/ns_vs_diffeq.jl
+++ b/docs/src/literate-tutorials/ns_vs_diffeq.jl
@@ -433,7 +433,7 @@ integrator = init(
     progress=true, progress_steps=1,
     saveat=Δt_save);
 
-pvd = paraview_collection("vortex-street.pvd");
+pvd = VTKFileCollection("vortex-street", grid);
 integrator = TimeChoiceIterator(integrator, 0.0:Δt_save:T)
 for (u_uc,t) in integrator
     # We ignored the Dirichlet constraints in the solution vector up to now,
@@ -442,13 +442,11 @@ for (u_uc,t) in integrator
     update!(ch, t)
     u = copy(u_uc)
     apply!(u, ch)
-    vtk_grid("vortex-street-$t.vtu", dh) do vtk
-        vtk_point_data(vtk,dh,u)
-        vtk_save(vtk)
-        pvd[t] = vtk
+    addstep!(pvd, t) do io
+        write_solution(io, dh, u)
     end
 end
-vtk_save(pvd);
+close(pvd);
 
 # Test the result for full proper development of the flow                   #src
 using Test                                                                  #hide
diff --git a/docs/src/literate-tutorials/plasticity.jl b/docs/src/literate-tutorials/plasticity.jl
index 33bd8fb5b8..6b34be0e2e 100644
--- a/docs/src/literate-tutorials/plasticity.jl
+++ b/docs/src/literate-tutorials/plasticity.jl
@@ -344,10 +344,10 @@ function solve()
         mises_values[el] /= length(cell_states) # average von Mises stress
         κ_values[el] /= length(cell_states)     # average drag stress
     end
-    vtk_grid("plasticity", dh) do vtkfile
-        vtk_point_data(vtkfile, dh, u) # displacement field
-        vtk_cell_data(vtkfile, mises_values, "von Mises [Pa]")
-        vtk_cell_data(vtkfile, κ_values, "Drag stress [Pa]")
+    VTKFile("plasticity", dh) do vtk
+        write_solution(vtk, dh, u) # displacement field
+        write_cell_data(vtk, mises_values, "von Mises [Pa]")
+        write_cell_data(vtk, κ_values, "Drag stress [Pa]")
     end
 
     return u_max, traction_magnitude
diff --git a/docs/src/literate-tutorials/porous_media.jl b/docs/src/literate-tutorials/porous_media.jl
index fb2a3c23df..335a1bef22 100644
--- a/docs/src/literate-tutorials/porous_media.jl
+++ b/docs/src/literate-tutorials/porous_media.jl
@@ -352,8 +352,8 @@ function solve(dh, ch, domains; Δt=0.025, t_total=1.0)
     r = zeros(ndofs(dh))
     a = zeros(ndofs(dh))
     a_old = copy(a)
-    pvd = paraview_collection("porous_media.pvd");
-    for (step, t) in enumerate(0:Δt:t_total)
+    pvd = VTKFileCollection("porous_media.pvd", dh);
+    for t in 0:Δt:t_total
         if t>0
             update!(ch, t)
             apply!(a, ch)    
@@ -364,13 +364,11 @@ function solve(dh, ch, domains; Δt=0.025, t_total=1.0)
             a .+= Δa
             copyto!(a_old, a)
         end
-        vtk_grid("porous_media-$step", dh) do vtk
-            vtk_point_data(vtk, dh, a)
-            vtk_save(vtk)
-            pvd[step] = vtk
+        addstep!(pvd, t) do io
+            write_solution(io, dh, a)
         end
     end
-    vtk_save(pvd);
+    close(pvd);
 end;
 
 # Finally we call the functions to actually run the code
diff --git a/docs/src/literate-tutorials/stokes-flow.jl b/docs/src/literate-tutorials/stokes-flow.jl
index a770a91a14..dfa0854c3a 100644
--- a/docs/src/literate-tutorials/stokes-flow.jl
+++ b/docs/src/literate-tutorials/stokes-flow.jl
@@ -528,8 +528,8 @@ function main()
     u = K \ f
     apply!(u, ch)
     ## Export the solution
-    vtk_grid("stokes-flow", grid) do vtk
-        vtk_point_data(vtk, dh, u)
+    VTKFile("stokes-flow", grid) do vtk
+        write_solution(vtk, dh, u)
     end
 
     ## Check the result                #src
diff --git a/docs/src/literate-tutorials/transient_heat_equation.jl b/docs/src/literate-tutorials/transient_heat_equation.jl
index 763e729a8d..2a865bed08 100644
--- a/docs/src/literate-tutorials/transient_heat_equation.jl
+++ b/docs/src/literate-tutorials/transient_heat_equation.jl
@@ -189,13 +189,11 @@ apply_analytical!(uₙ, dh, :u, x -> (x[1]^2 - 1) * (x[2]^2 - 1) * max_temp);
 # Here, we apply **once** the boundary conditions to the system matrix `A`.
 apply!(A, ch);
 
-# To store the solution, we initialize a `paraview_collection` (.pvd) file.
-pvd = paraview_collection("transient-heat.pvd");
+# To store the solution, we initialize a `VTKFileCollection` (.pvd) file.
+pvd = VTKFileCollection("transient-heat", grid);
 t = 0
-vtk_grid("transient-heat-$t", dh) do vtk
-    vtk_point_data(vtk, dh, uₙ)
-    vtk_save(vtk)
-    pvd[t] = vtk
+addstep!(pvd, t) do io
+    write_solution(io, dh, uₙ)
 end
 
 # At this point everything is set up and we can finally approach the time loop.
@@ -211,16 +209,14 @@ for t in Δt:Δt:T
     #Finally, we can solve the time step and save the solution afterwards.
     u = A \ b
 
-    vtk_grid("transient-heat-$t", dh) do vtk
-        vtk_point_data(vtk, dh, u)
-        vtk_save(vtk)
-        pvd[t] = vtk
+    addstep!(pvd, t) do io
+        write_solution(io, dh, u)
     end
     #At the end of the time loop, we set the previous solution to the current one and go to the next time step.
     uₙ .= u
 end
 # In order to use the .pvd file we need to store it to the disk, which is done by:
-vtk_save(pvd);
+close(pvd);
 
 #md # ## [Plain program](@id transient_heat_equation-plain-program)
 #md #
diff --git a/docs/src/reference/export.md b/docs/src/reference/export.md
index 2fc0b26a53..81a77e1cc9 100644
--- a/docs/src/reference/export.md
+++ b/docs/src/reference/export.md
@@ -24,10 +24,16 @@ evaluate_at_grid_nodes
 ```
 
 ## VTK Export
-
 ```@docs
-vtk_grid(filename::AbstractString, grid::Grid{dim,C,T}; compress::Bool) where {dim,C,T} 
-vtk_point_data
-vtk_cellset
-vtk_cell_data_colors
+VTKFile
+VTKFileCollection
+addstep!
+write_solution
+write_projection
+write_cell_data
+write_node_data
+Ferrite.write_cellset
+Ferrite.write_nodeset
+Ferrite.write_constraints
+Ferrite.write_cell_colors
 ```
diff --git a/docs/src/topics/export.md b/docs/src/topics/export.md
index 855b5d62c8..702561f522 100644
--- a/docs/src/topics/export.md
+++ b/docs/src/topics/export.md
@@ -1,7 +1,7 @@
 ```@setup export
 using Ferrite
 grid = generate_grid(Triangle, (2, 2))
-dh = DofHandler(grid); add!(dh, :u, Lagrange{2,RefTetrahedron,1}()); close!(dh)
+dh = DofHandler(grid); add!(dh, :u, Lagrange{RefTriangle,1}()); close!(dh)
 u = rand(ndofs(dh)); σ = rand(getncells(grid))
 ```
 
@@ -10,98 +10,50 @@ u = rand(ndofs(dh)); σ = rand(getncells(grid))
 When the problem is solved, and the solution vector `u` is known we typically
 want to visualize it. The simplest way to do this is to write the solution to a
 VTK-file, which can be viewed in e.g. [`Paraview`](https://www.paraview.org/).
-To write VTK-files, Ferrite uses, and extends, functions from the
-[`WriteVTK.jl`](https://github.com/jipolanco/WriteVTK.jl) package to simplify
+To write VTK-files, Ferrite comes with an export interface with a 
+[`WriteVTK.jl`](https://github.com/jipolanco/WriteVTK.jl) backend to simplify
 the exporting.
 
-First we need to create a file, based on the grid. This is done with the
-`vtk_grid` function:
-
-```@example export
-vtk = vtk_grid("my-solution", grid)
-# hide
-```
-
-Next we have to add data to the file. We may add different kinds of data;
-point data using `vtk_point_data` or cell data using
-`vtk_cell_data`. Point data is data for each nodal coordinate in the
-grid, for example our solution vector. Point data can be either scalars
-or vectors. Cell data is -- as the name suggests -- data for each cell. This
-can be for example the stress. As an example, lets add a solution vector `u`
-as point data, and a vector with stress for each cell, `σ`, as cell data:
-
-```@example export
-vtk_point_data(vtk, u, "my-point-data")
-vtk_cell_data(vtk,  σ, "my-cell-data")
-# hide
-```
-
-Finally, we need to save the file to disk, using `vtk_save`
-
-```@example export
-vtk_save(vtk)
-rm("my-solution.vtu") # hide
-```
-
-Alternatively, all of the above can be done using a `do` block:
-
-```@example export
-vtk_grid("my-solution", grid) do vtk
-    vtk_point_data(vtk, u, "my-point-data")
-    vtk_cell_data(vtk, σ, "my-cell-data")
-end
-rm("my-solution.vtu") # hide
-```
-
-For other functionality, and more information refer to the
-[`WriteVTK.jl` README](https://github.com/jipolanco/WriteVTK.jl/blob/master/README.md).
-In particular, for exporting the solution at multiple time steps, the
-[section on PVD files](https://github.com/jipolanco/WriteVTK.jl#paraview-data-pvd-file-format)
-is useful.
-
-## Exporting with `DofHandler`
-
-There is an even more convenient way to export a solution vector `u` -- using the
-`DofHandler`. The `DofHandler` already contains all of the information needed,
-such as the names of our fields and if they are scalar or vector fields. But most
-importantly the `DofHandler` knows about the numbering and distribution of
-degrees of freedom, and thus knows how to "distribute" the solution vector on
-the grid. For example, lets say we have a `DofHandler` `dh` and a solution
-vector `u`:
-
+The following structure can be used to write various output to a vtk-file:
 ```@example export
-vtk = vtk_grid("my-solution", dh)
-vtk_point_data(vtk, dh, u)
-vtk_save(vtk)
-rm("my-solution.vtu") # hide
-```
-
-or with a `do`-block:
-
+VTKFile("my_solution", grid) do vtk
+    write_solution(vtk, dh, u)
+end;
+```
+where `write_solution` is just one example of the following functions that can be used 
+
+* [`write_solution`](@ref)
+* [`write_cell_data`](@ref)
+* [`write_node_data`](@ref)
+* [`write_projection`](@ref)
+* [`Ferrite.write_cellset`](@ref)
+* [`Ferrite.write_nodeset`](@ref)
+* [`Ferrite.write_constraints`](@ref)
+* [`Ferrite.write_cell_colors`](@ref)
+
+Instead of using the `do`-block, it is also possible to do
 ```@example export
-vtk_grid("my-solution", dh) do vtk
-    vtk_point_data(vtk, dh, u)
-    vtk_cell_data(vtk, σ, "my-cell-data")
+vtk = VTKFile("my_solution", grid)
+write_solution(vtk, dh, u)
+# etc.
+close(vtk);
+```
+
+The data written by `write_solution`, `write_cell_data`, `write_node_data`, and `write_projection` may be either scalar (`Vector{<:Number}`) or tensor (`Vector{<:AbstractTensor}`) data. 
+
+For simulations with multiple time steps, typically one `VTK` (`.vtu`) file is written 
+for each time step. In order to connect the actual time with each of these files,
+a `VTKFileCollection` can be used, which will write one paraview datafile (`.pvd`)
+file and one `VTKFile` (`.vtu`) for each time step. 
+
+```@example export 
+pvd = VTKFileCollection("my_results", grid)
+for t in range(0, 1, 5)
+    # Do calculations to update u
+    addstep!(pvd, t) do vtk
+        write_solution(vtk, dh, u)
+    end
 end
-rm("my-solution.vtu") # hide
+close(pvd);
 ```
-
-When `vtk_point_data` is used with a `DofHandler` all of the fields will be
-written to the VTK file, and the names will be determined by the fieldname
-symbol that was used when the field was added to the `DofHandler`.
-
-## Exporting Boundary Conditions
-
-There is also a `vtk_point_data` which accepts a `ConstraintHandler`.
-This method is useful to verify that the boundary conditions are
-applied where they are supposed to. For a `ConstraintHandler` `ch`
-we can export the boundary conditions as
-
-```julia
-vtk_grid("boundary-conditions", grid) do vtk
-    vtk_point_data(vtk, ch)
-end
-```
-
-This will export zero-valued fields with ones on the parts where the
-boundary conditions are active.
+See [Transient heat equation](@ref tutorial-transient-heat-equation) for an example
\ No newline at end of file
diff --git a/src/Dofs/ConstraintHandler.jl b/src/Dofs/ConstraintHandler.jl
index d398bfeae0..c8529bc8f2 100644
--- a/src/Dofs/ConstraintHandler.jl
+++ b/src/Dofs/ConstraintHandler.jl
@@ -486,42 +486,6 @@ function _update!(inhomogeneities::Vector{T}, f::Function, ::Set{Int}, field::Sy
     end
 end
 
-# Saves the dirichlet boundary conditions to a vtkfile.
-# Values will have a 1 where bcs are active and 0 otherwise
-function WriteVTK.vtk_point_data(vtkfile, ch::ConstraintHandler)
-    unique_fields = []
-    for dbc in ch.dbcs
-        push!(unique_fields, dbc.field_name)
-    end
-    unique!(unique_fields)
-
-    for field in unique_fields
-        nd = getfielddim(ch.dh, field)
-        data = zeros(Float64, nd, getnnodes(get_grid(ch.dh)))
-        for dbc in ch.dbcs
-            dbc.field_name != field && continue
-            if eltype(dbc.facets) <: BoundaryIndex
-                functype = boundaryfunction(eltype(dbc.facets))
-                for (cellidx, facetidx) in dbc.facets
-                    for facetnode in functype(getcells(get_grid(ch.dh), cellidx))[facetidx]
-                        for component in dbc.components
-                            data[component, facetnode] = 1
-                        end
-                    end
-                end
-            else
-                for nodeidx in dbc.facets
-                    for component in dbc.components
-                        data[component, nodeidx] = 1
-                    end
-                end
-            end
-        end
-        vtk_point_data(vtkfile, data, string(field, "_bc"))
-    end
-    return vtkfile
-end
-
 """
 
     apply!(K::SparseMatrixCSC, rhs::AbstractVector, ch::ConstraintHandler)
diff --git a/src/Export/VTK.jl b/src/Export/VTK.jl
index 1b537f9da5..5df796035f 100644
--- a/src/Export/VTK.jl
+++ b/src/Export/VTK.jl
@@ -1,3 +1,135 @@
+
+"""
+    VTKFile(filename::AbstractString, grid::AbstractGrid; kwargs...)
+    VTKFile(filename::AbstractString, dh::DofHandler; kwargs...)
+
+Create a `VTKFile` that contains an unstructured VTK grid. 
+The keyword arguments are forwarded to `WriteVTK.vtk_grid`, see 
+[Data Formatting Options](https://juliavtk.github.io/WriteVTK.jl/stable/grids/syntax/#Data-formatting-options)
+
+This file handler can be used to to write data with
+
+* [`write_solution`](@ref)
+* [`write_cell_data`](@ref)
+* [`write_projection`](@ref)
+* [`write_node_data`](@ref).
+* [`Ferrite.write_cellset`](@ref)
+* [`Ferrite.write_nodeset`](@ref)
+* [`Ferrite.write_constraints`](@ref)
+
+It is necessary to call `close(::VTKFile)` to save the data after writing 
+to the file handler. Using the supported `do`-block does this automatically:
+```julia
+VTKFile(filename, grid) do vtk
+    write_solution(vtk, dh, u)
+    write_cell_data(vtk, celldata)
+end
+```
+"""
+struct VTKFile{VTK<:WriteVTK.VTKFile}
+    vtk::VTK
+end
+function VTKFile(filename::String, dh::DofHandler; kwargs...)
+    return VTKFile(filename, get_grid(dh); kwargs...)
+end
+function VTKFile(filename::String, grid::AbstractGrid; kwargs...)
+    vtk = create_vtk_grid(filename, grid; kwargs...)
+    return VTKFile(vtk)
+end
+# Makes it possible to use the `do`-block syntax
+function VTKFile(f::Function, args...; kwargs...)
+    vtk = VTKFile(args...; kwargs...)
+    try
+        f(vtk)
+    finally
+        close(vtk)
+    end
+end
+
+Base.close(vtk::VTKFile) = WriteVTK.vtk_save(vtk.vtk)
+
+function Base.show(io::IO, ::MIME"text/plain", vtk::VTKFile)
+    open_str = WriteVTK.isopen(vtk.vtk) ? "open" : "closed"
+    filename = vtk.vtk.path
+    print(io, "VTKFile for the $open_str file \"$(filename)\".")
+end
+
+"""
+    VTKFileCollection(name::String, grid::AbstractGrid; kwargs...)
+    VTKFileCollection(name::String, dh::DofHandler; kwargs...)
+
+Create a paraview data file (.pvd) that can be used to 
+save multiple vtk file along with a time stamp. The keyword arguments 
+are forwarded to `WriteVTK.paraview_collection`.
+
+See [`addstep!`](@ref) for examples for how to use `VTKFileCollection`.
+```
+"""
+mutable struct VTKFileCollection{P<:WriteVTK.CollectionFile,G_DH}
+    pvd::P
+    grid_or_dh::G_DH
+    name::String
+    step::Int
+end
+function VTKFileCollection(name::String, grid_or_dh::Union{AbstractGrid,AbstractDofHandler}; kwargs...)
+    pvd = WriteVTK.paraview_collection(name; kwargs...)
+    basename = string(first(split(pvd.path, ".pvd")))
+    return VTKFileCollection(pvd, grid_or_dh, basename, 0)
+end
+Base.close(pvd::VTKFileCollection) = WriteVTK.vtk_save(pvd.pvd)
+
+"""
+    addstep!(f::Function, pvd::VTKFileCollection, t::Real, [grid_or_dh]; kwargs...)
+
+Add a step at time `t` by writing a `VTKFile` to `pvd`. 
+The keyword arguments are forwarded to `WriteVTK.vtk_grid`.
+If required, a new grid can be used by supplying the grid or dofhandler as the last argument.
+Should be used in a do-block, e.g.
+```julia
+filename = "myoutput"
+pvd = VTKFileCollection(filename, grid)
+for (n, t) in pairs(timevector)
+    # Calculate, e.g., the solution `u` and the stress `σeff`
+    addstep!(pvd, t) do io 
+        write_cell_data(io, σeff, "Effective stress")
+        write_solution(io, dh, u)
+    end
+end
+close(pvd)
+```
+"""
+function addstep!(f::Function, pvd::VTKFileCollection, t, grid=pvd.grid_or_dh; kwargs...)
+    pvd.step += 1
+    VTKFile(string(pvd.name, "_", pvd.step), grid; kwargs...) do vtk
+        f(vtk)
+        pvd.pvd[t] = vtk.vtk # Add to collection
+    end
+end
+
+"""
+    addstep!(pvd::VTKFileCollection, vtk::VTKFile, t)
+
+As an alternative to using the `addstep!(pvd, t) do` block, it is 
+also possible to add a specific `vtk` at time `t` to `pvd`. 
+Note that this will close the `vtk`. Example
+```julia
+filename = "myoutput"
+pvd = VTKFileCollection(filename, grid)
+for (n, t) in pairs(timevector)
+    # Calculate, e.g., the solution `u` and the stress `σeff`
+    vtk = VTKFile(string(filename, "_", n), dh)
+    write_cell_data(vtk, σeff, "Effective stress")
+    write_solution(vtk, dh, u)
+    addstep!(pvd, vtk, t)
+end
+close(pvd)
+```
+"""
+function addstep!(pvd::VTKFileCollection, vtk::VTKFile, t)
+    pvd.step += 1
+    pvd.pvd[t] = vtk.vtk
+end
+
 cell_to_vtkcell(::Type{Line}) = VTKCellTypes.VTK_LINE
 cell_to_vtkcell(::Type{QuadraticLine}) = VTKCellTypes.VTK_QUADRATIC_EDGE
 
@@ -47,56 +179,46 @@ nodes_to_vtkorder(cell::QuadraticHexahedron) = [
     cell.nodes[27], # interior
 ]
 
-"""
-    vtk_grid(filename::AbstractString, grid::Grid; kwargs...)
-    vtk_grid(filename::AbstractString, dh::DofHandler; kwargs...)
-
-Create a unstructured VTK grid from `grid` (alternatively from the `grid` stored in `dh`). 
-Return a `DatasetFile` that data can be appended to, see 
-[`vtk_point_data`](@ref) and [`vtk_cell_data`](@ref).
-The keyword arguments are forwarded to `WriteVTK.vtk_grid`, see 
-[Data Formatting Options](https://juliavtk.github.io/WriteVTK.jl/stable/grids/syntax/#Data-formatting-options)
-"""
-function WriteVTK.vtk_grid(filename::AbstractString, grid::Grid{dim,C,T}; kwargs...) where {dim,C,T}
-    cls = MeshCell[]
+function create_vtk_griddata(grid::Grid{dim,C,T}) where {dim,C,T}
+    cls = WriteVTK.MeshCell[]
     for cell in getcells(grid)
         celltype = Ferrite.cell_to_vtkcell(typeof(cell))
-        push!(cls, MeshCell(celltype, nodes_to_vtkorder(cell)))
+        push!(cls, WriteVTK.MeshCell(celltype, nodes_to_vtkorder(cell)))
     end
     coords = reshape(reinterpret(T, getnodes(grid)), (dim, getnnodes(grid)))
-    return vtk_grid(filename, coords, cls; kwargs...)
+    return coords, cls
 end
-function WriteVTK.vtk_grid(filename::AbstractString, dh::AbstractDofHandler; kwargs...)
-    vtk_grid(filename, get_grid(dh); kwargs...)
+
+function create_vtk_grid(filename::AbstractString, grid::Grid{dim,C,T}; kwargs...) where {dim,C,T}
+    coords, cls = create_vtk_griddata(grid)
+    return WriteVTK.vtk_grid(filename, coords, cls; kwargs...)
 end
 
 function toparaview!(v, x::Vec{D}) where D
     v[1:D] .= x
 end
-function toparaview!(v, x::SecondOrderTensor{D}) where D
+function toparaview!(v, x::SecondOrderTensor)
     tovoigt!(v, x)
 end
 
-"""
-    vtk_point_data(vtk, data::Vector{<:AbstractTensor}, name)
-
-Write the tensor field `data` to the vtk file. Two-dimensional tensors are padded with zeros.
-
-For second order tensors the following indexing ordering is used:
-`[11, 22, 33, 23, 13, 12, 32, 31, 21]`. This is the default Voigt order in Tensors.jl.
-"""
-function WriteVTK.vtk_point_data(
+function _vtk_write_node_data(
     vtk::WriteVTK.DatasetFile,
-    data::Vector{S},
+    nodedata::Vector{S},
     name::AbstractString
     ) where {O, D, T, M, S <: Union{Tensor{O, D, T, M}, SymmetricTensor{O, D, T, M}}}
     noutputs = S <: Vec{2} ? 3 : M # Pad 2D Vec to 3D
-    npoints = length(data)
+    npoints = length(nodedata)
     out = zeros(T, noutputs, npoints)
     for i in 1:npoints
-        toparaview!(@view(out[:, i]), data[i])
+        toparaview!(@view(out[:, i]), nodedata[i])
     end
-    return vtk_point_data(vtk, out, name; component_names=component_names(S))
+    return WriteVTK.vtk_point_data(vtk, out, name; component_names=component_names(S))
+end
+function _vtk_write_node_data(vtk::WriteVTK.DatasetFile, nodedata::Vector{<:Real}, name::AbstractString)
+    return WriteVTK.vtk_point_data(vtk, nodedata, name)
+end
+function _vtk_write_node_data(vtk::WriteVTK.DatasetFile, nodedata::Matrix{<:Real}, name::AbstractString; component_names=nothing)
+    return WriteVTK.vtk_point_data(vtk, nodedata, name; component_names=component_names)
 end
 
 function component_names(::Type{S}) where S
@@ -113,46 +235,150 @@ function component_names(::Type{S}) where S
     return names
 end
 
-function vtk_nodeset(vtk::WriteVTK.DatasetFile, grid::Grid{dim}, nodeset::String) where {dim}
+"""
+    write_solution(vtk::VTKFile, dh::AbstractDofHandler, u::Vector, suffix="")
+
+Save the values at the nodes in the degree of freedom vector `u` to `vtk`.
+Each field in `dh` will be saved separately, and `suffix` can be used to append 
+to the fieldname.
+
+`u` can also contain tensorial values, but each entry in `u` must correspond to a 
+degree of freedom in `dh`, see [`write_node_data`](@ref write_node_data) for details. 
+Use `write_node_data` directly when exporting values that are already 
+sorted by the nodes in the grid. 
+"""
+function write_solution(vtk::VTKFile, dh::AbstractDofHandler, u::Vector, suffix="")
+    fieldnames = Ferrite.getfieldnames(dh)  # all primary fields
+    for name in fieldnames
+        data = _evaluate_at_grid_nodes(dh, u, name, #=vtk=# Val(true))
+        _vtk_write_node_data(vtk.vtk, data, string(name, suffix))
+    end
+    return vtk
+end
+
+"""
+    write_projection(vtk::VTKFile, proj::L2Projector, vals::Vector, name::AbstractString)
+
+Project `vals` to the grid nodes with `proj` and save to `vtk`.
+"""
+function write_projection(vtk::VTKFile, proj::L2Projector, vals, name)
+    data = _evaluate_at_grid_nodes(proj, vals, #=vtk=# Val(true))::Matrix
+    @assert size(data, 2) == getnnodes(get_grid(proj.dh))
+    _vtk_write_node_data(vtk.vtk, data, name; component_names=component_names(eltype(vals)))
+    return vtk
+end
+
+"""
+    write_cell_data(vtk::VTKFile, celldata::AbstractVector, name::String)
+
+Write the `celldata` that is ordered by the cells in the grid to the vtk file.
+"""
+function write_cell_data(vtk::VTKFile, celldata, name)
+    WriteVTK.vtk_cell_data(vtk.vtk, celldata, name)
+end
+
+"""
+    write_node_data(vtk::VTKFile, nodedata::Vector{Real}, name)
+    write_node_data(vtk::VTKFile, nodedata::Vector{<:AbstractTensor}, name)
+    
+Write the `nodedata` that is ordered by the nodes in the grid to `vtk`.
+
+When `nodedata` contains `Tensors.Vec`s, each component is exported. 
+Two-dimensional vectors are padded with zeros.
+
+When `nodedata` contains second order tensors, the index order, 
+`[11, 22, 33, 23, 13, 12, 32, 31, 21]`, follows the default Voigt order in Tensors.jl.
+"""
+function write_node_data(vtk::VTKFile, nodedata, name)
+    _vtk_write_node_data(vtk.vtk, nodedata, name)
+end
+
+
+"""
+    write_nodeset(vtk::VTKFile, grid::AbstractGrid, nodeset::String)
+
+Write nodal values of 1 for nodes in `nodeset`, and 0 otherwise
+"""
+function write_nodeset(vtk, grid::AbstractGrid, nodeset::String)
     z = zeros(getnnodes(grid))
     z[collect(getnodeset(grid, nodeset))] .= 1.0
-    vtk_point_data(vtk, z, nodeset)
+    write_node_data(vtk, z, nodeset)
+    return vtk
 end
 
 """
-    vtk_cellset(vtk, grid::Grid)
+    write_cellset(vtk, grid::AbstractGrid)
+    write_cellset(vtk, grid::AbstractGrid, cellset::String)
+    write_cellset(vtk, grid::AbstractGrid, cellsets::Union{AbstractVector{String},AbstractSet{String})
 
-Export all cell sets in the grid. Each cell set is exported with
-`vtk_cell_data` with value 1 if the cell is in the set, and 0 otherwise.
+Write all cell sets in the grid with name according to their keys and 
+celldata 1 if the cell is in the set, and 0 otherwise. It is also possible to 
+only export a single `cellset`, or multiple `cellsets`. 
 """
-function vtk_cellset(vtk::WriteVTK.DatasetFile, grid::AbstractGrid, cellsets=keys(grid.cellsets))
+function write_cellset(vtk, grid::AbstractGrid, cellsets=keys(getcellsets(getgrid(vtk))))
     z = zeros(getncells(grid))
     for cellset in cellsets
-        z .= 0.0
+        fill!(z, 0)
         z[collect(getcellset(grid, cellset))] .= 1.0
-        vtk_cell_data(vtk, z, cellset)
+        write_cell_data(vtk, z, cellset)
     end
     return vtk
 end
+write_cellset(vtk, grid::AbstractGrid, cellset::String) = write_cellset(vtk, grid, [cellset])
 
 """
-    vtk_cellset(vtk, grid::Grid, cellset::String)
+    write_constraints(vtk::VTKFile, ch::ConstraintHandler)
 
-Export the cell set specified by `cellset` as cell data with value 1 if
-the cell is in the set and 0 otherwise.
+Saves the dirichlet boundary conditions to a vtkfile.
+Values will have a 1 where bcs are active and 0 otherwise
 """
-vtk_cellset(vtk::WriteVTK.DatasetFile, grid::AbstractGrid, cellset::String) =
-    vtk_cellset(vtk, grid, [cellset])
+function write_constraints(vtk, ch::ConstraintHandler)    
+    unique_fields = []
+    for dbc in ch.dbcs
+        push!(unique_fields, dbc.field_name)
+    end
+    unique!(unique_fields)
 
+    for field in unique_fields
+        nd = getfielddim(ch.dh, field)
+        data = zeros(Float64, nd, getnnodes(get_grid(ch.dh)))
+        for dbc in ch.dbcs
+            dbc.field_name != field && continue
+            if eltype(dbc.facets) <: BoundaryIndex
+                functype = boundaryfunction(eltype(dbc.facets))
+                for (cellidx, facetidx) in dbc.facets
+                    for facetnode in functype(getcells(get_grid(ch.dh), cellidx))[facetidx]
+                        for component in dbc.components
+                            data[component, facetnode] = 1
+                        end
+                    end
+                end
+            else
+                for nodeidx in dbc.facets
+                    for component in dbc.components
+                        data[component, nodeidx] = 1
+                    end
+                end
+            end
+        end
+        write_node_data(vtk, data, string(field, "_bc"))
+    end
+    return vtk
+end
 
-function WriteVTK.vtk_point_data(vtkfile, dh::AbstractDofHandler, u::Vector, suffix="")
+"""
+    write_cell_colors(vtk::VTKFile, grid::AbstractGrid, cell_colors, name="coloring")
 
-    fieldnames = Ferrite.getfieldnames(dh)  # all primary fields
+Write cell colors (see [`create_coloring`](@ref)) to a VTK file for visualization.
 
-    for name in fieldnames
-        data = _evaluate_at_grid_nodes(dh, u, name, #=vtk=# Val(true))
-        vtk_point_data(vtkfile, data, string(name, suffix))
+In case of coloring a subset, the cells which are not part of the subset are represented as color 0.
+"""
+function write_cell_colors(vtk, grid::AbstractGrid, cell_colors::AbstractVector{<:AbstractVector{<:Integer}}, name="coloring")
+    color_vector = zeros(Int, getncells(grid))
+    for (i, cells_color) in enumerate(cell_colors)
+        for cell in cells_color
+            color_vector[cell] = i
+        end
     end
-
-    return vtkfile
+    write_cell_data(vtk, color_vector, name)
 end
diff --git a/src/Ferrite.jl b/src/Ferrite.jl
index 4ba6282a87..186e428f35 100644
--- a/src/Ferrite.jl
+++ b/src/Ferrite.jl
@@ -2,7 +2,6 @@ module Ferrite
 
 using Reexport: @reexport
 @reexport using Tensors
-@reexport using WriteVTK
 
 using Base:
     @propagate_inbounds
@@ -17,11 +16,12 @@ using SparseArrays:
     SparseArrays, SparseMatrixCSC, nonzeros, nzrange, rowvals, sparse, spzeros
 using StaticArrays:
     StaticArrays, MMatrix, SMatrix, SVector
+using WriteVTK:
+    WriteVTK, VTKCellTypes
 using Tensors:
     Tensors, AbstractTensor, SecondOrderTensor, SymmetricTensor, Tensor, Vec, gradient,
     rotation_tensor, symmetric, tovoigt!
-using WriteVTK:
-    WriteVTK, MeshCell, VTKCellTypes, vtk_cell_data, vtk_grid, vtk_point_data, vtk_save
+
 
 include("exports.jl")
 
diff --git a/src/Grid/coloring.jl b/src/Grid/coloring.jl
index db30efd5d0..9b4e18ffd8 100644
--- a/src/Grid/coloring.jl
+++ b/src/Grid/coloring.jl
@@ -205,20 +205,3 @@ function create_coloring(g::AbstractGrid, cellset=1:getncells(g); alg::ColoringA
         error("impossible")
     end
 end
-
-"""
-    vtk_cell_data_colors(vtkfile, cell_colors, name="coloring")
-
-Write cell colors (see [`create_coloring`](@ref)) to a VTK file for visualization.
-
-In case of coloring a subset, the cells which are not part of the subset are represented as color 0.
-"""
-function vtk_cell_data_colors(vtkfile, cell_colors::AbstractVector{<:AbstractVector{<:Integer}}, name="coloring")
-    color_vector = zeros(Int, vtkfile.Ncls)
-    for (i, cells_color) in enumerate(cell_colors)
-        for cell in cells_color
-            color_vector[cell] = i
-        end
-    end
-    vtk_cell_data(vtkfile, color_vector, name)
-end
diff --git a/src/L2_projection.jl b/src/L2_projection.jl
index d9bc37c939..e647c3cbf2 100644
--- a/src/L2_projection.jl
+++ b/src/L2_projection.jl
@@ -223,13 +223,6 @@ function _project(vars, proj::L2Projector, fe_values::AbstractValues, M::Integer
     return T[make_T(x) for x in eachrow(projected_vals)]
 end
 
-function WriteVTK.vtk_point_data(vtk::WriteVTK.DatasetFile, proj::L2Projector, vals::Vector{T}, name::AbstractString) where T
-    data = _evaluate_at_grid_nodes(proj, vals, #=vtk=# Val(true))::Matrix
-    @assert size(data, 2) == getnnodes(get_grid(proj.dh))
-    vtk_point_data(vtk, data, name; component_names=component_names(T))
-    return vtk
-end
-
 evaluate_at_grid_nodes(proj::L2Projector, vals::AbstractVector) =
     _evaluate_at_grid_nodes(proj, vals, Val(false))
 
diff --git a/src/deprecations.jl b/src/deprecations.jl
index 8f2561b846..9c9d9c7c17 100644
--- a/src/deprecations.jl
+++ b/src/deprecations.jl
@@ -59,6 +59,17 @@ Base.@deprecate_binding Line3D Line
 Base.@deprecate_binding Quadrilateral3D Quadrilateral
 export Line2D, Line3D, Quadrilateral3D
 
+using WriteVTK: vtk_grid 
+export vtk_grid # To give better error
+
+function WriteVTK.vtk_grid(::String, ::Union{AbstractGrid,AbstractDofHandler}; kwargs...)
+    error(join(("The vtk interface has been updated in Ferrite v1.0.",
+                "See https://github.com/Ferrite-FEM/Ferrite.jl/pull/679.",
+                "Use VTKFile to open a vtk file, and the functions",
+                "write_solution, write_cell_data, and write_projection to save data."), 
+            "\n"))
+end
+
 # Deprecation of auto-vectorized methods
 function add!(dh::DofHandler, name::Symbol, dim::Int)
     celltype = getcelltype(dh.grid)
diff --git a/src/exports.jl b/src/exports.jl
index c633f52880..d2daeb8eb6 100644
--- a/src/exports.jl
+++ b/src/exports.jl
@@ -106,7 +106,6 @@ export
 # Grid coloring
     create_coloring,
     ColoringAlgorithm,
-    vtk_cell_data_colors,
 
 # Dofs
     DofHandler,
@@ -159,14 +158,15 @@ export
     assemble!,
     finish_assemble,
 
-# VTK export
-    vtk_grid,
-    vtk_point_data,
-    vtk_cell_data,
-    vtk_nodeset,
-    vtk_cellset,
-    vtk_save,
-
+# exporting data
+    VTKFile,
+    write_solution,
+    write_cell_data,
+    write_projection,
+    write_node_data,
+    VTKFileCollection,
+    addstep!,
+    
 # L2 Projection
     project,
     L2Projector,
diff --git a/test/runtests.jl b/test/runtests.jl
index f864f31b9a..3d6498c011 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -42,6 +42,7 @@ include("test_assemble.jl")
 include("test_dofs.jl")
 include("test_constraints.jl")
 include("test_grid_dofhandler_vtk.jl")
+include("test_vtk_export.jl")
 include("test_abstractgrid.jl")
 include("test_grid_addboundaryset.jl")
 include("test_mixeddofhandler.jl")
diff --git a/test/test_constraints.jl b/test/test_constraints.jl
index 5bbb4468e7..94161ea981 100644
--- a/test/test_constraints.jl
+++ b/test/test_constraints.jl
@@ -1553,9 +1553,9 @@ end # testset
             @test norm(u_dbc) ≈ 3.8249286998373586
             @test norm(u_p) ≈ 3.7828270430540893
         end
-        # vtk_grid("local_application_azero_$(azero)", grid) do vtk
-        #     vtk_point_data(vtk, dh, u_dbc, "_dbc")
-        #     vtk_point_data(vtk, dh, u_p, "_p")
+        # VTKFile("local_application_azero_$(azero)", grid) do vtk
+        #     write_solution(vtk, dh, u_dbc, "_dbc")
+        #     write_solution(vtk, dh, u_p, "_p")
         # end
         @test K_dbc_standard \ f_dbc_standard ≈ K_dbc_ch \ f_dbc_ch ≈ K_dbc_local \ f_dbc_local ≈
               K_ac_standard \ f_ac_standard ≈ K_ac_ch \ f_ac_ch ≈ K_ac_local \ f_ac_local
diff --git a/test/test_deprecations.jl b/test/test_deprecations.jl
index 9547bfea8c..3e7a731eae 100644
--- a/test/test_deprecations.jl
+++ b/test/test_deprecations.jl
@@ -96,4 +96,9 @@ addfaceset!(grid, "right_face_explicit", Set(Ferrite.FaceIndex(fi[1], fi[2]) for
 @test getfacetset(grid, "right_face_explicit") == getfacetset(grid, "right")
 end
 
+@testset "vtk_grid" begin
+    # Ensure no MethodError on pre v1. 
+    @test_throws ErrorException vtk_grid("old", generate_grid(Line, (1,)))
+end
+
 end # testset deprecations
diff --git a/test/test_grid_dofhandler_vtk.jl b/test/test_grid_dofhandler_vtk.jl
index 776d0916a7..6ac0cba198 100644
--- a/test/test_grid_dofhandler_vtk.jl
+++ b/test/test_grid_dofhandler_vtk.jl
@@ -38,10 +38,10 @@ end
         addnodeset!(grid, "middle-nodes", x -> norm(x) < radius)
 
         gridfilename = "grid-$(repr(celltype))"
-        vtk_grid(gridfilename, grid) do vtk
-            vtk_cellset(vtk, grid, "cell-1")
-            vtk_cellset(vtk, grid, "middle-cells")
-            vtk_nodeset(vtk, grid, "middle-nodes")
+        VTKFile(gridfilename, grid) do vtk
+            Ferrite.write_cellset(vtk, grid, "cell-1")
+            Ferrite.write_cellset(vtk, grid, "middle-cells")
+            Ferrite.write_nodeset(vtk, grid, "middle-nodes")
         end
 
         # test the sha of the file
@@ -77,9 +77,9 @@ end
         apply!(u, ch)
 
         dofhandlerfilename = "dofhandler-$(repr(celltype))"
-        vtk_grid(dofhandlerfilename, dofhandler) do vtk
-            vtk_point_data(vtk, ch)
-            vtk_point_data(vtk, dofhandler, u)
+        VTKFile(dofhandlerfilename, grid) do vtk
+            Ferrite.write_constraints(vtk, ch)
+            write_solution(vtk, dofhandler, u)
         end
 
         # test the sha of the file
@@ -118,10 +118,10 @@ close(csio)
     vector_data = [Vec{3}(ntuple(i->i, 3)) for j=1:8]
 
     filename_3d = "test_vtk_3d"
-    vtk_grid(filename_3d, grid) do vtk_file
-        vtk_point_data(vtk_file, sym_tensor_data, "symmetric tensor")
-        vtk_point_data(vtk_file, tensor_data, "tensor")
-        vtk_point_data(vtk_file, vector_data, "vector")
+    VTKFile(filename_3d, grid) do vtk
+        write_node_data(vtk, sym_tensor_data, "symmetric tensor")
+        write_node_data(vtk, tensor_data, "tensor")
+        write_node_data(vtk, vector_data, "vector")
     end
 
     # 2D grid
@@ -133,11 +133,11 @@ close(csio)
     vector_data = [Vec{2}(ntuple(i->i, 2)) for j=1:4]
 
     filename_2d = "test_vtk_2d"
-    vtk_grid(filename_2d, grid) do vtk_file
-        vtk_point_data(vtk_file, sym_tensor_data, "symmetric tensor")
-        vtk_point_data(vtk_file, tensor_data, "tensor")
-        vtk_point_data(vtk_file, tensor_data_1D, "tensor_1d")
-        vtk_point_data(vtk_file, vector_data, "vector")
+    VTKFile(filename_2d, grid) do vtk
+        write_node_data(vtk, sym_tensor_data, "symmetric tensor")
+        write_node_data(vtk, tensor_data, "tensor")
+        write_node_data(vtk, tensor_data_1D, "tensor_1d")
+        write_node_data(vtk, vector_data, "vector")
     end
 
     # test the shas of the files
@@ -708,4 +708,53 @@ end
         close!(dh2)
         @test dh1.cell_dofs == dh2.cell_dofs
     end
+
+    @testset "VTKFileCollection" begin
+        @testset "equivalence of addstep methods" begin 
+            grid = generate_grid(Triangle, (2,2))
+            celldata = rand(getncells(grid))
+            fname = "addstep"
+            pvd1 = VTKFileCollection(fname, grid)
+            pvd2 = VTKFileCollection(fname, grid)
+            timesteps = 0:0.5:0.5
+            for (n, t) in pairs(timesteps)
+                addstep!(pvd1, t) do io
+                    write_cell_data(io, celldata*n, "celldata")
+                end
+                vtk = VTKFile(string(fname, "2_", n), grid)
+                write_cell_data(vtk, celldata*n, "celldata")
+                addstep!(pvd2, vtk, t)
+                @test !(isopen(vtk.vtk))
+            end
+            close.((pvd1, pvd2))
+            @test pvd1.step == pvd2.step # Same nr of steps added 
+            for (n, t) in pairs(timesteps)
+                fname1 = string(fname, "_", n, ".vtu")
+                fname2 = string(fname, "2_", n, ".vtu")
+                sha_vtk1 = bytes2hex(open(SHA.sha1, fname1))
+                sha_vtk2 = bytes2hex(open(SHA.sha1, fname2))
+                @test sha_vtk1 == sha_vtk2
+                rm.((fname1, fname2))
+            end
+            rm(string(fname, ".pvd"))
+            # Solving https://github.com/Ferrite-FEM/Ferrite.jl/issues/397
+            # would allow checking the pvd files as well. 
+        end
+        @testset "kwargs forwarding" begin
+            grid = generate_grid(Quadrilateral, (10,10))
+            file_sizes = Int[]
+            fname = "test_collection_kwargs"
+            for compress in (true, false)
+                pvd = VTKFileCollection(fname, grid)
+                addstep!(pvd, 0.0; compress) do io
+                    nothing
+                end
+                close(pvd)
+                push!(file_sizes, stat(string(fname, "_1.vtu")).size)
+                rm(string(fname, "_1.vtu"))
+            end
+            rm(string(fname, ".pvd"))
+            @test file_sizes[1] < file_sizes[2] # Check that compress=true gives smaller file size
+        end
+    end
 end
diff --git a/test/test_l2_projection.jl b/test/test_l2_projection.jl
index 6a2a979198..7f046306cb 100644
--- a/test/test_l2_projection.jl
+++ b/test/test_l2_projection.jl
@@ -281,16 +281,17 @@ function test_export(;subset::Bool)
     end
 
     mktempdir() do tmp
-        fname = vtk_grid(joinpath(tmp, "projected"), grid) do vtk
-            vtk_point_data(vtk, p, p_scalar, "p_scalar")
-            vtk_point_data(vtk, p, p_vec, "p_vec")
-            vtk_point_data(vtk, p, p_tens, "p_tens")
-            vtk_point_data(vtk, p, p_stens, "p_stens")
+        fname = joinpath(tmp, "projected")
+        VTKFile(fname, grid) do vtk
+            write_projection(vtk, p, p_scalar, "p_scalar")
+            write_projection(vtk, p, p_vec, "p_vec")
+            write_projection(vtk, p, p_tens, "p_tens")
+            write_projection(vtk, p, p_stens, "p_stens")
         end
         # The following test may fail due to floating point inaccuracies
         # These could occur due to e.g. changes in system architecture.
         if Sys.islinux() && Sys.ARCH === :x86_64
-            @test bytes2hex(open(SHA.sha1, fname[1], "r")) == (
+            @test bytes2hex(open(SHA.sha1, fname*".vtu", "r")) == (
                 subset ? "b3fef3de9f38ca9ddd92f2f67a1606d07ca56d67" :
                          "bc2ec8f648f9b8bccccf172c1fc48bf03340329b"
             )
diff --git a/test/test_mixeddofhandler.jl b/test/test_mixeddofhandler.jl
index 8e2d548e52..c9a0b50603 100644
--- a/test/test_mixeddofhandler.jl
+++ b/test/test_mixeddofhandler.jl
@@ -357,11 +357,11 @@ function test_2_element_heat_eq()
     gridfilename = "mixed_grid"
     addcellset!(grid, "cell-1", [1,])
     addcellset!(grid, "cell-2", [2,])
-    vtk_grid(gridfilename, grid) do vtk
-        vtk_cellset(vtk, grid, "cell-1")
-        vtk_cellset(vtk, grid, "cell-2")
-        vtk_point_data(vtk, dh, u)
-        # vtk_point_data(vtk, ch)  #FIXME
+    VTKFile(gridfilename, grid) do vtk
+        Ferrite.write_cellset(vtk, grid, "cell-1")
+        Ferrite.write_cellset(vtk, grid, "cell-2")
+        write_solution(vtk, dh, u)
+        # Ferrite.write_constraints(vtk, ch)  #FIXME
     end
     sha = bytes2hex(open(SHA.sha1, gridfilename*".vtu"))
     @test sha in ("e96732c000b0b385db7444f002461468b60b3b2c", "7b26edc27b5e59a2f60907374cd5a5790cc37a6a")
@@ -635,8 +635,8 @@ function test_vtk_export()
     close!(dh)
     u = collect(1:ndofs(dh))
     filename = "mixed_2d_grid"
-    vtk_grid(filename, dh) do vtk
-        vtk_point_data(vtk, dh, u)
+    VTKFile(filename, dh) do vtk
+        write_solution(vtk, dh, u)
     end
     sha = bytes2hex(open(SHA.sha1, filename*".vtu"))
     @test sha == "339ab8a8a613c2f38af684cccd695ae816671607"
diff --git a/test/test_vtk_export.jl b/test/test_vtk_export.jl
new file mode 100644
index 0000000000..97c88fe20f
--- /dev/null
+++ b/test/test_vtk_export.jl
@@ -0,0 +1,44 @@
+@testset "VTKFile" begin #TODO: Move all vtk tests here 
+    @testset "show(::VTKFile)" begin
+        mktempdir() do tmp    
+            grid = generate_grid(Quadrilateral, (2,2))
+            vtk = VTKFile(joinpath(tmp, "showfile"), grid)
+            showstring_open = sprint(show, MIME"text/plain"(), vtk)
+            @test startswith(showstring_open, "VTKFile for the open file")
+            @test contains(showstring_open, "showfile.vtu")
+            close(vtk)
+            showstring_closed = sprint(show, MIME"text/plain"(), vtk)
+            @test startswith(showstring_closed, "VTKFile for the closed file")
+            @test contains(showstring_closed, "showfile.vtu")
+        end
+    end
+    @testset "cellcolors" begin
+        mktempdir() do tmp
+            grid = generate_grid(Quadrilateral, (4, 4))
+            colors = create_coloring(grid)
+            fname = joinpath(tmp, "colors")
+            VTKFile(fname, grid) do vtk
+                Ferrite.write_cell_colors(vtk, grid, colors)
+            end
+            @test bytes2hex(open(SHA.sha1, fname*".vtu")) == "b804d0b064121b672d8e35bcff8446eda361cac3"
+        end
+    end
+    @testset "constraints" begin
+        mktempdir() do tmp
+            grid = generate_grid(Tetrahedron, (4, 4, 4))
+            dh = DofHandler(grid)
+            add!(dh, :u, Lagrange{RefTetrahedron, 1}())
+            close!(dh)
+            ch = ConstraintHandler(dh)
+            add!(ch, Dirichlet(:u, getfacetset(grid, "left"), x -> 0.0))
+            addnodeset!(grid, "nodeset", x -> x[1] ≈ 1.0)
+            add!(ch, Dirichlet(:u, getnodeset(grid, "nodeset"), x -> 0.0))
+            close!(ch)
+            fname = joinpath(tmp, "constraints")
+            VTKFile(fname, grid) do vtk
+                Ferrite.write_constraints(vtk, ch)
+            end
+            @test bytes2hex(open(SHA.sha1, fname*".vtu")) == "31b506bd9729b11992f8bcb79a2191eb65d223bf"
+        end
+    end
+end