diff --git a/deps/build.jl b/deps/build.jl
index 3032ec4..829dcf0 100644
--- a/deps/build.jl
+++ b/deps/build.jl
@@ -36,10 +36,11 @@ files = [
 
 Sys.rm(notebooks_dir;recursive=true,force=true)
 for (i,(title,filename)) in enumerate(files)
-  notebook_prefix = string("t",@sprintf "%03d_" i)
-  notebook = string(notebook_prefix,splitext(filename)[1])
+  notebook_prefix = string("t",@sprintf "%03d" i)
+  notebook = string(notebook_prefix,"_",splitext(filename)[1])
   notebook_title = string("# # Tutorial ", i, ": ", title)
   function preprocess_notebook(content)
+    content = replace(content, "output_path" => notebook_prefix)
     return string(notebook_title, "\n\n", content)
   end
   Literate.notebook(joinpath(repo_src,filename), notebooks_dir; name=notebook, preprocess=preprocess_notebook, documenter=false, execute=false)
diff --git a/src/TopOptEMFocus.jl b/src/TopOptEMFocus.jl
index 2686996..c6b7725 100644
--- a/src/TopOptEMFocus.jl
+++ b/src/TopOptEMFocus.jl
@@ -421,12 +421,14 @@ function gf_p(p0::Vector, grad::Vector; r, β, η, phys_params, fem_params)
         grad[:] = dgdp
     end
     gvalue = gf_p(p0::Vector; r, β, η, phys_params, fem_params)
-    open("gvalue.txt", "a") do io
+    open("output_path/gvalue.txt", "a") do io
         write(io, "$gvalue \n")
     end
     gvalue
 end
 
+mkpath("output_path")
+
 # Using the following codes, we can check if we can get the derivatives correctly from the adjoint method by comparing it with the finite difference results.
 #
 
@@ -494,7 +496,7 @@ ax.aspect = AxisAspect(1)
 ax.title = "Design Shape"
 rplot = 110 # Region for plot
 limits!(ax, -rplot, rplot, (h1)/2-rplot, (h1)/2+rplot)
-save("shape.png", fig)
+save("output_path/shape.png", fig)
 
 
 # ![](../assets/TopOptEMFocus/shape.png)
@@ -512,7 +514,7 @@ Colorbar(fig[1,2], plt)
 ax.title = "|E|"
 ax.aspect = AxisAspect(1)
 limits!(ax, -rplot, rplot, (h1)/2-rplot, (h1)/2+rplot)
-save("Field.png", fig)
+save("output_path/Field.png", fig)
 
 # ![](../assets/TopOptEMFocus/Field.png)
 #
diff --git a/src/advection_diffusion.jl b/src/advection_diffusion.jl
index 6900ef3..bdf2631 100644
--- a/src/advection_diffusion.jl
+++ b/src/advection_diffusion.jl
@@ -162,9 +162,11 @@ uh_non_zero = Gridap.Algebra.solve(op_non_zero)
 
 # Ouput the result as a `.vtk` file.
 
-writevtk(Ω,"results_zero",cellfields=["uh_zero"=>uh_zero])
+mkpath("output_path")
 
-writevtk(Ω,"results_non_zero",cellfields=["uh_non_zero"=>uh_non_zero])
+writevtk(Ω,"output_path/results_zero",cellfields=["uh_zero"=>uh_zero])
+
+writevtk(Ω,"output_path/results_non_zero",cellfields=["uh_non_zero"=>uh_non_zero])
 
 # ## Visualization
 
diff --git a/src/darcy.jl b/src/darcy.jl
index a99e7e7..69421ab 100644
--- a/src/darcy.jl
+++ b/src/darcy.jl
@@ -136,7 +136,8 @@ uh, ph = xh
 
 # Since this is a multi-field example, the `solve` function returns a multi-field solution `xh`, which can be unpacked in order to finally recover each field of the problem. The resulting single-field objects can be visualized as in previous tutorials (see next figure).
 
-writevtk(trian,"darcyresults",cellfields=["uh"=>uh,"ph"=>ph])
+mkpath("output_path")
+writevtk(trian,"output_path/darcyresults",cellfields=["uh"=>uh,"ph"=>ph])
 
 # ![](../assets/darcy/darcy_results.png)
 #
diff --git a/src/dg_discretization.jl b/src/dg_discretization.jl
index 7b444a7..7fb31dc 100644
--- a/src/dg_discretization.jl
+++ b/src/dg_discretization.jl
@@ -162,7 +162,8 @@ model = CartesianDiscreteModel(domain,partition)
 # to be generated in each direction (here $4\times4\times4$ cells). You can
 # write the model in vtk format to visualize it (see next figure).
 
- writevtk(model,"model")
+mkpath("output_path")
+writevtk(model,"output_path/model")
 
 # ![](../assets/dg_discretization/model.png)
 # 
@@ -240,7 +241,7 @@ U = TrialFESpace(V)
 # written into a vtk file for its visualization (see next figure, where the
 # interior facets $\mathcal{F}_\Lambda$ are clearly observed).
 
-writevtk(Λ,"strian")
+writevtk(Λ,"output_path/strian")
 
 # ![](../assets/dg_discretization/skeleton_trian.png)
 # 
@@ -329,7 +330,7 @@ uh = solve(op)
 # faces. We compute and visualize the jump of these values as follows (see next
 # figure):
 
-writevtk(Λ,"jumps",cellfields=["jump_u"=>jump(uh)])
+writevtk(Λ,"output_path/jumps",cellfields=["jump_u"=>jump(uh)])
 
 # Note that the jump of the numerical solution is very small, close to the
 # machine precision (as expected in this example with manufactured solution).
diff --git a/src/elasticity.jl b/src/elasticity.jl
index 0069ecc..ed22c9f 100644
--- a/src/elasticity.jl
+++ b/src/elasticity.jl
@@ -49,7 +49,8 @@ model = DiscreteModelFromFile("../models/solid.json")
 
 # In order to inspect it, write the model to vtk
 
-writevtk(model,"model")
+mkpath("output_path")
+writevtk(model,"output_path/model")
 
 # and open the resulting files with Paraview. The boundaries $\Gamma_{\rm B}$ and $\Gamma_{\rm G}$ are identified  with the names `"surface_1"` and `"surface_2"` respectively.  For instance, if you visualize the faces of the model and color them by the field `"surface_2"` (see next figure), you will see that only the faces on $\Gamma_{\rm G}$ have a value different from zero.
 #
@@ -118,7 +119,7 @@ uh = solve(op)
 #
 # Finally, we write the results to a file. Note that we also include the strain and stress tensors into the results file.
 
-writevtk(Ω,"results",cellfields=["uh"=>uh,"epsi"=>ε(uh),"sigma"=>σ∘ε(uh)])
+writevtk(Ω,"output_path/results",cellfields=["uh"=>uh,"epsi"=>ε(uh),"sigma"=>σ∘ε(uh)])
 
 # It can be clearly observed (see next figure) that the surface  $\Gamma_{\rm B}$ is pulled in $x_1$-direction and that the solid deforms accordingly.
 #
@@ -186,7 +187,7 @@ uh = solve(op)
 
 # Once the solution is computed, we can store the results in a file for visualization. Note that, we are including the stress tensor in the file (computed with the bi-material law).
 
-writevtk(Ω,"results_bimat",cellfields=
+writevtk(Ω,"output_path/results_bimat",cellfields=
   ["uh"=>uh,"epsi"=>ε(uh),"sigma"=>σ_bimat∘(ε(uh),tags)])
 
 # #### Constant multi-material law
@@ -200,6 +201,6 @@ tags_field = CellField(tags, Ω)
 # `tags_field` is a field which value at $x$ is the tag of the cell containing $x$. `σ_bimat_cst` is used like a constant in (bi)linear form definition and solution export:
 
 a(u,v) = ∫( σ_bimat_cst * ∇(u)⋅∇(v))*dΩ
-writevtk(Ω,"const_law",cellfields= ["sigma"=>σ_bimat_cst])
+writevtk(Ω,"output_path/const_law",cellfields= ["sigma"=>σ_bimat_cst])
 
 #  Tutorial done!
diff --git a/src/emscatter.jl b/src/emscatter.jl
index 570ee3b..2310f1a 100644
--- a/src/emscatter.jl
+++ b/src/emscatter.jl
@@ -291,7 +291,8 @@ uh_t = CellField(x->H_t(x,xc,r,ϵ₁,λ),Ω)
 # ![](../assets/emscatter/Results.png)
 
 # ### Save to file and view
-writevtk(Ω,"demo",cellfields=["Real"=>real(uh),
+mkpath("output_path")
+writevtk(Ω,"output_path/demo",cellfields=["Real"=>real(uh),
         "Imag"=>imag(uh),
         "Norm"=>abs2(uh),
         "Real_t"=>real(uh_t),
diff --git a/src/fsi_tutorial.jl b/src/fsi_tutorial.jl
index b1ba2a4..076ded1 100644
--- a/src/fsi_tutorial.jl
+++ b/src/fsi_tutorial.jl
@@ -76,7 +76,8 @@ using Gridap
 model = DiscreteModelFromFile("../models/elasticFlag.json")
 
 # We can inspect the loaded geometry and associated parts by printing to a `vtk` file:
-writevtk(model,"model")
+mkpath("output_path")
+writevtk(model,"output_path/model")
 
 # This will produce an output in which we can identify the different parts of the domain, with the associated labels and tags.
 #
@@ -357,12 +358,12 @@ uhs, uhf, ph = solve(op)
 # ```
 # ### Visualization
 # The solution fields $[\mathbf{U}^h_{\rm S},\mathbf{U}^h_{\rm F},\mathbf{P}^h_{\rm F}]^T$ are defined over all the domain, extended with zeros on the inactive part. Calling the function `writevtk` passing the global triangulation, we will output the global fields.
-writevtk(Ω,"trian", cellfields=["uhs" => uhs, "uhf" => uhf, "ph" => ph])
+writevtk(Ω,"output_path/trian", cellfields=["uhs" => uhs, "uhf" => uhf, "ph" => ph])
 # ![](../assets/fsi/Global_solution.png)
 
 # However, we can also restrict the fields to the active part by calling the function `restrict` with the field along with the respective active triangulation.
-writevtk(Ω_s,"trian_solid",cellfields=["uhs"=>uhs])
-writevtk(Ω_f,"trian_fluid",cellfields=["ph"=>ph,"uhf"=>uhf])
+writevtk(Ω_s,"output_path/trian_solid",cellfields=["uhs"=>uhs])
+writevtk(Ω_f,"output_path/trian_fluid",cellfields=["ph"=>ph,"uhf"=>uhf])
 # ![](../assets/fsi/Local_solution.png)
 
 # ```@raw HTML
diff --git a/src/geometry_dev.jl b/src/geometry_dev.jl
index 68f7c7b..09e7b1c 100644
--- a/src/geometry_dev.jl
+++ b/src/geometry_dev.jl
@@ -225,7 +225,8 @@ reffes, cell_types = compress_cell_data(cell_reffes)
 new_grid = UnstructuredGrid(new_node_coordinates,cell_to_nodes,reffes,cell_types)
 
 # Save for visualization:
-writevtk(new_grid,"half_cylinder_linear")
+mkpath("output_path")
+writevtk(new_grid,"output_path/half_cylinder_linear")
 
 #
 # If we visualize the result, we'll notice that despite applying a curved mapping,
@@ -267,7 +268,7 @@ new_node_coordinates = map(F,new_node_coordinates)
 
 # Create the high-order grid:
 new_grid = UnstructuredGrid(new_node_coordinates,new_cell_to_nodes,new_reffes,cell_types)
-writevtk(new_grid,"half_cylinder_quadratic")
+writevtk(new_grid,"output_path/half_cylinder_quadratic")
 
 # The resulting mesh now accurately represents the curved geometry of the half-cylinder,
 # with quadratic elements properly capturing the curvature (despite paraview still showing 
@@ -427,7 +428,7 @@ node_to_entity = get_face_entity(labels,0) # For each node, its associated entit
 
 # It is usually more convenient to visualise it in Paraview by exporting to vtk: 
 
-writevtk(model,"labels_basic",labels=labels)
+writevtk(model,"output_path/labels_basic",labels=labels)
 
 # Another useful way to create a `FaceLabeling` is by providing a coloring for the mesh cells, 
 # where each color corresponds to a different tag.
@@ -436,21 +437,21 @@ writevtk(model,"labels_basic",labels=labels)
 cell_to_tag = [1,1,1,2,2,3,2,2,3]
 tag_to_name = ["A","B","C"]
 labels_cw = Geometry.face_labeling_from_cell_tags(topo,cell_to_tag,tag_to_name)
-writevtk(model,"labels_cellwise",labels=labels_cw)
+writevtk(model,"output_path/labels_cellwise",labels=labels_cw)
 
 # We can also create a `FaceLabeling` from a vertex filter. The resulting `FaceLabeling` will have
 # only one tag, gathering the d-faces whose vertices ALL fullfill `filter(x) == true`.
 
 vfilter(x) = abs(x[1]- 1.0) < 1.e-5
 labels_vf = Geometry.face_labeling_from_vertex_filter(topo, "top", vfilter)
-writevtk(model,"labels_filter",labels=labels_vf)
+writevtk(model,"output_path/labels_filter",labels=labels_vf)
 
 # `FaceLabeling` objects can also be merged together. The resulting `FaceLabeling` will have 
 # the union of the tags and entities of the original ones.
 # Note that this modifies the first `FaceLabeling` in place.
 
 labels = merge!(labels, labels_cw, labels_vf)
-writevtk(model,"labels_merged",labels=labels)
+writevtk(model,"output_path/labels_merged",labels=labels)
 
 # ### Creating new tags from existing ones
 #
@@ -474,7 +475,7 @@ Geometry.add_tag_from_tags_complementary!(labels,"!A",["A"])
 # and creates a new tag that contains all the d-faces that are in the first list but not in the second.
 Geometry.add_tag_from_tags_setdiff!(labels,"A-B",["A"],["B"]) # set difference
 
-writevtk(model,"labels_setops",labels=labels)
+writevtk(model,"output_path/labels_setops",labels=labels)
 
 # ### FaceLabeling queries
 #
diff --git a/src/hyperelasticity.jl b/src/hyperelasticity.jl
index f1edc31..bf83480 100644
--- a/src/hyperelasticity.jl
+++ b/src/hyperelasticity.jl
@@ -96,7 +96,8 @@ function run(x0,disp_x,step,nsteps,cache)
 
   uh, cache = solve!(uh,solver,op,cache)
 
-  writevtk(Ω,"results_$(lpad(step,3,'0'))",cellfields=["uh"=>uh,"sigma"=>σ∘∇(uh)])
+  mkpath("output_path")
+  writevtk(Ω,"output_path/results_$(lpad(step,3,'0'))",cellfields=["uh"=>uh,"sigma"=>σ∘∇(uh)])
 
   return get_free_dof_values(uh), cache
 
diff --git a/src/inc_navier_stokes.jl b/src/inc_navier_stokes.jl
index 11117b1..5ce7f9f 100644
--- a/src/inc_navier_stokes.jl
+++ b/src/inc_navier_stokes.jl
@@ -151,7 +151,8 @@ uh, ph = solve(solver,op)
 
 # Finally, we write the results for visualization (see next figure).
 
-writevtk(Ωₕ,"ins-results",cellfields=["uh"=>uh,"ph"=>ph])
+mkpath("output_path")
+writevtk(Ωₕ,"output_path/ins-results",cellfields=["uh"=>uh,"ph"=>ph])
 
 # ![](../assets/inc_navier_stokes/ins_solution.png)
 #
diff --git a/src/interpolation_fe.jl b/src/interpolation_fe.jl
index 36c7267..b1b0d7a 100644
--- a/src/interpolation_fe.jl
+++ b/src/interpolation_fe.jl
@@ -131,8 +131,9 @@ g̃ₕ.cell_dof_values
 
 # We can visualize the results using Paraview
 
-writevtk(get_triangulation(fₕ), "source", cellfields=["fₕ"=>fₕ])
-writevtk(get_triangulation(gₕ), "target", cellfields=["gₕ"=>gₕ])
+mkpath("output_path")
+writevtk(get_triangulation(fₕ), "output_path/source", cellfields=["fₕ"=>fₕ])
+writevtk(get_triangulation(gₕ), "output_path/target", cellfields=["gₕ"=>gₕ])
 
 # which produces the following output
 
diff --git a/src/isotropic_damage.jl b/src/isotropic_damage.jl
index cd90e3c..98bfd48 100644
--- a/src/isotropic_damage.jl
+++ b/src/isotropic_damage.jl
@@ -120,6 +120,8 @@ function main(;n,nsteps)
   uh = zero(V)
   cache = nothing
 
+  mkpath("output_path")
+  
   for (istep,factor) in enumerate(factors)
 
     println("\n+++ Solving for load factor $factor in step $istep of $nsteps +++\n")
@@ -129,7 +131,7 @@ function main(;n,nsteps)
     rh = project(r,model,dΩ,order)
 
     writevtk(
-      Ω,"results_$(lpad(istep,3,'0'))",
+      Ω,"output_path/results_$(lpad(istep,3,'0'))",
       cellfields=["uh"=>uh,"epsi"=>ε(uh),"damage"=>dh,
                   "threshold"=>rh,"sigma_elast"=>σe∘ε(uh)])
 
diff --git a/src/lagrange_multipliers.jl b/src/lagrange_multipliers.jl
index a9af103..841506c 100644
--- a/src/lagrange_multipliers.jl
+++ b/src/lagrange_multipliers.jl
@@ -124,4 +124,5 @@ l2_error = sqrt(sum(∫(eh⋅eh)*dΩ))
 #
 # We can visualize the solution and error by writing them to a VTK file:
 
-writevtk(Ω, "results", cellfields=["uh"=>uh, "error"=>eh])
+mkpath("output_path")
+writevtk(Ω, "output_path/results", cellfields=["uh"=>uh, "error"=>eh])
diff --git a/src/p_laplacian.jl b/src/p_laplacian.jl
index d8abff8..1b968be 100644
--- a/src/p_laplacian.jl
+++ b/src/p_laplacian.jl
@@ -48,7 +48,8 @@ model = DiscreteModelFromFile("../models/model.json")
 
 # As stated before, we want to impose Dirichlet boundary conditions on $\Gamma_0$ and $\Gamma_g$,  but none of these boundaries is identified in the model. E.g., you can easily see by writing the model in vtk format
 
-writevtk(model,"model")
+mkpath("output_path")
+writevtk(model,"output_path/model")
 
 # and by opening the file `"model_0"` in Paraview that the boundary identified as `"sides"` only includes the vertices in the interior of $\Gamma_0$, but here we want to impose Dirichlet boundary conditions in the closure of $\Gamma_0$, i.e., also on the vertices on the contour of $\Gamma_0$. Fortunately, the objects on the contour of $\Gamma_0$ are identified  with the tag `"sides_c"` (see next figure). Thus, the Dirichlet boundary $\Gamma_0$ can be built as the union of the objects identified as `"sides"` and `"sides_c"`.
 #
@@ -134,7 +135,7 @@ uh, = solve!(uh0,solver,op)
 
 # We finish this tutorial by writing the computed solution for visualization (see next figure).
 
-writevtk(Ω,"results",cellfields=["uh"=>uh])
+writevtk(Ω,"output_path/results",cellfields=["uh"=>uh])
 
 # ![](../assets/p_laplacian/sol-plap.png)
 #
diff --git a/src/poisson.jl b/src/poisson.jl
index d06a22e..7d882bf 100644
--- a/src/poisson.jl
+++ b/src/poisson.jl
@@ -54,7 +54,8 @@ model = DiscreteModelFromFile("../models/model.json")
 #
 # You can easily inspect the generated discrete model in [Paraview](https://www.paraview.org/) by writing it in `vtk` format.
 
-writevtk(model,"model")
+mkpath("output_path")
+writevtk(model,"output_path/model")
 
 # The previous line generates four different files `model_0.vtu`, `model_1.vtu`, `model_2.vtu`, and `model_3.vtu` containing the vertices, edges, faces, and cells present in the discrete model. Moreover, you can easily inspect which boundaries are defined within the model.
 #
@@ -137,7 +138,7 @@ uh = solve(solver,op)
 
 # The `solve` function returns the computed numerical solution `uh`. This object is an instance of `FEFunction`, the type used to represent a function in a FE space. We can inspect the result by writing it into a `vtk` file:
 
-writevtk(Ω,"results",cellfields=["uh"=>uh])
+writevtk(Ω,"output_path/results",cellfields=["uh"=>uh])
 
 #  which will generate a file named `results.vtu` having a nodal field named `"uh"` containing the solution of our problem (see next figure).
 #
diff --git a/src/poisson_amr.jl b/src/poisson_amr.jl
index 564aab2..f86d5f9 100644
--- a/src/poisson_amr.jl
+++ b/src/poisson_amr.jl
@@ -164,6 +164,7 @@ end
 nsteps = 5
 order = 1
 model = LShapedModel(10)
+mkpath("output_path")
 
 last_error = Inf
 for i in 1:nsteps
@@ -173,7 +174,7 @@ for i in 1:nsteps
   
   Ω = Triangulation(model)
   writevtk(
-    Ω,"model_$(i-1)",append=false,
+    Ω,"output_path/model_$(i-1)",append=false,
     cellfields = [
       "uh" => uh,                    # Computed solution
       "η" => CellField(η,Ω),        # Error indicators
diff --git a/src/poisson_distributed.jl b/src/poisson_distributed.jl
index 6d4549f..22dbbb5 100644
--- a/src/poisson_distributed.jl
+++ b/src/poisson_distributed.jl
@@ -33,11 +33,12 @@ function main_ex1(rank_partition,distribute)
   l(v) = ∫( v*f )dΩ
   op = AffineFEOperator(a,l,U,V)
   uh = solve(op)
-  writevtk(Ω,"results_ex1",cellfields=["uh"=>uh,"grad_uh"=>∇(uh)])
+  writevtk(Ω,"output_path/results_ex1",cellfields=["uh"=>uh,"grad_uh"=>∇(uh)])
 end
 
 # Once the `main_ex1` function has been defined, we have to trigger its execution on the different parts. To this end, one calls the `with_mpi` function of [`PartitionedArrays.jl`](https://github.com/fverdugo/PartitionedArrays.jl) right at the beginning of the program.
 
+mkpath("output_path")
 rank_partition = (2,2)
 with_mpi() do distribute
   main_ex1(rank_partition,distribute)
@@ -73,7 +74,7 @@ function main_ex2(rank_partition,distribute)
     op = AffineFEOperator(a,l,U,V)
     solver = PETScLinearSolver()
     uh = solve(solver,op)
-    writevtk(Ω,"results_ex2",cellfields=["uh"=>uh,"grad_uh"=>∇(uh)])
+    writevtk(Ω,"output_path/results_ex2",cellfields=["uh"=>uh,"grad_uh"=>∇(uh)])
   end
 end
 
@@ -112,7 +113,7 @@ function main_ex3(nparts,distribute)
     op = AffineFEOperator(a,l,U,V)
     solver = PETScLinearSolver()
     uh = solve(solver,op)
-    writevtk(Ω,"results_ex3",cellfields=["uh"=>uh,"grad_uh"=>∇(uh)])
+    writevtk(Ω,"output_path/results_ex3",cellfields=["uh"=>uh,"grad_uh"=>∇(uh)])
   end
 end
 
@@ -144,7 +145,7 @@ function main_ex4(nparts,distribute)
     op = AffineFEOperator(a,l,U,V)
     solver = PETScLinearSolver()
     uh = solve(solver,op)
-    writevtk(Ω,"results_ex4",cellfields=["uh"=>uh,"grad_uh"=>∇(uh)])
+    writevtk(Ω,"output_path/results_ex4",cellfields=["uh"=>uh,"grad_uh"=>∇(uh)])
   end
 end
 
diff --git a/src/poisson_hdg.jl b/src/poisson_hdg.jl
index aa2b9b5..97593ea 100644
--- a/src/poisson_hdg.jl
+++ b/src/poisson_hdg.jl
@@ -154,7 +154,8 @@ dΩ = Measure(Ω,degree)
 eh = uh - u
 l2_uh = sqrt(sum(∫(eh⋅eh)*dΩ))
 
-writevtk(Ω,"results",cellfields=["uh"=>uh,"qh"=>qh,"eh"=>eh])
+mkpath("output_path")
+writevtk(Ω,"output_path/results",cellfields=["uh"=>uh,"qh"=>qh,"eh"=>eh])
 
 # ## Going Further
 #
diff --git a/src/poisson_unfitted.jl b/src/poisson_unfitted.jl
index e2525cd..371cec1 100644
--- a/src/poisson_unfitted.jl
+++ b/src/poisson_unfitted.jl
@@ -104,8 +104,9 @@ cutgeo = cut(bgmodel,geo3)
 # To illustrate this concept, we can plot both the background and active triangulations to compare them.
 
 Ω_bg = Triangulation(bgmodel)
-writevtk(Ω_bg,"bg_trian")
-writevtk(Ω_act,"act_trian")
+mkpath("output_path")
+writevtk(Ω_bg,"output_path/bg_trian")
+writevtk(Ω_act,"output_path/act_trian")
 
 # In the picture below of the background grid, white cells are _inactive_, whereas gray cells are _active_.
 #
@@ -116,7 +117,7 @@ writevtk(Ω_act,"act_trian")
 # An `EmbeddedDiscretization` instance (here, `cutgeo`) also generates subtriangulations on each cut cells to represent the portion of the cell which is inside the domain of analysis. We use these subtriangulations to generate the so called _physical_ triangulations. Physical triangulations are nothing other than a body-fitted mesh of our domain $\Omega$, but _we only use them to integrate the weak form_ of the problem in $\Omega$, we won't define FE spaces and assign DoFs on top of them. In [GridapEmbedded](https://github.com/gridap/GridapEmbedded.jl) we build physical triangulations using the `PHYSICAL` keyword.
 
 Ω = Triangulation(cutgeo,PHYSICAL)
-writevtk(Ω,"phys_trian")
+writevtk(Ω,"output_path/phys_trian")
 
 # Once again, we can combine plots of the physical and active triangulations to illustrate these concepts. In the first plot, we show the physical triangulation within the background one. 
 #
@@ -172,7 +173,7 @@ aggregates = aggregate(strategy,cutgeo)
 
 colors = color_aggregates(aggregates,bgmodel)
 Ω_bg = Triangulation(bgmodel)
-writevtk(Ω_bg,"aggs_on_bg_trian",celldata=["aggregate"=>aggregates,"color"=>colors])
+writevtk(Ω_bg,"output_path/aggs_on_bg_trian",celldata=["aggregate"=>aggregates,"color"=>colors])
 
 # Finally, we use the aggregates to constrain the exterior DoFs of `Vstd` in terms of the interior ones. This leads to the AgFEM space.
 
@@ -250,7 +251,7 @@ using Test
 @test el2/ul2 < 1.e-8
 @test eh1/uh1 < 1.e-7
 
-writevtk(Ω,"results.vtu",cellfields=["uh"=>uh])
+writevtk(Ω,"output_path/results.vtu",cellfields=["uh"=>uh])
 
 # ![fig10](../assets/unfitted_poisson/fig_results.png)
 
diff --git a/src/stokes.jl b/src/stokes.jl
index 15dd61a..81f4228 100644
--- a/src/stokes.jl
+++ b/src/stokes.jl
@@ -53,4 +53,5 @@ op = AffineFEOperator(a,l,X,Y)
 uh, ph = solve(op)
 
 # Export results to vtk
-writevtk(Ωₕ,"results",order=2,cellfields=["uh"=>uh,"ph"=>ph])
+mkpath("output_path")
+writevtk(Ωₕ,"output_path/results",order=2,cellfields=["uh"=>uh,"ph"=>ph])
diff --git a/src/transient_linear.jl b/src/transient_linear.jl
index 1b70d93..8291809 100644
--- a/src/transient_linear.jl
+++ b/src/transient_linear.jl
@@ -108,14 +108,12 @@ uh = solve(solver, op, t0, tF, uh0)
 
 # We highlight that `uh` is an iterable function and the result at each time step is only computed lazily when iterating over it. We can post-process the results and generate the corresponding `vtk` files using the `createpvd` and `createvtk` functions. The former will create a `.pvd` file with the collection of `.vtu` files saved at each time step by `createvtk`. The computation of the problem solutions will be triggered in the following loop:
 
-if !isdir("tmp")
-  mkdir("tmp")
-end
+mkpath("output_path/results")
 
-createpvd("results") do pvd
-  pvd[0] = createvtk(Ω, "tmp/results_0" * ".vtu", cellfields=["u" => uh0])
+createpvd("output_path/results") do pvd
+  pvd[0] = createvtk(Ω, "output_path/results/results_0" * ".vtu", cellfields=["u" => uh0])
   for (tn, uhn) in uh
-    pvd[tn] = createvtk(Ω, "tmp/results_$tn" * ".vtu", cellfields=["u" => uhn])
+    pvd[tn] = createvtk(Ω, "output_path/results/results_$tn" * ".vtu", cellfields=["u" => uhn])
   end
 end
 
diff --git a/src/transient_nonlinear.jl b/src/transient_nonlinear.jl
index f7baeb2..14430c8 100644
--- a/src/transient_nonlinear.jl
+++ b/src/transient_nonlinear.jl
@@ -126,14 +126,12 @@ uh = solve(solver, op_sl, t0, tF, uh0)
 
 # Here again, we export the solution at each time step as follows
 
-if !isdir("tmp_nl")
-  mkdir("tmp_nl")
-end
+mkpath("output_path/results")
 
-createpvd("results_nl") do pvd
-  pvd[0] = createvtk(Ω, "tmp_nl/results_0" * ".vtu", cellfields=["u" => uh0])
+createpvd("output_path/results_nl") do pvd
+  pvd[0] = createvtk(Ω, "output_path/results/results_0" * ".vtu", cellfields=["u" => uh0])
   for (tn, uhn) in uh
-    pvd[tn] = createvtk(Ω, "tmp_nl/results_$tn" * ".vtu", cellfields=["u" => uhn])
+    pvd[tn] = createvtk(Ω, "output_path/results/results_$tn" * ".vtu", cellfields=["u" => uhn])
   end
 end
 
diff --git a/src/validation.jl b/src/validation.jl
index 5c52be9..96e7f9c 100644
--- a/src/validation.jl
+++ b/src/validation.jl
@@ -66,7 +66,8 @@ model3d = CartesianDiscreteModel(domain3d,partition3d)
 
 # Let us return to the 2D `CartesianDiscreteModel` that we have already constructed. You can inspect it by writing it into vtk format. Note that you can also print a 3D model, but not a 4D one. In the future, it would be cool to generate a movie from a 4D model, but this functionality is not yet implemented.
 
-writevtk(model,"model")
+mkpath("output_path")
+writevtk(model,"output_path/model")
 
 
 # If you open the generated files, you will see that the boundary vertices and facets are identified with the name "boundary". This is just what we need to impose the Dirichlet boundary conditions in this example.
@@ -113,7 +114,7 @@ e = u - uh
 
 # Once the error is defined, you can, e.g., visualize it.
 
-writevtk(Ω,"error",cellfields=["e" => e])
+writevtk(Ω,"output_path/error",cellfields=["e" => e])
 
 # This generates a file called `error.vtu`. Open it with Paraview to check that the error is of the order of the machine precision.
 #
diff --git a/src/validation_DrWatson.jl b/src/validation_DrWatson.jl
index 92a4b40..9981197 100644
--- a/src/validation_DrWatson.jl
+++ b/src/validation_DrWatson.jl
@@ -117,7 +117,7 @@ end
 
 function run_or_load(case::Dict)
   produce_or_load(
-    projectdir("assets","validation_DrWatson"),
+    projectdir("notebooks","output_path"),
     case,
     run,
     prefix="res",
@@ -129,7 +129,7 @@ end
 
 map(run_or_load,dicts)
 
-# Note that the results of each case are stored in a binary database file in the `projectdir("assets","validation_DrWatson")` folder. Each result file stores the output dictionary that returns from `run(case)`.
+# Note that the results of each case are stored in a binary database file in the `projectdir("notebooks","output_path")` folder. Each result file stores the output dictionary that returns from `run(case)`.
 
 # We also observe that we set `tag=true` in `produce_or_load`. This option is *key to preserve reproducibility*. It adds to the output dictionary the field `:gitcommit`, thus allowing us to trace the status of the code, at which we obtained those results. Furthermore, if the git repo is dirty, one more field `:gitpatch` is added, storing the difference string.
 
@@ -143,7 +143,7 @@ using DataFrames
 
 # To collect all simulation results, it suffices to use the `collect_results!` function from `DrWatson.jl` from the folder where the results are stored.
 
-df = collect_results(projectdir("assets","validation_DrWatson"))
+df = collect_results(projectdir("notebooks","output_path"))
 
 # We order next the database by (ascending) mesh size and we extract the arrays of mesh sizes and errors