Abstraction layer around VTK export (#692)

This patch introduces a Ferrite-specific abstraction layer on top of
WriteVTK.jl instead of re-exporting and extending methods.

VTK files are created with the `VTKFile` constructor (instead of
`WriteVTK.vtk_grid`) and closed with `Base.close`. To add data to the
file the following functions, which replaces methods of
`WriteVTK.vtk_point_data` and `WriteVTK.vtk_cell_data`, should be used:

 - `write_solution(::VTKFile, ::DofHandler, ::Vector)`: write a solution
   vector corresponding to DoFs in the DofHandler.
 - `write_projection(::VTKFile, ::L2Projector, ::Vector)`: write a
   projected vector corresponding the DoFs in the DofHandler of the
   projector.
 - `write_cell_data(::VTKFile, ::Vector)`: write data per cell.
 - `write_node_data(::VTKFile, ::Vector)`: write data per node.
 - `Ferrite.write_cellset`, `Ferrite.write_nodeset`: visualize cellsets
   and nodesets.
 - `Ferrite.write_constraints`: visualize constrained degrees of
   freedom.
 - `Ferrite.write_cell_colors`: visualize the result of grid coloring.

To collect multiple time steps the `VTKFileCollection` struct is used
(instead of `WriteVTK.paraview_collection`). VTK files are added to the
collection with `addstep!`.

Fixes #278, fixes #686. Closes #691.

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
-

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

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,

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;

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)

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

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
 

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;
 

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

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

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
 

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

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

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

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

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

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

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 #

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
 ```