diff --git a/src/PhysicalModels/ElectroMechanicalModels.jl b/src/PhysicalModels/ElectroMechanicalModels.jl index 45344e7..cf8d21b 100644 --- a/src/PhysicalModels/ElectroMechanicalModels.jl +++ b/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 diff --git a/src/WeakForms/ElectroMechanics.jl b/src/WeakForms/ElectroMechanics.jl index c4fbe86..4416e1a 100644 --- a/src/WeakForms/ElectroMechanics.jl +++ b/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 # ------------------- diff --git a/src/WeakForms/WeakForms.jl b/src/WeakForms/WeakForms.jl index 9bbcbeb..151be58 100644 --- a/src/WeakForms/WeakForms.jl +++ b/src/WeakForms/WeakForms.jl @@ -36,10 +36,10 @@ 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 """ @@ -47,22 +47,10 @@ end 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 diff --git a/test/TestConstitutiveModels/ElectroMechanicalTests.jl b/test/TestConstitutiveModels/ElectroMechanicalTests.jl index 8dcf6be..34c168f 100644 --- a/test/TestConstitutiveModels/ElectroMechanicalTests.jl +++ b/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 diff --git a/test/TestWeakForms/WeakForms.jl b/test/TestWeakForms/WeakForms.jl index d000845..79c16ca 100644 --- a/test/TestWeakForms/WeakForms.jl +++ b/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 diff --git a/test/data/ViscoElasticSimulation.jl b/test/data/ViscoElasticSimulation.jl index 4c6aba2..1385776 100644 --- a/test/data/ViscoElasticSimulation.jl +++ b/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,10 +60,11 @@ 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) @@ -70,7 +72,6 @@ function visco_elastic_simulation(;t_end=15, writevtk=true, verbose=true) λx = Float64[] σΓ = Float64[] - F,_,_ = get_Kinematics(k) function driverpost(post) σh11, _... = Cauchy(cons_model, Kinematics(Mechano,Solid),uh, unh, state_vars, Ω, dΩ, 0.0, Δt)