Merge pull request #57 from MultiSimOLab/visco-dispatch

Author: jmartfrut

src/PhysicalModels/ElectroMechanicalModels.jl

@@ -36,12 +36,12 @@ struct ElectroMechModel{E<:Electro,M<:Mechano} <: ElectroMechano
   function (obj::ElectroMechModel{<:Electro,<:ViscoElastic})(Λ::Float64=1.0; Δt)
     Ψm, ∂Ψm_u, ∂Ψm_uu = obj.mechano(Λ, Δt=Δt)
     Ψem, ∂Ψem_u, ∂Ψem_φ, ∂Ψem_uu, ∂Ψem_φu, ∂Ψem_φφ = _getCoupling(obj.electro, obj.mechano, Λ)
-    Ψ(F, Fn, E, A...) = Ψm(F, Fn, A...) + Ψem(F, E)
-    ∂Ψu(F, Fn, E, A...) = ∂Ψm_u(F, Fn, A...) + ∂Ψem_u(F, E)
-    ∂Ψφ(F, Fn, E, A...) = ∂Ψem_φ(F, E)
-    ∂Ψuu(F, Fn, E, A...) = ∂Ψm_uu(F, Fn, A...) + ∂Ψem_uu(F, E)
-    ∂Ψφu(F, Fn, E, A...) = ∂Ψem_φu(F, E)
-    ∂Ψφφ(F, Fn, E, A...) = ∂Ψem_φφ(F, E)
+    Ψ(F, E, Fn, A...) = Ψm(F, Fn, A...) + Ψem(F, E)
+    ∂Ψu(F, E, Fn, A...) = ∂Ψm_u(F, Fn, A...) + ∂Ψem_u(F, E)
+    ∂Ψφ(F, E, Fn, A...) = ∂Ψem_φ(F, E)
+    ∂Ψuu(F, E, Fn, A...) = ∂Ψm_uu(F, Fn, A...) + ∂Ψem_uu(F, E)
+    ∂Ψφu(F, E, Fn, A...) = ∂Ψem_φu(F, E)
+    ∂Ψφφ(F, E, Fn, A...) = ∂Ψem_φφ(F, E)
     return (Ψ, ∂Ψu, ∂Ψφ, ∂Ψuu, ∂Ψφu, ∂Ψφφ)
   end
 

src/WeakForms/ElectroMechanics.jl

@@ -7,49 +7,49 @@
 # Stagered residual
 # -----------------
 
-function residual(physicalmodel::ElectroMechano, ::Type{Mechano},  kine::NTuple{2,KinematicModel}, (u, φ), v, dΩ, Λ=1.0)
-    DΨ    = physicalmodel(Λ)
+function residual(physicalmodel::ElectroMechano, ::Type{Mechano}, kine::NTuple{2,KinematicModel}, (u, φ), v, dΩ, Λ=1.0, vars...; kwargs...)
+    DΨ    = physicalmodel(; kwargs...)
     F,_,_ = get_Kinematics(kine[1]; Λ=Λ)
     E     = get_Kinematics(kine[2]; Λ=Λ)
     ∂Ψu   = DΨ[2]
-    ∫(∇(v)' ⊙ (∂Ψu ∘ (F∘∇(u)', E∘∇(φ))))dΩ
+    ∫(∇(v)' ⊙ (∂Ψu ∘ (F∘∇(u)', E∘∇(φ), vars...)))dΩ
 end
 
-function residual(physicalmodel::ElectroMechano, ::Type{Electro},  kine::NTuple{2,KinematicModel}, (u, φ), vφ, dΩ, Λ=1.0)
-    DΨ    = physicalmodel(Λ)
+function residual(physicalmodel::ElectroMechano, ::Type{Electro}, kine::NTuple{2,KinematicModel}, (u, φ), vφ, dΩ, Λ=1.0, vars...; kwargs...)
+    DΨ    = physicalmodel(; kwargs...)
     F,_,_ = get_Kinematics(kine[1]; Λ=Λ)
     E     = get_Kinematics(kine[2]; Λ=Λ)
     ∂Ψφ   = DΨ[3]
-    -1.0*∫(∇(vφ) ⋅ (∂Ψφ ∘ (F∘∇(u)', E∘∇(φ))))dΩ
+    -1.0*∫(∇(vφ) ⋅ (∂Ψφ ∘ (F∘∇(u)', E∘∇(φ), vars...)))dΩ
 end
 
 # -----------------
 # Stagered jacobian
 # -----------------
 
-function jacobian(physicalmodel::ElectroMechano, ::Type{Mechano},  kine::NTuple{2,KinematicModel}, (u, φ), du, v, dΩ, Λ=1.0)
-    DΨ    = physicalmodel(Λ)
+function jacobian(physicalmodel::ElectroMechano, ::Type{Mechano}, kine::NTuple{2,KinematicModel}, (u, φ), du, v, dΩ, Λ=1.0, vars...; kwargs...)
+    DΨ    = physicalmodel(; kwargs...)
     F,_,_ = get_Kinematics(kine[1]; Λ=Λ)
     E     = get_Kinematics(kine[2]; Λ=Λ)
     ∂Ψuu  = DΨ[4]
-    ∫(∇(v)' ⊙ ((∂Ψuu ∘ (F∘∇(u)', E∘∇(φ))) ⊙ ∇(du)'))dΩ
+    ∫(∇(v)' ⊙ ((∂Ψuu ∘ (F∘∇(u)', E∘∇(φ), vars...)) ⊙ ∇(du)'))dΩ
 end
 
-function jacobian(physicalmodel::ElectroMechano, ::Type{Electro},  kine::NTuple{2,KinematicModel}, (u, φ), dφ, vφ, dΩ, Λ=1.0)
-    DΨ    = physicalmodel(Λ)
+function jacobian(physicalmodel::ElectroMechano, ::Type{Electro}, kine::NTuple{2,KinematicModel}, (u, φ), dφ, vφ, dΩ, Λ=1.0, vars...; kwargs...)
+    DΨ    = physicalmodel(; kwargs...)
     F,_,_ = get_Kinematics(kine[1]; Λ=Λ)
     E     = get_Kinematics(kine[2]; Λ=Λ)
     ∂Ψφφ  = DΨ[6]
-    ∫(∇(vφ)' ⋅ ((∂Ψφφ ∘ (F∘∇(u)', E∘∇(φ))) ⋅ ∇(dφ)))dΩ
+    ∫(∇(vφ)' ⋅ ((∂Ψφφ ∘ (F∘∇(u)', E∘∇(φ), vars...)) ⋅ ∇(dφ)))dΩ
 end
 
-function jacobian(physicalmodel::ElectroMechano, ::Type{ElectroMechano},  kine::NTuple{2,KinematicModel},(u, φ), (du, dφ), (v, vφ), dΩ, Λ=1.0)
-    DΨ    = physicalmodel(Λ)
+function jacobian(physicalmodel::ElectroMechano, ::Type{ElectroMechano}, kine::NTuple{2,KinematicModel}, (u, φ), (du, dφ), (v, vφ), dΩ, Λ=1.0, vars...; kwargs...)
+    DΨ    = physicalmodel(; kwargs...)
     F,_,_ = get_Kinematics(kine[1]; Λ=Λ)
     E     = get_Kinematics(kine[2]; Λ=Λ)
     ∂Ψφu  = DΨ[5]
-    -1.0*∫(∇(dφ) ⋅ ((∂Ψφu ∘ (F∘∇(u)', E∘∇(φ))) ⊙ ∇(v)'))dΩ -
-        ∫(∇(vφ) ⋅ ((∂Ψφu ∘ (F∘∇(u)', E∘∇(φ))) ⊙ ∇(du)'))dΩ 
+    -1.0*∫(∇(dφ) ⋅ ((∂Ψφu ∘ (F∘∇(u)', E∘∇(φ), vars...)) ⊙ ∇(v)'))dΩ -
+         ∫(∇(vφ) ⋅ ((∂Ψφu ∘ (F∘∇(u)', E∘∇(φ), vars...)) ⊙ ∇(du)'))dΩ 
 end
 
 # -------------------
@@ -68,60 +68,6 @@ function jacobian(physicalmodel::ElectroMechano,  kine::NTuple{2,KinematicModel}
     jacobian(physicalmodel, ElectroMechano, kine, (u, φ), (du, dφ), (v, vφ), dΩ, Λ)
 end
 
-
-# =================
-# Visco-elasticity
-# =================
-
-# -----------------
-# Stagered residual
-# -----------------
-
-function residual(physicalmodel::ViscoElectricModel, ::Type{Mechano},  kine::NTuple{2,KinematicModel}, (u, φ), v, dΩ, Λ, Δt, un, A)
-    DΨ    = physicalmodel(Λ, Δt=Δt)
-    F,_,_ = get_Kinematics(kine[1]; Λ=Λ)
-    E     = get_Kinematics(kine[2]; Λ=Λ)
-    ∂Ψu   = DΨ[2]
-    ∫(∇(v)' ⊙ (∂Ψu ∘ (F∘∇(u)', F∘∇(un)', E∘∇(φ), A...)))dΩ
-end
-
-function residual(physicalmodel::ViscoElectricModel, ::Type{Electro}, kine::NTuple{2,KinematicModel},  (u, φ), vφ, dΩ, Λ, Δt, un, A)
-    DΨ    = physicalmodel(Λ, Δt=Δt)
-    F,_,_ = get_Kinematics(kine[1]; Λ=Λ)
-    E     = get_Kinematics(kine[2]; Λ=Λ)
-    ∂Ψφ   = DΨ[3]
-    -1.0*∫(∇(vφ) ⋅ (∂Ψφ ∘ (F∘∇(u)', F∘∇(un)', E∘∇(φ), A...)))dΩ
-end
-
-# -----------------
-# Stagered jacobian
-# -----------------
-
-function jacobian(physicalmodel::ViscoElectricModel, ::Type{Mechano},  kine::NTuple{2,KinematicModel}, (u, φ), du, v, dΩ, Λ, Δt, un, A)
-    DΨ    = physicalmodel(Λ, Δt=Δt)
-    F,_,_ = get_Kinematics(kine[1]; Λ=Λ)
-    E     = get_Kinematics(kine[2]; Λ=Λ)
-    ∂Ψuu  = DΨ[4]
-    ∫(∇(v)' ⊙ ((∂Ψuu ∘ (F∘∇(u)', F∘∇(un)', E∘∇(φ), A...)) ⊙ (∇(du)')))dΩ
-end
-
-function jacobian(physicalmodel::ViscoElectricModel, ::Type{Electro},  kine::NTuple{2,KinematicModel}, (u, φ), dφ, vφ, dΩ, Λ, Δt, un, A)
-    DΨ    = physicalmodel(Λ, Δt=Δt)
-    F,_,_ = get_Kinematics(kine[1]; Λ=Λ)
-    E     = get_Kinematics(kine[2]; Λ=Λ)
-    ∂Ψφφ  = DΨ[6]
-    ∫(∇(vφ)' ⋅ ((∂Ψφφ ∘ (F∘∇(u)', F∘∇(un)', E∘∇(φ), A...)) ⋅ ∇(dφ)))dΩ
-end
-
-function jacobian(physicalmodel::ViscoElectricModel, ::Type{ElectroMechano},  kine::NTuple{2,KinematicModel}, (u, φ), (du, dφ), (v, vφ), dΩ, Λ, Δt, un, A)
-    DΨ    = physicalmodel(Λ, Δt=Δt)
-    F,_,_ = get_Kinematics(kine[1]; Λ=Λ)
-    E     = get_Kinematics(kine[2]; Λ=Λ)
-    ∂Ψφu  = DΨ[5]
-    -1.0*∫(∇(dφ) ⋅ ((∂Ψφu ∘ (F∘∇(u)', F∘∇(un)', E∘∇(φ), A...)) ⊙ (∇(v)')))dΩ -
-        ∫(∇(vφ) ⋅ ((∂Ψφu ∘ (F∘∇(u)', F∘∇(un)', E∘∇(φ), A...)) ⊙ (∇(du)')))dΩ 
-end
-
 # -------------------
 # Monolithic strategy
 # -------------------

src/WeakForms/WeakForms.jl

@@ -36,33 +36,21 @@ end
 
 Calculate the residual using the given constitutive model and finite element functions.
 """
-function residual(physicalmodel::Mechano, km::KinematicModel ,u, v, dΩ, Λ=1.0, Δt=0.0, vars...)
-  _, ∂Ψu, _ = physicalmodel(Λ)
+function residual(physicalmodel::Mechano, km::KinematicModel, u, v, dΩ, Λ=1.0, vars...; kwargs...)
+  _, ∂Ψu, _ = physicalmodel(; kwargs...)
   F, _, _   = get_Kinematics(km; Λ=Λ)
-  ∫(∇(v)' ⊙ (∂Ψu ∘ (F∘∇(u)')))dΩ
+  ∫(∇(v)' ⊙ (∂Ψu ∘ (F∘∇(u)', vars...)))dΩ
 end
 
 """
     jacobian(...)::Gridap.CellData.Integrand
 
 Calculate the jacobian using the given constitutive model and finite element functions.
 """
-function jacobian(physicalmodel::Mechano, km::KinematicModel, u, du, v, dΩ, Λ=1.0, Δt=0.0, vars...)
-  _, _, ∂Ψuu = physicalmodel(Λ)
+function jacobian(physicalmodel::Mechano, km::KinematicModel, u, du, v, dΩ, Λ=1.0, vars...; kwargs...)
+  _, _, ∂Ψuu = physicalmodel(; kwargs...)
   F, _, _    = get_Kinematics(km; Λ=Λ)
-  ∫(∇(v)' ⊙ ((∂Ψuu ∘ (F∘∇(u)')) ⊙ ∇(du)'))dΩ
-end
-
-function residual(physicalmodel::ViscoElastic, km::KinematicModel,u, v, dΩ, t, Δt, un, A)
-  _, ∂Ψu, _ = physicalmodel(t, Δt=Δt)
-  F, _, _   = get_Kinematics(km, Λ=t)
-  ∫(∇(v)' ⊙ (∂Ψu ∘ (F∘∇(u)', F∘∇(un)', A...)))dΩ
-end
-
-function jacobian(physicalmodel::ViscoElastic, km::KinematicModel, u, du, v, dΩ, t, Δt, un, A)
-  _, _, ∂Ψuu = physicalmodel(t, Δt=Δt)
-  F, _, _   = get_Kinematics(km, Λ=t)
-  ∫(∇(v)' ⊙ (inner ∘ (∂Ψuu∘(F∘∇(u)', F∘∇(un)', A...), ∇(du)')))dΩ
+  ∫(∇(v)' ⊙ ((∂Ψuu ∘ (F∘∇(u)', vars...)) ⊙ ∇(du)'))dΩ
 end
 
 

test/TestConstitutiveModels/ElectroMechanicalTests.jl

@@ -115,8 +115,8 @@ end
   λvn = 1e-3
   Avn = VectorValue(Uvn..., λvn)
   Ψ, ∂Ψu, ∂Ψφ, ∂Ψuu, ∂Ψφu, ∂Ψφφ = model(Δt=0.01)
-  @test norm(∂Ψu(F(∇u), F(∇un), E(∇φ), Avn)) ≈ 25.049301121178615
-  @test norm(∂Ψuu(F(∇u), F(∇un), E(∇φ), Avn)) ≈ 3110.7607787445168
+  @test norm(∂Ψu(F(∇u), E(∇φ), F(∇un), Avn)) ≈ 25.049301121178615
+  @test norm(∂Ψuu(F(∇u), E(∇φ), F(∇un), Avn)) ≈ 3110.7607787445168
 end
 
 
@@ -137,6 +137,6 @@ end
   λvn = 1e-3
   Avn = VectorValue(Uvn..., λvn)
   Ψ, ∂Ψu, ∂Ψφ, ∂Ψuu, ∂Ψφu, ∂Ψφφ = model(Δt=0.01)
-  @test norm(∂Ψu(F(∇u), F(∇un), E(∇φ), Avn, Avn)) ≈ 25.102080194257017
-  @test norm(∂Ψuu(F(∇u), F(∇un), E(∇φ), Avn, Avn)) ≈ 3110.9722775475557
+  @test norm(∂Ψu(F(∇u), E(∇φ), F(∇un), Avn, Avn)) ≈ 25.102080194257017
+  @test norm(∂Ψuu(F(∇u), E(∇φ), F(∇un), Avn, Avn)) ≈ 3110.9722775475557
 end

test/TestWeakForms/WeakForms.jl

@@ -99,22 +99,22 @@ end
 
 
     ke=Kinematics(Electro, Solid)
-    km=Kinematics(Mechano,Solid)
+    km=Kinematics(Mechano, Solid)
     k=(km,ke)
     function jach(uh, φh)
-        jac((du, dφ), (v, vφ)) = jacobian(modelelectro,k, (uh, φh), (du, dφ), (v, vφ), dΩ)
+        jac((du, dφ), (v, vφ)) = jacobian(modelelectro, k, (uh, φh), (du, dφ), (v, vφ), dΩ)
     end
     function jac_mech(uh, φh)
-        jac(du, v) = jacobian(modelelectro, Mechano, k,(uh, φh), du, v, dΩ)
+        jac(du, v) = jacobian(modelelectro, Mechano, k, (uh, φh), du, v, dΩ)
     end
     function jac_elech(uh, φh)
-        jac(dφ, vφ) = jacobian(modelelectro, Electro, k,(uh, φh), dφ, vφ, dΩ)
+        jac(dφ, vφ) = jacobian(modelelectro, Electro, k, (uh, φh), dφ, vφ, dΩ)
     end
 
     jac_ = assemble_matrix(jach(uh, φh), V, V)
     jac_m = assemble_matrix(jac_mech(uh, φh), Vu, Vu)
     jac_e = assemble_matrix(jac_elech(uh, φh), Vφ, Vφ)
-    ≈
+    
     @test norm(jac_) ≈ 18.934585248125135
     @test jac_[1] ≈ 0.7777777777777775
     @test jac_[end] ≈ -1.3333333333333326

test/data/ViscoElasticSimulation.jl

@@ -29,6 +29,7 @@ function visco_elastic_simulation(;t_end=15, writevtk=true, verbose=true)
   hyper_elastic_model = NeoHookean3D(λ=λ, μ=μ)
   viscous_branch = ViscousIncompressible(IncompressibleNeoHookean3D(λ=0., μ=μv₁), τ=τv₁)
   cons_model = GeneralizedMaxwell(hyper_elastic_model, viscous_branch)
+  k=Kinematics(Mechano,Solid)
 
   # Dirichlet boundary conditions
   strain = 0.5
@@ -59,18 +60,18 @@ function visco_elastic_simulation(;t_end=15, writevtk=true, verbose=true)
   uh = FEFunction(Uu, zero_free_values(Uu))
   unh = FEFunction(Uun, zero_free_values(Uun))
   state_vars = initializeStateVariables(cons_model, dΩ)
+  F,_,_ = get_Kinematics(k)
+  Fnh = F∘∇(unh)'
 
-  k=Kinematics(Mechano,Solid)
-  res(Λ) = (u,v)->residual(cons_model, k, u, v, dΩ, t_end * Λ, Δt, unh, state_vars)
-  jac(Λ) = (u,du,v)->jacobian(cons_model, k, u, du, v, dΩ, t_end * Λ, Δt, unh, state_vars)
+  res(Λ) = (u,v)->residual(cons_model, k, u, v, dΩ, t_end * Λ, Fnh, state_vars...; Δt=Δt)
+  jac(Λ) = (u,du,v)->jacobian(cons_model, k, u, du, v, dΩ, t_end * Λ, Fnh, state_vars...; Δt=Δt)
 
   ls = LUSolver()
   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)
 
   λx = Float64[]
   σΓ = Float64[]
-  F,_,_ = get_Kinematics(k)
 
   function driverpost(post)
     σh11, _... = Piola(cons_model, Kinematics(Mechano,Solid),uh, unh, state_vars, Ω, dΩ, 0.0, Δt)