From e39b58000e4a3409c0f068d3f0fdb3cd2b08fe25 Mon Sep 17 00:00:00 2001
From: JordiManyer <jordi.manyer@monash.edu>
Date: Mon, 11 Nov 2024 18:06:08 +1100
Subject: [PATCH] Added printing for tests

---
 test/seq/EmbeddedTests.jl | 67 ++++++++++++++++++++++++++++++++++++---
 1 file changed, 63 insertions(+), 4 deletions(-)

diff --git a/test/seq/EmbeddedTests.jl b/test/seq/EmbeddedTests.jl
index e1c373ea..c7990e86 100644
--- a/test/seq/EmbeddedTests.jl
+++ b/test/seq/EmbeddedTests.jl
@@ -6,7 +6,7 @@ using Gridap, Gridap.Geometry, Gridap.Adaptivity
 using GridapEmbedded, GridapEmbedded.LevelSetCutters
 
 using Gridap.Arrays: Operation
-using GridapTopOpt: get_conormal_vector
+using GridapTopOpt: get_conormal_vector,get_subfacet_normal_vector,get_ghost_normal_vector
 
 function generate_model(D,n)
   domain = (D==2) ? (0,1,0,1) : (0,1,0,1,0,1)
@@ -41,8 +41,11 @@ function level_set(shape::Symbol)
   end
 end
 
-function main(model,φ::Function,f::Function)
-
+function main(
+  model,φ::Function,f::Function;
+  vtk=false,
+  name="embedded"
+)
   order = 1
   reffe = ReferenceFE(lagrangian,Float64,order)
   V_φ = TestFESpace(model,reffe)
@@ -109,6 +112,62 @@ function main(model,φ::Function,f::Function)
   dJ_int_exact_vec = assemble_vector(dJ_int_exact,V_φ)
 
   @test norm(dJ_int_AD_vec - dJ_int_exact_vec) < 1e-10
+
+  if vtk
+    path = "results/$(name)"
+    Ω_bg = Triangulation(model)
+    writevtk(
+      Ω_bg,"$(path)_results",
+      cellfields = [
+        "φh" => φh,"∇φh" => ∇(φh),
+        "dJ_bulk_AD" => FEFunction(V_φ,dJ_bulk_AD_vec),
+        "dJ_bulk_exact" => FEFunction(V_φ,dJ_bulk_exact_vec),
+        "dJ_int_AD" => FEFunction(V_φ,dJ_int_AD_vec),
+        "dJ_int_exact" => FEFunction(V_φ,dJ_int_exact_vec)
+      ],
+      celldata = [
+        "inoutcut" => GridapEmbedded.Interfaces.compute_bgcell_to_inoutcut(model,geo)
+      ];
+      append = false
+    )
+  
+    writevtk(
+      Ω, "$(path)_bulk"; append = false
+    )
+    writevtk(
+      Γ, "$(path)_boundary"; append = false
+    )
+
+    n_∂Ω_Λ = get_subfacet_normal_vector(Λ)
+    n_k_Λ  = get_ghost_normal_vector(Λ)
+    writevtk(
+      Λ, "$(path)_skeleton",
+      cellfields = [
+        "n_∂Ω.plus" => n_∂Ω_Λ.plus,"n_∂Ω.minus" => n_∂Ω_Λ.minus,
+        "n_k.plus" => n_k_Λ.plus,"n_k.minus" => n_k_Λ.minus,
+        "n_S" => n_S_Λ,
+        "m_k.plus" => m_k_Λ.plus,"m_k.minus" => m_k_Λ.minus,
+        "∇ˢφ" => ∇ˢφ_Λ,
+        "∇φh_Γs_plus" => ∇(φh).plus,"∇φh_Γs_minus" => ∇(φh).minus,
+        "jump(fh*m_k)" => jump(fh*m_k_Λ)
+      ];
+      append = false
+    )
+
+    if num_cells(Σ) > 0
+      n_∂Ω_Σ = get_subfacet_normal_vector(Σ)
+      n_k_Σ  = get_ghost_normal_vector(Σ)
+      writevtk(
+        Σ, "$(path)_skeleton",
+        cellfields = [
+          "n_∂Ω" => n_∂Ω_Σ, "n_k" => n_k_Σ,
+          "n_S" => n_S_Σ, "m_k" => m_k_Σ,
+          "∇ˢφ" => ∇ˢφ_Σ, "∇φh_Γs" => ∇(φh),
+        ];
+        append = false
+      )
+    end
+  end
 end
 
 D = 2
@@ -116,6 +175,6 @@ n = 10
 model = generate_model(D,n)
 φ = level_set(:circle)
 f = x -> 1.0
-main(model,φ,f)
+main(model,φ,f;vtk=true)
 
 end
\ No newline at end of file