tutorials 14 to 27

Author: miguelmaso

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

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)
 

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