[Ready to merge] Electro-mechanic coupling for visco-elastic constitutive models

benchmark/ConstitutiveModelsBenchmark/ViscousModelsBenchmarks.jl

@@ -14,7 +14,7 @@ function benchmark_viscous_model()
   J = det(F)
   Uvn *= J^(-1/3)
   λvn = 1e-3
-  Avn = VectorValue(Uvn.data..., λvn)
+  Avn = VectorValue(Uvn..., λvn)
   SUITE["Constitutive models"]["Visco-elastic Ψ"] = @benchmarkable $Ψ($F, $Fn, $Avn)
   SUITE["Constitutive models"]["Visco-elastic ∂Ψu"] = @benchmarkable $∂Ψu($F, $Fn, $Avn)
   SUITE["Constitutive models"]["Visco-elastic ∂Ψuu"] = @benchmarkable $∂Ψuu($F, $Fn, $Avn)

benchmark/Simulations/ViscoElasticBenchmark.jl

@@ -56,7 +56,7 @@ function viscousbenchmark()
 
   # residual and jacobian
   uh = FEFunction(Uu, zero_free_values(Uu))
-  unh = FEFunction(Uu, zero_free_values(Uu))
+  unh = FEFunction(Uun, zero_free_values(Uun))
   state_vars = initializeStateVariables(cons_model, dΩ)
 
   res(Λ) = (u,v)->residual(cons_model, u, v, dΩ, t_end * Λ, Δt, unh, state_vars)
@@ -67,7 +67,7 @@ function viscousbenchmark()
   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!(cons_model, Δt, uh, unh, A)
+    updateStateVariables!(A, cons_model, Δt, uh, unh)
   end
   post_model = PostProcessor(comp_model, driverpost; is_vtk=false, filepath="")
 

src/PhysicalModels/ElectroMechanicalModels.jl

@@ -0,0 +1,81 @@
+
+struct ElectroMechModel{M<:Mechano, E<:Electro} <: ElectroMechano
+  mechano::M
+  electro::E
+
+  function ElectroMechModel(mechano::M, electro::E) where {M<:Mechano, E<:Electro}
+    new{M,E}(mechano,electro)
+  end
+
+  function ElectroMechModel(;mechano::M, electro::E) where {M<:Mechano, E<:Electro}
+    new{M,E}(mechano,electro)
+  end
+
+  function (obj::ElectroMechModel)(Λ::Float64=1.0)
+    Ψm, ∂Ψm_u, ∂Ψm_uu = obj.mechano(Λ)
+    Ψem, ∂Ψem_u, ∂Ψem_φ, ∂Ψem_uu, ∂Ψem_φu, ∂Ψem_φφ = _getCoupling(obj.mechano, obj.electro, Λ)
+    Ψ(F, E) = Ψm(F) + Ψem(F, E)
+    ∂Ψu(F, E) = ∂Ψm_u(F) + ∂Ψem_u(F, E)
+    ∂Ψφ(F, E) = ∂Ψem_φ(F, E)
+    ∂Ψuu(F, E) = ∂Ψm_uu(F) + ∂Ψem_uu(F, E)
+    ∂Ψφu(F, E) = ∂Ψem_φu(F, E)
+    ∂Ψφφ(F, E) = ∂Ψem_φφ(F, E)
+    return (Ψ, ∂Ψu, ∂Ψφ, ∂Ψuu, ∂Ψφu, ∂Ψφφ)
+  end
+
+  function (obj::ElectroMechModel{<:ViscoElastic,<:Electro})(Λ::Float64=1.0; Δt)
+    Ψm, ∂Ψm_u, ∂Ψm_uu = obj.mechano(Λ, Δt=Δt)
+    Ψem, ∂Ψem_u, ∂Ψem_φ, ∂Ψem_uu, ∂Ψem_φu, ∂Ψem_φφ = _getCoupling(obj.mechano, obj.electro, Λ)
+    Ψ(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)
+    return (Ψ, ∂Ψu, ∂Ψφ, ∂Ψuu, ∂Ψφu, ∂Ψφφ)
+  end
+end
+
+ViscoElectricModel = ElectroMechModel{<:ViscoElastic,<:Electro}
+
+function initializeStateVariables(obj::ElectroMechano, points::Measure)
+  initializeStateVariables(obj.mechano, points)
+end
+
+function updateStateVariables!(state, obj::ElectroMechano, args...)
+  updateStateVariables!(state, obj.mechano, args)
+end
+
+
+function _getCoupling(mec::Mechano, elec::Electro, Λ::Float64)
+  _, H, J = get_Kinematics(mec.Kinematic; Λ=Λ)
+
+  # Energy #
+  HE(F, E) = H(F) * E
+  HEHE(F, E) = HE(F, E) ⋅ HE(F, E)
+  Ψem(F, E) = (-elec.ε / (2.0 * J(F))) * HEHE(F, E)
+  # First Derivatives #
+  ∂Ψem_∂H(F, E) = (-elec.ε / (J(F))) * (HE(F, E) ⊗ E)
+  ∂Ψem_∂J(F, E) = (+elec.ε / (2.0 * J(F)^2.0)) * HEHE(F, E)
+  ∂Ψem_∂E(F, E) = (-elec.ε / (J(F))) * (H(F)' * HE(F, E))
+  ∂Ψem_u(F, E) = ∂Ψem_∂H(F, E) × F + ∂Ψem_∂J(F, E) * H(F)
+  ∂Ψem_φ(F, E) = ∂Ψem_∂E(F, E)
+
+  # Second Derivatives #
+  ∂Ψem_HH(F, E) = (-elec.ε / (J(F))) * (I3 ⊗₁₃²⁴ (E ⊗ E))
+  ∂Ψem_HJ(F, E) = (+elec.ε / (J(F))^2.0) * (HE(F, E) ⊗ E)
+  ∂Ψem_JJ(F, E) = (-elec.ε / (J(F))^3.0) * HEHE(F, E)
+  ∂Ψem_uu(F, E) = (F × (∂Ψem_HH(F, E) × F)) +
+                  H(F) ⊗₁₂³⁴ (∂Ψem_HJ(F, E) × F) +
+                  (∂Ψem_HJ(F, E) × F) ⊗₁₂³⁴ H(F) +
+                  ∂Ψem_JJ(F, E) * (H(F) ⊗₁₂³⁴ H(F)) +
+                  ×ᵢ⁴(∂Ψem_∂H(F, E) + ∂Ψem_∂J(F, E) * F)
+
+  ∂Ψem_EH(F, E) = (-elec.ε / (J(F))) * ((I3 ⊗₁₃² HE(F, E)) + (H(F)' ⊗₁₂³ E))
+  ∂Ψem_EJ(F, E) = (+elec.ε / (J(F))^2.0) * (H(F)' * HE(F, E))
+  ∂Ψem_φu(F, E) = (∂Ψem_EH(F, E) × F) + (∂Ψem_EJ(F, E) ⊗₁²³ H(F))
+
+  ∂Ψem_φφ(F, E) = (-elec.ε / (J(F))) * (H(F)' * H(F))
+
+  return (Ψem, ∂Ψem_u, ∂Ψem_φ, ∂Ψem_uu, ∂Ψem_φu, ∂Ψem_φφ)
+end

src/PhysicalModels/PhysicalModels.jl

@@ -97,8 +97,13 @@ abstract type FlexoElectro <: MultiPhysicalModel end
 abstract type MagnetoMechano <: MultiPhysicalModel end
 
 include("KinematicModels.jl")
+
 include("ViscousModels.jl")
+
+include("ElectroMechanicalModels.jl")
+
 include("PINNs.jl")
+
 export Kinematics
 export KinematicModel
 export EvolutiveKinematics
@@ -117,7 +122,7 @@ function initializeStateVariables(::PhysicalModel, points::Measure)
   return nothing
 end
 
-function updateStateVariables!(::PhysicalModel, vars...)
+function updateStateVariables!(::Any, ::PhysicalModel, vars...)
 end
 
 # ============================================
@@ -1319,8 +1324,8 @@ end
 struct FlexoElectroModel{A} <: FlexoElectro
   ElectroMechano::A
   κ::Float64
-  function FlexoElectroModel(; Mechano::Mechano, Electro::Electro, κ=1.0)
-    physmodel = ElectroMechModel(Mechano=Mechano, Electro=Electro)
+  function FlexoElectroModel(; mechano::Mechano, electro::Electro, κ=1.0)
+    physmodel = ElectroMechModel(mechano=mechano, electro=electro)
     A = typeof(physmodel)
     new{A}(physmodel, κ)
   end
@@ -1340,26 +1345,6 @@ struct FlexoElectroModel{A} <: FlexoElectro
 
 end
 
-struct ElectroMechModel{A,B} <: ElectroMechano
-  Mechano::A
-  Electro::B
-  function ElectroMechModel(; Mechano::Mechano, Electro::Electro)
-    A, B = typeof(Mechano), typeof(Electro)
-    new{A,B}(Mechano, Electro)
-  end
-  function (obj::ElectroMechModel)(Λ::Float64=1.0)
-    Ψm, ∂Ψm_u, ∂Ψm_uu = obj.Mechano(Λ)
-    Ψem, ∂Ψem_u, ∂Ψem_φ, ∂Ψem_uu, ∂Ψem_φu, ∂Ψem_φφ = _getCoupling(obj.Mechano, obj.Electro, Λ)
-    Ψ(F, E) = Ψm(F) + Ψem(F, E)
-    ∂Ψu(F, E) = ∂Ψm_u(F) + ∂Ψem_u(F, E)
-    ∂Ψφ(F, E) = ∂Ψem_φ(F, E)
-    ∂Ψuu(F, E) = ∂Ψm_uu(F) + ∂Ψem_uu(F, E)
-    ∂Ψφu(F, E) = ∂Ψem_φu(F, E)
-    ∂Ψφφ(F, E) = ∂Ψem_φφ(F, E)
-    return (Ψ, ∂Ψu, ∂Ψφ, ∂Ψuu, ∂Ψφu, ∂Ψφφ)
-  end
-
-end
 
 struct ThermoMechModel{A,B,C,D} <: ThermoMechano
   Thermo::A
@@ -1625,43 +1610,6 @@ end
 # Coupling terms for multiphysic
 # ===============================
 
-function _getCoupling(mec::Mechano, elec::IdealDielectric, Λ::Float64)
-  _, H, J = get_Kinematics(mec.Kinematic; Λ=Λ)
-
-  # Energy #
-  HE(F, E) = H(F) * E
-  HEHE(F, E) = HE(F, E) ⋅ HE(F, E)
-  Ψem(F, E) = (-elec.ε / (2.0 * J(F))) * HEHE(F, E)
-  # First Derivatives #
-  ∂Ψem_∂H(F, E) = (-elec.ε / (J(F))) * (HE(F, E) ⊗ E)
-  ∂Ψem_∂J(F, E) = (+elec.ε / (2.0 * J(F)^2.0)) * HEHE(F, E)
-  ∂Ψem_∂E(F, E) = (-elec.ε / (J(F))) * (H(F)' * HE(F, E))
-  ∂Ψem_u(F, E) = ∂Ψem_∂H(F, E) × F + ∂Ψem_∂J(F, E) * H(F)
-  # ∂Ψem_φ(F, E) = -∂Ψem_∂E(F, E)
-  ∂Ψem_φ(F, E) = ∂Ψem_∂E(F, E)
-
-  # Second Derivatives #
-  ∂Ψem_HH(F, E) = (-elec.ε / (J(F))) * (I3 ⊗₁₃²⁴ (E ⊗ E))
-  ∂Ψem_HJ(F, E) = (+elec.ε / (J(F))^2.0) * (HE(F, E) ⊗ E)
-  ∂Ψem_JJ(F, E) = (-elec.ε / (J(F))^3.0) * HEHE(F, E)
-  ∂Ψem_uu(F, E) = (F × (∂Ψem_HH(F, E) × F)) +
-                  H(F) ⊗₁₂³⁴ (∂Ψem_HJ(F, E) × F) +
-                  (∂Ψem_HJ(F, E) × F) ⊗₁₂³⁴ H(F) +
-                  ∂Ψem_JJ(F, E) * (H(F) ⊗₁₂³⁴ H(F)) +
-                  ×ᵢ⁴(∂Ψem_∂H(F, E) + ∂Ψem_∂J(F, E) * F)
-
-  ∂Ψem_EH(F, E) = (-elec.ε / (J(F))) * ((I3 ⊗₁₃² HE(F, E)) + (H(F)' ⊗₁₂³ E))
-  ∂Ψem_EJ(F, E) = (+elec.ε / (J(F))^2.0) * (H(F)' * HE(F, E))
-
-  # ∂Ψem_φu(F, E) = -(∂Ψem_EH(F, E) × F) - (∂Ψem_EJ(F, E) ⊗₁²³ H(F))
-  ∂Ψem_φu(F, E) = (∂Ψem_EH(F, E) × F) + (∂Ψem_EJ(F, E) ⊗₁²³ H(F))
-
-  ∂Ψem_φφ(F, E) = (-elec.ε / (J(F))) * (H(F)' * H(F))
-
-  return (Ψem, ∂Ψem_u, ∂Ψem_φ, ∂Ψem_uu, ∂Ψem_φu, ∂Ψem_φφ)
-
-end
-
 function _getCoupling(mec::Mechano, term::Thermo, Λ::Float64)
   _, H, J = get_Kinematics(mec.Kinematic; Λ=Λ)
 

src/PhysicalModels/ViscousModels.jl

@@ -10,7 +10,7 @@ struct ViscousIncompressible{T} <: Visco
   ShortTerm::Elasto
   τ::Float64
   Kinematic::T
-  function ViscousIncompressible(shortTerm, τ::Float64; kinematic::KinematicModel=Kinematics(Visco))
+  function ViscousIncompressible(shortTerm, τ::Float64; kinematic::KinematicModel=Kinematics(Visco))  # TODO: Make τ keyword argument
     new{typeof(kinematic)}(shortTerm, τ, kinematic)
   end
   function (obj::ViscousIncompressible)(Λ::Float64=1.0; Δt::Float64)
@@ -27,7 +27,7 @@ function initializeStateVariables(::ViscousIncompressible, points::Measure)
   CellState(v, points)
 end
 
-function updateStateVariables!(obj::ViscousIncompressible, Δt, u, un, stateVar)
+function updateStateVariables!(stateVar, obj::ViscousIncompressible, Δt, u, un)
   F, _, _ = get_Kinematics(obj.Kinematic)
   _, Se, ∂Se∂Ce = obj.ShortTerm(KinematicDescription{:SecondPiola}())
   return_mapping(A, F, Fn) = ReturnMapping(obj, Δt, Se, ∂Se∂Ce, F, Fn, A)
@@ -56,10 +56,10 @@ function initializeStateVariables(model::GeneralizedMaxwell, points::Measure)
   map(b -> initializeStateVariables(b, points), model.Branches)
 end
 
-function updateStateVariables!(model::GeneralizedMaxwell, Δt, u, un, stateVars)
+function updateStateVariables!(stateVars, model::GeneralizedMaxwell, Δt, u, un)
   @assert length(model.Branches) == length(stateVars)
   for (branch, state) in zip(model.Branches, stateVars)
-    updateStateVariables!(branch, Δt, u, un, state)
+    updateStateVariables!(state, branch, Δt, u, un)
   end
 end
 

src/WeakForms/ElectroMechanics.jl

@@ -0,0 +1,138 @@
+
+# =================
+# Hyper-elasticity
+# =================
+
+# -----------------
+# Stagered residual
+# -----------------
+
+function residual(physicalmodel::ElectroMechano, ::Type{Mechano}, (u, φ), v, dΩ, Λ=1.0)
+    DΨ    = physicalmodel(Λ)
+    F,_,_ = get_Kinematics(physicalmodel.mechano.Kinematic; Λ=Λ)
+    E     = get_Kinematics(physicalmodel.electro.Kinematic; Λ=Λ)
+    ∂Ψu   = DΨ[2]
+    ∫(∇(v)' ⊙ (∂Ψu ∘ (F∘∇(u)', E∘∇(φ))))dΩ
+end
+
+function residual(physicalmodel::ElectroMechano, ::Type{Electro}, (u, φ), vφ, dΩ, Λ=1.0)
+    DΨ    = physicalmodel(Λ)
+    F,_,_ = get_Kinematics(physicalmodel.mechano.Kinematic; Λ=Λ)
+    E     = get_Kinematics(physicalmodel.electro.Kinematic; Λ=Λ)
+    ∂Ψφ   = DΨ[3]
+    -1.0*∫(∇(vφ) ⋅ (∂Ψφ ∘ (F∘∇(u)', E∘∇(φ))))dΩ
+end
+
+# -----------------
+# Stagered jacobian
+# -----------------
+
+function jacobian(physicalmodel::ElectroMechano, ::Type{Mechano}, (u, φ), du, v, dΩ, Λ=1.0)
+    DΨ    = physicalmodel(Λ)
+    F,_,_ = get_Kinematics(physicalmodel.mechano.Kinematic; Λ=Λ)
+    E     = get_Kinematics(physicalmodel.electro.Kinematic; Λ=Λ)
+    ∂Ψuu  = DΨ[4]
+    ∫(∇(v)' ⊙ ((∂Ψuu ∘ (F∘∇(u)', E∘∇(φ))) ⊙ ∇(du)'))dΩ
+end
+
+function jacobian(physicalmodel::ElectroMechano, ::Type{Electro}, (u, φ), dφ, vφ, dΩ, Λ=1.0)
+    DΨ    = physicalmodel(Λ)
+    F,_,_ = get_Kinematics(physicalmodel.mechano.Kinematic; Λ=Λ)
+    E     = get_Kinematics(physicalmodel.electro.Kinematic; Λ=Λ)
+    ∂Ψφφ  = DΨ[6]
+    ∫(∇(vφ)' ⋅ ((∂Ψφφ ∘ (F∘∇(u)', E∘∇(φ))) ⋅ ∇(dφ)))dΩ
+end
+
+function jacobian(physicalmodel::ElectroMechano, ::Type{ElectroMechano}, (u, φ), (du, dφ), (v, vφ), dΩ, Λ=1.0)
+    DΨ    = physicalmodel(Λ)
+    F,_,_ = get_Kinematics(physicalmodel.mechano.Kinematic; Λ=Λ)
+    E     = get_Kinematics(physicalmodel.electro.Kinematic; Λ=Λ)
+    ∂Ψφu  = DΨ[5]
+    -1.0*∫(∇(dφ) ⋅ ((∂Ψφu ∘ (F∘∇(u)', E∘∇(φ))) ⊙ ∇(v)'))dΩ -
+        ∫(∇(vφ) ⋅ ((∂Ψφu ∘ (F∘∇(u)', E∘∇(φ))) ⊙ ∇(du)'))dΩ 
+end
+
+# -------------------
+# Monolithic strategy
+# -------------------
+
+function residual(physicalmodel::ElectroMechano, (u, φ), (v, vφ), dΩ, Λ=1.0)
+    residual(physicalmodel, Mechano, (u, φ), v, dΩ, Λ) +
+    residual(physicalmodel, Electro, (u, φ), vφ, dΩ, Λ)
+end
+
+
+function jacobian(physicalmodel::ElectroMechano, (u, φ), (du, dφ), (v, vφ), dΩ, Λ=1.0)
+    jacobian(physicalmodel, Mechano, (u, φ), du, v, dΩ, Λ)+
+    jacobian(physicalmodel, Electro, (u, φ), dφ, vφ, dΩ, Λ)+
+    jacobian(physicalmodel, ElectroMechano, (u, φ), (du, dφ), (v, vφ), dΩ, Λ)
+end
+
+
+# =================
+# Visco-elasticity
+# =================
+
+# -----------------
+# Stagered residual
+# -----------------
+
+function residual(physicalmodel::ViscoElectricModel, ::Type{Mechano}, (u, φ), v, dΩ, Λ, Δt, un, A)
+    DΨ    = physicalmodel(Λ, Δt=Δt)
+    F,_,_ = get_Kinematics(physicalmodel.mechano.Kinematic; Λ=Λ)
+    E     = get_Kinematics(physicalmodel.electro.Kinematic; Λ=Λ)
+    ∂Ψu   = DΨ[2]
+    ∫(∇(v)' ⊙ (∂Ψu ∘ (F∘∇(u)', F∘∇(un)', E∘∇(φ), A...)))dΩ
+end
+
+function residual(physicalmodel::ViscoElectricModel, ::Type{Electro}, (u, φ), vφ, dΩ, Λ, Δt, un, A)
+    DΨ    = physicalmodel(Λ, Δt=Δt)
+    F,_,_ = get_Kinematics(physicalmodel.mechano.Kinematic; Λ=Λ)
+    E     = get_Kinematics(physicalmodel.electro.Kinematic; Λ=Λ)
+    ∂Ψφ   = DΨ[3]
+    -1.0*∫(∇(vφ) ⋅ (∂Ψφ ∘ (F∘∇(u)', F∘∇(un)', E∘∇(φ), A...)))dΩ
+end
+
+# -----------------
+# Stagered jacobian
+# -----------------
+
+function jacobian(physicalmodel::ViscoElectricModel, ::Type{Mechano}, (u, φ), du, v, dΩ, Λ, Δt, un, A)
+    DΨ    = physicalmodel(Λ, Δt=Δt)
+    F,_,_ = get_Kinematics(physicalmodel.mechano.Kinematic; Λ=Λ)
+    E     = get_Kinematics(physicalmodel.electro.Kinematic; Λ=Λ)
+    ∂Ψuu  = DΨ[4]
+    ∫(∇(v)' ⊙ ((∂Ψuu ∘ (F∘∇(u)', F∘∇(un)', E∘∇(φ), A...)) ⊙ (∇(du)')))dΩ
+end
+
+function jacobian(physicalmodel::ViscoElectricModel, ::Type{Electro}, (u, φ), dφ, vφ, dΩ, Λ, Δt, un, A)
+    DΨ    = physicalmodel(Λ, Δt=Δt)
+    F,_,_ = get_Kinematics(physicalmodel.mechano.Kinematic; Λ=Λ)
+    E     = get_Kinematics(physicalmodel.electro.Kinematic; Λ=Λ)
+    ∂Ψφφ  = DΨ[6]
+    ∫(∇(vφ)' ⋅ ((∂Ψφφ ∘ (F∘∇(u)', F∘∇(un)', E∘∇(φ), A...)) ⋅ ∇(dφ)))dΩ
+end
+
+function jacobian(physicalmodel::ViscoElectricModel, ::Type{ElectroMechano}, (u, φ), (du, dφ), (v, vφ), dΩ, Λ, Δt, un, A)
+    DΨ    = physicalmodel(Λ, Δt=Δt)
+    F,_,_ = get_Kinematics(physicalmodel.mechano.Kinematic; Λ=Λ)
+    E     = get_Kinematics(physicalmodel.electro.Kinematic; Λ=Λ)
+    ∂Ψφ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
+# -------------------
+
+function residual(physicalmodel::ViscoElectricModel, (u, φ), (v, vφ), dΩ, Λ, Δt, un, A)
+    residual(physicalmodel, Mechano, (u, φ), v, dΩ, Λ, Δt, un, A) +
+    residual(physicalmodel, Electro, (u, φ), vφ, dΩ, Λ, Δt, un, A)
+end
+
+function jacobian(physicalmodel::ViscoElectricModel, (u, φ), (du, dφ), (v, vφ), dΩ, Λ, Δt, un, A)
+    jacobian(physicalmodel, Mechano, (u, φ), du, v, dΩ, Λ, Δt, un, A) +
+    jacobian(physicalmodel, Electro, (u, φ), dφ, vφ, dΩ, Λ, Δt, un, A) +
+    jacobian(physicalmodel, ElectroMechano, (u, φ), (du, dφ), (v, vφ), dΩ, Λ, Δt, un, A)
+end