Merge pull request #49 from MultiSimOLab/sims-tests

jmartfrut · web-flow · commit 0bc8b3ba043a · 2025-11-06T10:22:43.000+01:00
Added more simulations tests
diff --git a/benchmark/Project.toml b/benchmark/Project.toml
@@ -3,3 +3,4 @@ BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 Gridap = "56d4f2e9-7ea1-5844-9cf6-b9c51ca7ce8e"
 GridapGmsh = "3025c34a-b394-11e9-2a55-3fee550c04c8"
 GridapSolvers = "6d3209ee-5e3c-4db7-a716-942eb12ed534"
+TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
diff --git a/benchmark/SimulationsBenchmarks/StaticMechanicalDirichletBenchmark.jl b/benchmark/SimulationsBenchmarks/StaticMechanicalDirichletBenchmark.jl
@@ -0,0 +1,5 @@
+
+filename = projdir("test/data/StaticMechanicalDirichletSimulation.jl")
+include(filename)
+
+SUITE["Simulations"]["StaticMechanicalDirichlet"] = @benchmarkable static_mechanical_dirichlet_simulation(writevtk=false, verbose=false)
diff --git a/benchmark/SimulationsBenchmarks/StaticMechanicalNeumannBenchmark.jl b/benchmark/SimulationsBenchmarks/StaticMechanicalNeumannBenchmark.jl
@@ -0,0 +1,5 @@
+
+filename = projdir("test/data/StaticMechanicalNeumannSimulation.jl")
+include(filename)
+
+SUITE["Simulations"]["StaticMechanicalNeumann"] = @benchmarkable static_mechanical_neumann_simulation(writevtk=false, verbose=false)
diff --git a/benchmark/SimulationsBenchmarks/ViscoElasticSimulationBenchmark.jl b/benchmark/SimulationsBenchmarks/ViscoElasticSimulationBenchmark.jl
@@ -1,6 +1,5 @@
 
-BASE_FOLDER = dirname(dirname(pathof(HyperFEM)))
-filename = joinpath(BASE_FOLDER, "test/data/ViscoElasticSimulation.jl")
+filename = projdir("test/data/ViscoElasticSimulation.jl")
 include(filename)
 
-SUITE["Simulations"]["ViscoElastic"] = @benchmarkable visco_elastic_simulation(write_vtk=false)
+SUITE["Simulations"]["ViscoElastic"] = @benchmarkable visco_elastic_simulation(writevtk=false, verbose=false)
diff --git a/benchmark/SimulationsBenchmarks/benchmarks.jl b/benchmark/SimulationsBenchmarks/benchmarks.jl
@@ -2,3 +2,7 @@
 SUITE["Simulations"] = BenchmarkGroup()
 
 include("ViscoElasticSimulationBenchmark.jl")
+
+include("StaticMechanicalDirichletBenchmark.jl")
+
+include("StaticMechanicalNeumannBenchmark.jl")
diff --git a/src/Exports.jl b/src/Exports.jl
@@ -136,8 +136,6 @@ end
 @publish Solvers Roman_LS
 @publish Solvers update_cellstate!
 
- 
-
-
-# @publish LinearSolvers solve
-# @publish LinearSolvers solve!
+export setupfolder
+export projdir
+export filedir
diff --git a/src/HyperFEM.jl b/src/HyperFEM.jl
@@ -1,16 +1,12 @@
 module HyperFEM
 
-using TimerOutputs
-
-
 include("TensorAlgebra/TensorAlgebra.jl")
 include("PhysicalModels/PhysicalModels.jl")
 include("WeakForms/WeakForms.jl")
 include("Solvers/Solvers.jl")
 include("ComputationalModels/ComputationalModels.jl")
 
 include("Io.jl")
-export setupfolder
-
 include("Exports.jl")
+
 end
diff --git a/src/Io.jl b/src/Io.jl
@@ -31,3 +31,15 @@ end
 
 function cleandir(::String, ::Nothing)
 end
+
+"""
+Return the path to the specified folders relative to the HyperFEM path.
+
+# Examples
+    folder = projdir("data", "sims")
+    folder = projdir("test/data/mesh.msh")
+"""
+function projdir(folders::String...)
+  base_folder = dirname(dirname(pathof(HyperFEM)))
+  joinpath(base_folder, folders...)
+end
diff --git a/test/Project.toml b/test/Project.toml
@@ -5,3 +5,4 @@ GridapSolvers = "6d3209ee-5e3c-4db7-a716-942eb12ed534"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
diff --git a/test/SimulationsTests/StaticMechanicalDirichletTest.jl b/test/SimulationsTests/StaticMechanicalDirichletTest.jl
@@ -0,0 +1,7 @@
+
+filename = projdir("test/data/StaticMechanicalDirichletSimulation.jl")
+include(filename)
+
+x = static_mechanical_dirichlet_simulation(writevtk=false, verbose=false)
+
+@test norm(x) ≈ 0.27148722276
diff --git a/test/SimulationsTests/StaticMechanicalNeumannTest.jl b/test/SimulationsTests/StaticMechanicalNeumannTest.jl
@@ -0,0 +1,7 @@
+
+filename = projdir("test/data/StaticMechanicalNeumannSimulation.jl")
+include(filename)
+
+x = static_mechanical_neumann_simulation(writevtk=false, verbose=false)
+
+@test norm(x) ≈ 1.000148588846
diff --git a/test/SimulationsTests/ViscoElasticSimulationTest.jl b/test/SimulationsTests/ViscoElasticSimulationTest.jl
@@ -1,8 +1,7 @@
 
-BASE_FOLDER = dirname(dirname(pathof(HyperFEM)))
-filename = joinpath(BASE_FOLDER, "test/data/ViscoElasticSimulation.jl")
+filename = projdir("test/data/ViscoElasticSimulation.jl")
 include(filename)
 
-λx, σΓ = visco_elastic_simulation(t_end=5, write_vtk=false, verbose=false)
+λx, σΓ = visco_elastic_simulation(t_end=5, writevtk=false, verbose=false)
 
 @test σΓ[end] ≈ 152821.386
diff --git a/test/SimulationsTests/runtests.jl b/test/SimulationsTests/runtests.jl
@@ -6,4 +6,8 @@ using Test
    
     include("ViscoElasticSimulationTest.jl")
 
+    include("StaticMechanicalDirichletTest.jl")
+
+    include("StaticMechanicalNeumannTest.jl")
+
 end
diff --git a/test/TestConstitutiveModels/PhysicalModelTests.jl b/test/TestConstitutiveModels/PhysicalModelTests.jl
@@ -527,8 +527,7 @@ end
     return n_layers, Weights, Biases, ϵ, β
   end
 
-  BASE_FOLDER = dirname(dirname(pathof(HyperFEM)))
-  data_filename = joinpath(BASE_FOLDER, "test/models/test_NN_TEM.json")
+  data_filename = projdir("test/models/test_NN_TEM.json")
   n_layers, Weights, Biases, ϵ, β = ExtractingInfo(data_filename)
 
   model = ThermoElectroMech_PINNs(; W=Weights, b=Biases, ϵ=ϵ, β=β, nLayer=n_layers)
diff --git a/test/data/StaticMechanicalDirichletSimulation.jl b/test/data/StaticMechanicalDirichletSimulation.jl
@@ -0,0 +1,85 @@
+using Gridap, GridapSolvers
+using GridapSolvers.NonlinearSolvers
+using GridapSolvers.LinearSolvers
+using TimerOutputs
+using Gridap.FESpaces
+using HyperFEM
+using HyperFEM.ComputationalModels.CartesianTags
+
+
+function static_mechanical_dirichlet_simulation(;writevtk=true, verbose=true)
+
+  pname = "Stretch"
+  simdir = projdir("data", "sims", pname)
+  setupfolder(simdir)
+
+  long   = 0.05   # m
+  width  = 0.005  # m
+  thick  = 0.002  # m
+  geometry = CartesianDiscreteModel((0, long, 0, width, 0, thick), (5,2,2))
+  labels = get_face_labeling(geometry)
+  add_tag_from_tags!(labels, "fixed", CartesianTags.faceX0)
+  add_tag_from_tags!(labels, "moving", CartesianTags.faceX1)
+
+  physmodel = MooneyRivlin3D(λ=3.0, μ1=1.0, μ2=0.0, ρ=1.0)
+
+  # Setup integration
+  order = 1
+  degree = 2 * order
+  Ω = Triangulation(geometry)
+  dΩ = Measure(Ω, degree)
+
+  # Dirichlet boundary conditions
+  dir_u_tags = ["fixed", "moving"]
+  dir_u_values = [[0.0, 0.0, 0.0], [0.08, 0.0, 0.0]]
+  dir_u_timesteps = [Λ -> 1.0, Λ -> Λ]
+  D_bc = DirichletBC(dir_u_tags, dir_u_values, dir_u_timesteps)
+
+  #  FE spaces
+  reffeu = ReferenceFE(lagrangian, VectorValue{3,Float64}, order)
+
+  V = TestFESpace(Ω, reffeu, D_bc, conformity=:H1)
+  U = TrialFESpace(V, D_bc, 0.0)
+
+  #  residual and jacobian function of load factor
+  k = Kinematics(Mechano, Solid)
+  res(Λ) = (u, v) -> residual(physmodel, k, u, v, dΩ)
+  jac(Λ) = (u, du, v) -> jacobian(physmodel, k, u, du, v, dΩ)
+
+  #computational model
+  ls = LUSolver()
+  nls = NewtonSolver(ls; maxiter=15, rtol=1.e-12, verbose=verbose)
+
+  comp_model = StaticNonlinearModel(res, jac, U, V, D_bc; nls=nls)
+
+  function driverpost_mech(post)
+    if writevtk
+      state = post.comp_model.caches[3]
+      Λ_ = post.iter
+      Λ = post.Λ[Λ_]
+      xh = FEFunction(U, state)
+      pvd = post.cachevtk[3]
+      filePath = post.cachevtk[2]
+      if post.cachevtk[1]  
+        Λstring = replace(string(round(Λ, digits=2)), "." => "_")
+        pvd[Λ_] = createvtk(Ω,
+          filePath *  "/" * pname * "_Λ_" * Λstring * ".vtu",
+          cellfields=["u" => xh])
+      end
+    end
+  end
+
+  post_model = PostProcessor(comp_model, driverpost_mech; is_vtk=writevtk, filepath=simdir)
+
+  @timeit pname begin
+    x, flag = solve!(comp_model; stepping=(nsteps=10, maxbisec=10), post=post_model)
+  end
+  return x
+end
+
+
+if abspath(PROGRAM_FILE) == @__FILE__
+  reset_timer!()
+  static_mechanical_dirichlet_simulation()
+  print_timer()
+end
diff --git a/test/data/StaticMechanicalNeumannSimulation.jl b/test/data/StaticMechanicalNeumannSimulation.jl
@@ -0,0 +1,90 @@
+using Gridap, GridapSolvers
+using GridapSolvers.NonlinearSolvers
+using GridapSolvers.LinearSolvers
+using TimerOutputs
+using Gridap.FESpaces
+using HyperFEM
+using HyperFEM.ComputationalModels.CartesianTags
+
+
+function static_mechanical_neumann_simulation(;writevtk=true, verbose=true)
+
+  pname = "StaticMechanical"
+  simdir = projdir("data", "sims", pname)
+  setupfolder(simdir)
+
+  long   = 0.05   # m
+  width  = 0.005  # m
+  thick  = 0.002  # m
+  geometry = CartesianDiscreteModel((0, long, 0, width, 0, thick), (5,2,2))
+  labels = get_face_labeling(geometry)
+  add_tag_from_tags!(labels, "fixed", CartesianTags.faceX0)
+  add_tag_from_tags!(labels, "force", CartesianTags.faceX1)
+  
+  physmodel = MooneyRivlin3D(λ=3.0, μ1=1.0, μ2=0.0, ρ=1.0)
+
+  # Setup integration
+  order = 1
+  degree = 2 * order + 1
+  Ω = Triangulation(geometry)
+  dΩ = Measure(Ω, degree)
+
+  # Dirichlet conditions 
+  dir_u_tags = ["fixed"]
+  dir_u_values = [[0.0, 0.0, 0.0]]
+  dir_u_timesteps = [Λ -> 1.0]
+  D_bc = DirichletBC(dir_u_tags, dir_u_values, dir_u_timesteps)
+
+  # Neumann conditions 
+  neu_F_tags = ["force"]
+  neu_F_values = [[0.0, 0.0, -1e-3]]
+  neu_F_timesteps = [Λ -> Λ]
+  N_bc = NeumannBC(neu_F_tags, neu_F_values, neu_F_timesteps)
+  dΓ = get_Neumann_dΓ(geometry, N_bc, degree)
+
+  #  FE spaces
+  reffeu = ReferenceFE(lagrangian, VectorValue{3,Float64}, order)
+  V = TestFESpace(Ω, reffeu, D_bc, conformity=:H1)
+  U = TrialFESpace(V, D_bc, 0.0)
+
+  #  residual and jacobian function of load factor
+  k = Kinematics(Mechano, Solid)
+  res(Λ) = (u, v) -> residual(physmodel, k, u, v, dΩ) + residual_Neumann(N_bc, v, dΓ, Λ)
+  jac(Λ) = (u, du, v) -> jacobian(physmodel, k, u, du, v, dΩ)
+
+  ls = LUSolver()
+  nls = NewtonSolver(ls; maxiter=20, atol=1.e-10, rtol=1.e-8, verbose=verbose)
+
+  comp_model = StaticNonlinearModel(res, jac, U, V, D_bc; nls=nls)
+
+  function driverpost_mech(post)
+    if writevtk
+      state = post.comp_model.caches[3]
+      Λ_ = post.iter
+      Λ = post.Λ[Λ_]
+      xh = FEFunction(U, state)
+      pvd = post.cachevtk[3]
+      filePath = post.cachevtk[2]
+      if post.cachevtk[1]
+        Λstring = replace(string(round(Λ, digits=2)), "." => "_")
+        pvd[Λ_] = createvtk(Ω,
+          filePath * "/" * pname * "_Λ_" * Λstring * ".vtu",
+          cellfields=["u" => xh])
+      end
+    end
+  end
+
+  post_model = PostProcessor(comp_model, driverpost_mech; is_vtk=writevtk, filepath=simdir)
+
+  @timeit pname begin
+    x = solve!(comp_model; stepping=(nsteps=8, maxbisec=0), post=post_model)
+  end
+  return x
+end
+
+
+if abspath(PROGRAM_FILE) == @__FILE__
+  reset_timer!()
+  static_mechanical_neumann_simulation()
+  print_timer()
+end
diff --git a/test/data/ViscoElasticSimulation.jl b/test/data/ViscoElasticSimulation.jl
@@ -8,7 +8,7 @@ using HyperFEM.ComputationalModels:constant
 using HyperFEM.ComputationalModels:triangular
 using HyperFEM.ComputationalModels.PostMetrics
 
-function visco_elastic_simulation(;t_end=15, write_vtk=true, verbose=true)
+function visco_elastic_simulation(;t_end=15, writevtk=true, verbose=true)
   # Domain and tessellation
   long   = 0.05   # m
   width  = 0.005  # m
@@ -81,7 +81,7 @@ function visco_elastic_simulation(;t_end=15, write_vtk=true, verbose=true)
     updateStateVariables!(state_vars, cons_model, Δt, F∘(∇(uh)'), F∘(∇(unh)'))
   end
 
-  post_model = PostProcessor(comp_model, driverpost; is_vtk=false, filepath="")
+  post_model = PostProcessor(comp_model, driverpost; is_vtk=writevtk, filepath="")
   solve!(comp_model; stepping=(nsteps=Int(t_end/Δt), maxbisec=1), post=post_model, ProjectDirichlet=true)
   (λx, σΓ)
 end
