Writing output of tutorials to a folder

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)

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)
 #

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
 

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)
 #

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).

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!

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),

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

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
 #

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
 

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)
 #

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
 

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)])
 

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])

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)
 #

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).
 #

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

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
 

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
 #