Post-processing reactions, testing, benchmarking

README.md

@@ -9,7 +9,7 @@ src="https://github.com/jmartfrut/HyperFEM/blob/main/docs/imgs/logo.png?raw=true
 
 [![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://jmartfrut.github.io/HyperFEM.jl/stable/)
 [![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://jmartfrut.github.io/HyperFEM.jl/dev/)
-[![Build Status](https://github.com/jmartfrut/HyperFEM.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/jmartfrut/HyperFEM.jl/actions/workflows/CI.yml?query=branch%3Amain)
+[![Build Status](https://github.com/MultiSimOLab/HyperFEM/actions/workflows/ci.yml/badge.svg?branch=main)](https://github.com/MultiSimOLab/HyperFEM/actions/workflows/ci.yml?branch=main)
 [![Coverage](https://codecov.io/gh/jmartfrut/HyperFEM.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/jmartfrut/HyperFEM.jl)
 
 # **M**ultiphysics-informed **D**esign of **T**unable **S**mart **M**aterials
@@ -48,7 +48,7 @@ geomodel = GmshDiscreteModel("./test/models/test_static_EM.msh")
 # Constitutive model
 physmodel_mec = NeoHookean3D(λ=10.0, μ=1.0)
 physmodel_elec = IdealDielectric(ε=1.0)
-physmodel = ElectroMechModel(Mechano=physmodel_mec, Electro=physmodel_elec)
+physmodel = ElectroMechModel(mechano=physmodel_mec, electro=physmodel_elec)
 
 # Setup integration
 order = 1

benchmark/SimulationsBenchmarks/ViscoElasticSimulationBenchmark.jl

@@ -0,0 +1,6 @@
+
+BASE_FOLDER = dirname(dirname(pathof(HyperFEM)))
+filename = joinpath(BASE_FOLDER, "test/data/ViscoElasticSimulation.jl")
+include(filename)
+
+SUITE["Simulations"]["ViscoElastic"] = @benchmarkable visco_elastic_simulation(write_vtk=false)

benchmark/SimulationsBenchmarks/benchmarks.jl

@@ -0,0 +1,4 @@
+
+SUITE["Simulations"] = BenchmarkGroup()
+
+include("ViscoElasticSimulationBenchmark.jl")

benchmark/benchmarks.jl

@@ -8,4 +8,4 @@ include("TensorAlgebraBenchmarks/benchmarks.jl")
 
 include("ConstitutiveModelsBenchmark/benchmarks.jl")
 
-include("Simulations/benchmarks.jl")
+include("SimulationsBenchmarks/benchmarks.jl")

src/ComputationalModels/ComputationalModels.jl

@@ -87,4 +87,8 @@ export Jacobian
 export Entropy
 export D0
 export reset!
+export interpolate_L2_tensor
+export interpolate_L2_vector
+export interpolate_L2_scalar
+
 end

src/ComputationalModels/PostMetrics.jl

@@ -3,6 +3,9 @@ module PostMetrics
 using Gridap
 using HyperFEM.PhysicalModels
 
+export component_LInf
+export volume_diff
+
 """
     component_LInf(::FEFunction, ::Symbol, ::Triangulation)::Float64
 

src/ComputationalModels/PostProcessors.jl

@@ -102,37 +102,28 @@ end
 
 
 function Cauchy(model::Elasto, uh, unh, state_vars, Ω, dΩ, t, Δt)
-    _, ∂Ψu, _ = model(t)
-    F, _, _ = get_Kinematics(model.Kinematic)
-    σ = ∂Ψu ∘ (F∘∇(uh))
-    return interpolate_σ_everywhere(σ, Ω, dΩ)
+    σh = Cauchy(model, uh)
+    interpolate_L2_tensor(σh, Ω, dΩ)
 end
 
 
 function Cauchy(model::ViscoElastic, uh, unh, state_vars, Ω, dΩ, t, Δt)
-    _, ∂Ψu, _ = model(t, Δt=Δt)
+    σh = Cauchy(model, uh, unh, state_vars, Δt)
+    interpolate_L2_tensor(σh, Ω, dΩ)
+end
+
+
+function Cauchy(model::Elasto, uh, vars...)
+    _, ∂Ψu, _ = model()
     F, _, _ = get_Kinematics(model.Kinematic)
-    σ = ∂Ψu ∘ (F∘∇(uh), F∘∇(unh), state_vars...)
-    return interpolate_σ_everywhere(σ, Ω, dΩ)
+    ∂Ψu ∘ (F∘∇(uh))
 end
 
 
-function interpolate_σ_everywhere(σ, Ω, dΩ)
-    ref_L2 = ReferenceFE(lagrangian, Float64, 0)
-    ref_fe = ReferenceFE(lagrangian, Float64, 1)
-    VL2 = FESpace(Ω, ref_L2, conformity=:L2)
-    V = FESpace(Ω, ref_fe, conformity=:H1)
-    n1 = VectorValue(1.0, 0.0, 0.0)
-    n2 = VectorValue(0.0, 1.0, 0.0)
-    n3 = VectorValue(0.0, 0.0, 1.0)
-    σ11h = interpolate_everywhere(L2_Projection(n1 ⋅ σ ⋅ n1, dΩ, VL2), V)
-    σ12h = interpolate_everywhere(L2_Projection(n1 ⋅ σ ⋅ n2, dΩ, VL2), V)
-    σ13h = interpolate_everywhere(L2_Projection(n1 ⋅ σ ⋅ n3, dΩ, VL2), V)
-    σ22h = interpolate_everywhere(L2_Projection(n2 ⋅ σ ⋅ n2, dΩ, VL2), V)
-    σ23h = interpolate_everywhere(L2_Projection(n2 ⋅ σ ⋅ n3, dΩ, VL2), V)
-    σ33h = interpolate_everywhere(L2_Projection(n3 ⋅ σ ⋅ n3, dΩ, VL2), V)
-    ph   = interpolate_everywhere(L2_Projection(tr ∘ σ, dΩ, VL2), V)
-    return (σ11h, σ12h, σ13h, σ22h, σ23h, σ33h, ph)
+function Cauchy(model::ViscoElastic, uh, unh, states, Δt)
+    _, ∂Ψu, _ = model(Δt=Δt)
+    F, _, _ = get_Kinematics(model.Kinematic)
+    ∂Ψu ∘ (F∘∇(uh), F∘∇(unh), states...)
 end
 
 
@@ -172,6 +163,49 @@ function D0(physmodel::ThermoElectroMechano, uh, φh, θh, Ω, dΩ, Λ=1.0)
 end
 
 
+function interpolate_L2_tensor(A, Ω, dΩ, Γ=Ω)
+    refL2 = ReferenceFE(lagrangian, Float64, 0)
+    reffe = ReferenceFE(lagrangian, Float64, 1)
+    VL2 = FESpace(Ω, refL2, conformity=:L2)
+    VH1 = FESpace(Γ, reffe, conformity=:H1)
+    n1 = VectorValue(1.0, 0.0, 0.0)
+    n2 = VectorValue(0.0, 1.0, 0.0)
+    n3 = VectorValue(0.0, 0.0, 1.0)
+    A11 = interpolate_everywhere(L2_Projection(n1 ⋅ A ⋅ n1, dΩ, VL2), VH1)
+    A12 = interpolate_everywhere(L2_Projection(n1 ⋅ A ⋅ n2, dΩ, VL2), VH1)
+    A13 = interpolate_everywhere(L2_Projection(n1 ⋅ A ⋅ n3, dΩ, VL2), VH1)
+    A22 = interpolate_everywhere(L2_Projection(n2 ⋅ A ⋅ n2, dΩ, VL2), VH1)
+    A23 = interpolate_everywhere(L2_Projection(n2 ⋅ A ⋅ n3, dΩ, VL2), VH1)
+    A33 = interpolate_everywhere(L2_Projection(n3 ⋅ A ⋅ n3, dΩ, VL2), VH1)
+    trA = interpolate_everywhere(L2_Projection(tr ∘ A,     dΩ, VL2), VH1)
+    (A11, A12, A13, A22, A23, A33, trA)
+end
+
+
+function interpolate_L2_vector(b, Ω, dΩ, Γ=Ω)
+    refL2 = ReferenceFE(lagrangian, Float64, 0)
+    reffe = ReferenceFE(lagrangian, Float64, 1)
+    VL2 = FESpace(Ω, refL2, conformity=:L2)
+    VH1 = FESpace(Γ, reffe, conformity=:H1)
+    n1 = VectorValue(1.0, 0.0, 0.0)
+    n2 = VectorValue(0.0, 1.0, 0.0)
+    n3 = VectorValue(0.0, 0.0, 1.0)
+    b1 = interpolate_everywhere(L2_Projection(n1 ⋅ b, dΩ, VL2), VH1)
+    b1 = interpolate_everywhere(L2_Projection(n2 ⋅ b, dΩ, VL2), VH1)
+    b1 = interpolate_everywhere(L2_Projection(n3 ⋅ b, dΩ, VL2), VH1)
+    (b1, b2, b3)
+end
+
+
+function interpolate_L2_scalar(x, Ω, dΩ, Γ=Ω)
+    refL2 = ReferenceFE(lagrangian, Float64, 0)
+    reffe = ReferenceFE(lagrangian, Float64, 1)
+    VL2 = FESpace(Ω, refL2, conformity=:L2)
+    VH1 = FESpace(Γ, reffe, conformity=:H1)
+    interpolate_everywhere(L2_Projection(x, dΩ, VL2), VH1)
+end
+
+
 function L2_Projection(u, dΩ, V)
     a(w, v) = ∫(w * v) * dΩ
     l(v)    = ∫(v * u) * dΩ

src/Exports.jl

@@ -111,12 +111,16 @@ end
 @publish ComputationalModels  GmshDiscreteModel
 @publish ComputationalModels  updateBC!
 @publish ComputationalModels  PostProcessor
+@publish ComputationalModels  PostMetrics
 @publish ComputationalModels  StaggeredModel
 @publish ComputationalModels  Cauchy
 @publish ComputationalModels  Jacobian
 @publish ComputationalModels  Entropy
 @publish ComputationalModels  D0
 @publish ComputationalModels  reset!
+@publish ComputationalModels  interpolate_L2_tensor
+@publish ComputationalModels  interpolate_L2_vector
+@publish ComputationalModels  interpolate_L2_scalar
 @publish ComputationalModels  DirichletCoupling
 @publish ComputationalModels  evaluate!
 @publish ComputationalModels  InterpolableBC

test/Project.toml

@@ -1,6 +1,7 @@
 [deps]
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Gridap = "56d4f2e9-7ea1-5844-9cf6-b9c51ca7ce8e"
+GridapSolvers = "6d3209ee-5e3c-4db7-a716-942eb12ed534"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

test/SimulationsTests/ViscoElasticSimulationTest.jl

@@ -0,0 +1,8 @@
+
+BASE_FOLDER = dirname(dirname(pathof(HyperFEM)))
+filename = joinpath(BASE_FOLDER, "test/data/ViscoElasticSimulation.jl")
+include(filename)
+
+λx, σΓ = visco_elastic_simulation(t_end=5, write_vtk=false, verbose=false)
+
+@test σΓ[end] ≈ 152821.386

test/SimulationsTests/runtests.jl

@@ -0,0 +1,9 @@
+using Gridap
+using HyperFEM
+using Test
+
+@testset "SimulationsTests" begin
+
+    include("ViscoElasticSimulationTest.jl")
+
+end

benchmark/Simulations/ViscoElasticBenchmark.jl → test/data/ViscoElasticSimulation.jl

@@ -6,9 +6,10 @@ using HyperFEM
 using HyperFEM.ComputationalModels.CartesianTags
 using HyperFEM.ComputationalModels:constant
 using HyperFEM.ComputationalModels:triangular
+using HyperFEM.ComputationalModels.PostMetrics
 
-function viscousbenchmark()
-  # Domain   and Tessellation
+function visco_elastic_simulation(;t_end=15, write_vtk=true, verbose=true)
+  # Domain and tessellation
   long   = 0.05   # m
   width  = 0.005  # m
   thick  = 0.001  # m
@@ -24,37 +25,38 @@ function viscousbenchmark()
   # Constitutive model
   μ = 1.367e4    # Pa
   λ = 1000μ      # Pa
-  hyper_elastic_model = NeoHookean3D(λ=λ, μ=μ)
   μv₁ = 3.153e5  # Pa
   τv₁ = 10.72    # s
+  hyper_elastic_model = NeoHookean3D(λ=λ, μ=μ)
   viscous_branch = ViscousIncompressible(IncompressibleNeoHookean3D(λ=0., μ=μv₁), τv₁)
   cons_model = GeneralizedMaxwell(hyper_elastic_model, viscous_branch)
 
+  # Dirichlet boundary conditions
+  strain = 0.5
+  D_bc = DirichletBC(
+    ["corner1", "corner2", "fixed", "moving"],
+    [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [long * strain, 0.0, 0.0]],
+    [constant(), constant(), constant(), triangular(10/t_end)])
+  dirichlet_masks = [
+    [true, true, true], [true, false, true], [true, false, false], [true, false, false]]
+
   # Setup integration
   order = 2
   degree = 2 * order
   Ω = Triangulation(model)
   dΩ = Measure(Ω, degree)
   dΓ = get_Neumann_dΓ(model, NothingBC(), degree)
+  Γ1  = BoundaryTriangulation(model, tags=D_bc.tags[4])
+  dΓ1 = Measure(Γ1, degree)
   Δt = 0.05
-  t_end = 5
-
-  # Dirichlet boundary conditions 
-  strain = 1.5
-  D_bc = DirichletBC(
-    ["corner1", "corner2", "fixed", "moving"],
-    [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [long * strain, 0.0, 0.0]],
-    [constant(), constant(), constant(), triangular(10/t_end)])
-  dirichlet_masks = [
-    [true, true, true], [true, false, true], [true, false, false], [true, false, false]]
 
   # FE spaces
   reffe = ReferenceFE(lagrangian, VectorValue{3,Float64}, order)
   Vu = TestFESpace(model, reffe, D_bc, dirichlet_masks=dirichlet_masks, conformity=:H1)
+  VL2 = FESpace(Ω, reffe, conformity=:L2)
   Uu = TrialFESpace(Vu, D_bc, 0.0)
   Uun = TrialFESpace(Vu, D_bc, 0.0)
 
-  # residual and jacobian
   uh = FEFunction(Uu, zero_free_values(Uu))
   unh = FEFunction(Uun, zero_free_values(Uun))
   state_vars = initializeStateVariables(cons_model, dΩ)
@@ -63,15 +65,27 @@ function viscousbenchmark()
   jac(Λ) = (u,du,v)->jacobian(cons_model, u, du, v, dΩ, t_end * Λ, Δt, unh, state_vars)
 
   ls = LUSolver()
-  nls_ = NewtonSolver(ls; maxiter=20, atol=1.e-6, rtol=1.e-6, verbose=true)
-  comp_model = StaticNonlinearModel(res, jac, Uu, Vu, D_bc; nls=nls_, xh=uh, xh⁻=unh)
+  nls = NewtonSolver(ls; maxiter=20, atol=1.e-6, rtol=1.e-6, verbose=verbose)
+  comp_model = StaticNonlinearModel(res, jac, Uu, Vu, D_bc; nls=nls, xh=uh, xh⁻=unh)
 
-  function driverpost(post; cons_model=cons_model, Δt=Δt, uh=uh, unh=unh, A=state_vars, Ω=Ω, dΩ=dΩ)
-    updateStateVariables!(A, cons_model, Δt, uh, unh)
+  λx = Float64[]
+  σΓ = Float64[]
+  function driverpost(post)
+    σh11, _... = Cauchy(cons_model, uh, unh, state_vars, Ω, dΩ, 0.0, Δt)
+    σΓ1 = sum(∫(σh11)dΓ1) / sum(∫(1.0)dΓ1)
+    push!(σΓ, σΓ1)
+    push!(λx, 1.0 + component_LInf(uh, :x, Ω) / long)
+    updateStateVariables!(state_vars, cons_model, Δt, uh, unh)
   end
-  post_model = PostProcessor(comp_model, driverpost; is_vtk=false, filepath="")
 
+  post_model = PostProcessor(comp_model, driverpost; is_vtk=false, filepath="")
   solve!(comp_model; stepping=(nsteps=Int(t_end/Δt), maxbisec=1), post=post_model, ProjectDirichlet=true)
+  (λx, σΓ)
 end
 
-SUITE["Simulations"]["ViscoElastic"] = @benchmarkable viscousbenchmark()
+
+if abspath(PROGRAM_FILE) == @__FILE__
+  using Plots
+  λx, σΓ = visco_elastic_simulation()
+  plot(λx, σΓ)
+end

test/runtests.jl

@@ -1,18 +1,15 @@
+using Gridap
 using HyperFEM
 using Test
-using Gridap
 
 
 @testset "HyperFEMTests" verbose = true begin
 
-    @time begin
-        include("../test/TestConstitutiveModels/runtests.jl")
-    end
-    @time begin
-        include("./TestTensorAlgebra/runtests.jl")
-    end
-    @time begin
-        include("../test/TestWeakForms/runtests.jl")
-    end
+    include("TestConstitutiveModels/runtests.jl")
+
+    include("TestTensorAlgebra/runtests.jl")
+
+    include("TestWeakForms/runtests.jl")
 
+    include("SimulationsTests/runtests.jl")
 end;