Merge pull request #47 from MultiSimOLab/test

Adds HGO_1Fiber model and tests

Author: jmartfrut

README.md

@@ -39,7 +39,6 @@ Example of a Monolithic implementation for the simulation of Nonlinear ElectroMe
 
 ```julia
 using HyperFEM
-using HyperFEM: jacobian, solve!
 using Gridap, GridapGmsh, GridapSolvers, DrWatson
 using GridapSolvers.NonlinearSolvers
 using Gridap.FESpaces

src/Exports.jl

@@ -65,6 +65,7 @@ end
 @publish PhysicalModels ViscousIncompressible
 @publish PhysicalModels GeneralizedMaxwell
 @publish PhysicalModels HGO_4Fibers
+@publish PhysicalModels HGO_1Fiber
 
 @publish PhysicalModels Mechano
 @publish PhysicalModels Thermo

src/PhysicalModels/ElectroMechanicalModels.jl

@@ -1,17 +1,17 @@
 
-struct ElectroMechModel{E<:Electro, M<:Mechano} <: ElectroMechano
+struct ElectroMechModel{E<:Electro,M<:Mechano} <: ElectroMechano
   electro::E
   mechano::M
 
-  function ElectroMechModel(electro::E, mechano::M) where {E<:Electro, M<:Mechano}
+  function ElectroMechModel(electro::E, mechano::M) where {E<:Electro,M<:Mechano}
     new{E,M}(electro, mechano)
   end
 
-  function ElectroMechModel(; electro::E, mechano::M) where {E<:Electro, M<:Mechano}
+  function ElectroMechModel(; electro::E, mechano::M) where {E<:Electro,M<:Mechano}
     new{E,M}(electro, mechano)
   end
 
-  function (obj::ElectroMechModel)(Λ::Float64=1.0)
+  function (obj::ElectroMechModel{<:Electro,<:IsoElastic})(Λ::Float64=1.0)
     Ψm, ∂Ψm_u, ∂Ψm_uu = obj.mechano(Λ)
     Ψem, ∂Ψem_u, ∂Ψem_φ, ∂Ψem_uu, ∂Ψem_φu, ∂Ψem_φφ = _getCoupling(obj.electro, obj.mechano, Λ)
     Ψ(F, E) = Ψm(F) + Ψem(F, E)
@@ -22,7 +22,17 @@ struct ElectroMechModel{E<:Electro, M<:Mechano} <: ElectroMechano
     ∂Ψφφ(F, E) = ∂Ψem_φφ(F, E)
     return (Ψ, ∂Ψu, ∂Ψφ, ∂Ψuu, ∂Ψφu, ∂Ψφφ)
   end
-  
+  function (obj::ElectroMechModel{<:Electro,<:AnisoElastic})(Λ::Float64=1.0)
+    Ψm, ∂Ψm_u, ∂Ψm_uu = obj.mechano(Λ)
+    Ψem, ∂Ψem_u, ∂Ψem_φ, ∂Ψem_uu, ∂Ψem_φu, ∂Ψem_φφ = _getCoupling(obj.electro, obj.mechano, Λ)
+    Ψ(F, E, N) = Ψm(F, N) + Ψem(F, E)
+    ∂Ψu(F, E, N) = ∂Ψm_u(F, N) + ∂Ψem_u(F, E)
+    ∂Ψφ(F, E, N) = ∂Ψem_φ(F, E)
+    ∂Ψuu(F, E, N) = ∂Ψm_uu(F, N) + ∂Ψem_uu(F, E)
+    ∂Ψφu(F, E, N) = ∂Ψem_φu(F, E)
+    ∂Ψφφ(F, E, N) = ∂Ψem_φφ(F, E)
+    return (Ψ, ∂Ψu, ∂Ψφ, ∂Ψuu, ∂Ψφu, ∂Ψφφ)
+  end
   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, Λ)
@@ -34,6 +44,7 @@ struct ElectroMechModel{E<:Electro, M<:Mechano} <: ElectroMechano
     ∂Ψφφ(F, Fn, E, A...) = ∂Ψem_φφ(F, E)
     return (Ψ, ∂Ψu, ∂Ψφ, ∂Ψuu, ∂Ψφu, ∂Ψφφ)
   end
+
 end
 
 ViscoElectricModel = ElectroMechModel{<:Electro,<:ViscoElastic}
@@ -92,12 +103,12 @@ struct FlexoElectroModel{EM<:ElectroMechano} <: FlexoElectro
   electromechano::EM
   κ::Float64
 
-  function FlexoElectroModel(electro::E, mechano::M; κ=1.0) where {E <: Electro, M <: Mechano}
+  function FlexoElectroModel(electro::E, mechano::M; κ=1.0) where {E<:Electro,M<:Mechano}
     physmodel = ElectroMechModel(electro, mechano)
     new{ElectroMechModel{E,M}}(physmodel, κ)
   end
 
-  function FlexoElectroModel(; electro::E, mechano::M, κ=1.0) where {E <: Electro, M <: Mechano}
+  function FlexoElectroModel(; electro::E, mechano::M, κ=1.0) where {E<:Electro,M<:Mechano}
     physmodel = ElectroMechModel(electro, mechano)
     new{ElectroMechModel{E,M}}(physmodel, κ)
   end

src/PhysicalModels/MechanicalModels.jl

@@ -757,6 +757,36 @@ end
 
 
 
+struct HGO_1Fiber <: AnisoElastic
+  c1::Float64
+  c2::Float64
+  Kinematic::Kinematics{Mechano}
+  function HGO_1Fiber(; c1::Float64, c2::Float64, Kinematic::KinematicModel=Kinematics(Mechano))
+    new(c1, c2, Kinematic)
+  end
+
+  function (obj::HGO_1Fiber)(Λ::Float64=1.0; Threshold=0.01)
+    c1, c2 = obj.c1, obj.c2
+
+    Ψ(F, N1) = begin
+      M1 = N1 / norm(N1)
+      c1 / (4 * c2) * (exp(c2 * ((F * M1) ⋅ (F * M1) - 1.0)^2.0) - 1.0)
+    end
+
+    ∂Ψ∂F(F, N1) = begin
+      M1 = N1 / norm(N1)
+      c1 * exp(c2 * ((F * M1) ⋅ (F * M1) - 1.0)^2.0) * ((F * M1) ⋅ (F * M1) - 1.0) * ((F * M1) ⊗ M1) 
+    end
+
+    ∂Ψ2∂F∂F(F, N1) = begin
+      M1 = N1 / norm(N1)
+      c1 * exp(c2 * ((F * M1) ⋅ (F * M1) - 1.0)^2.0) * ((4 * c2 * (((F * M1) ⋅ (F * M1) - 1.0)^2.0) + 2.0) * (((F * M1) ⊗ M1) ⊗ ((F * M1) ⊗ M1)) + ((F * M1) ⋅ (F * M1) - 1.0) * (I3 ⊗₁₃²⁴ (M1 ⊗ M1))) 
+    end
+
+    return (Ψ, ∂Ψ∂F, ∂Ψ2∂F∂F)
+  end
+end
+
 
 struct IncompressibleNeoHookean3D{A} <: IsoElastic
   λ::Float64

src/PhysicalModels/PhysicalModels.jl

@@ -55,6 +55,7 @@ export MagnetoMechModel
 export GeneralizedMaxwell
 export ViscousIncompressible
 export HGO_4Fibers
+export HGO_1Fiber
 
 export PhysicalModel
 export Mechano

test/TestConstitutiveModels/ElectroMechanicalTests.jl

@@ -7,6 +7,38 @@ const ∇u = TensorValue(1.0:9.0...) * 1e-3
 const ∇un = TensorValue(1.0:9.0...) * 5e-4
 
 
+@testset "Electro+4*HGO_1Fiber" begin
+  c1  =  [0.6639232500447778, 0.5532987701062146, 0.9912576142028674, 0.4951942011240962]
+  c2  =  [0.800583033264982, 0.3141082734275339, 0.8063905248474006, 0.5850486948450955]
+  M1  =  [ 0.36799630150742724, 0.8353476258002335, 0.57704047269419]
+  M1  =  VectorValue(M1/norm(M1))
+  M2  =  [ 0.3857610953527303, 0.024655018338846868, 0.6133770006613235]
+  M2  =  VectorValue(M2/norm(M2))
+  M3  =  [ 0.4516424464747618, 0.6609741557924332, 0.8441070681368911]
+  M3  =  VectorValue(M3/norm(M3))
+  M4  =  [0.7650638541897623, 0.8268625401770648, 0.3103412304991431]
+  M4  =  VectorValue(M4/norm(M4))
+  
+  model1 = HGO_1Fiber(c1=c1[1], c2=c2[1])
+  model2 = HGO_1Fiber(c1=c1[2], c2=c2[2])
+  model3 = HGO_1Fiber(c1=c1[3], c2=c2[3])
+  model4 = HGO_1Fiber(c1=c1[4], c2=c2[4])
+  modelMR = MooneyRivlin3D(λ=3.0, μ1=1.0, μ2=2.0)
+  modelID = IdealDielectric(ε=4.0)
+  modelelectro=modelMR+[model1 model2 model3 model4]+modelID
+ 
+
+  Ψ, ∂Ψ∂F, ∂Ψ∂F∂F = modelelectro()
+ 
+  F, _, _ = get_Kinematics(model1.Kinematic)
+  E = get_Kinematics(modelID.Kinematic)
+  @test Ψ(F(∇u),E(∇φ) ,(M1,M2,M3,M4)) == -27.513663654827152
+  @test isapprox(norm(∂Ψ∂F(F(∇u),E(∇φ) ,(M1,M2,M3,M4))) ,  47.45420724735932,rtol=1e-14)
+  @test isapprox(norm(∂Ψ∂F∂F(F(∇u),E(∇φ) ,(M1,M2,M3,M4))), 14.707913034885005,rtol=1e-14)
+end
+
+
+
 @testset "ElectroMechano" begin
   modelMR = MooneyRivlin3D(λ=3.0, μ1=1.0, μ2=2.0)
   modelID = IdealDielectric(ε=4.0)

test/TestConstitutiveModels/PhysicalModelTests.jl

@@ -39,6 +39,38 @@ function test_derivatives_3D_(model::PhysicalModel; rtol=1e-14, kwargs...)
 end
 
 
+
+
+@testset "composition of HGO_1Fiber" begin
+  c1  =  [0.6639232500447778, 0.5532987701062146, 0.9912576142028674, 0.4951942011240962]
+  c2  =  [0.800583033264982, 0.3141082734275339, 0.8063905248474006, 0.5850486948450955]
+  M1  =  [ 0.36799630150742724, 0.8353476258002335, 0.57704047269419]
+  M1  =  VectorValue(M1/norm(M1))
+  M2  =  [ 0.3857610953527303, 0.024655018338846868, 0.6133770006613235]
+  M2  =  VectorValue(M2/norm(M2))
+  M3  =  [ 0.4516424464747618, 0.6609741557924332, 0.8441070681368911]
+  M3  =  VectorValue(M3/norm(M3))
+  M4  =  [0.7650638541897623, 0.8268625401770648, 0.3103412304991431]
+  M4  =  VectorValue(M4/norm(M4))
+  
+  model1 = HGO_1Fiber(c1=c1[1], c2=c2[1])
+  model2 = HGO_1Fiber(c1=c1[2], c2=c2[2])
+  model3 = HGO_1Fiber(c1=c1[3], c2=c2[3])
+  model4 = HGO_1Fiber(c1=c1[4], c2=c2[4])
+
+  model=[model1 model2 model3 model4]
+
+  Ψ, ∂Ψ∂F, ∂Ψ∂F∂F = model()
+  ∂Ψ∂F_(F,M1,M2,M3,M4) = TensorValue(ForwardDiff.gradient(F -> Ψ(F,(get_array(M1),get_array(M2),get_array(M3),get_array(M4))), get_array(F)))
+  ∂Ψ∂F∂F_(F,M1,M2,M3,M4) = TensorValue(ForwardDiff.hessian(F -> Ψ(F,(get_array(M1),get_array(M2),get_array(M3),get_array(M4))), get_array(F)))
+
+
+  F, _, _ = get_Kinematics(model1.Kinematic)
+  @test Ψ(F(∇u3), (M1,M2,M3,M4)) == 0.0005561006012767033 
+  @test isapprox(norm(∂Ψ∂F(F(∇u3), (M1,M2,M3,M4))) ,  norm(∂Ψ∂F_(F(∇u3), M1,M2,M3,M4)),rtol=1e-14)
+  @test isapprox(norm(∂Ψ∂F∂F(F(∇u3), (M1,M2,M3,M4))), norm(∂Ψ∂F∂F_(F(∇u3), M1,M2,M3,M4)),rtol=1e-14)
+end
+
 @testset "HGO_4Fibers" begin
   c1  =  [0.6639232500447778, 0.5532987701062146, 0.9912576142028674, 0.4951942011240962]
   c2  =  [0.800583033264982, 0.3141082734275339, 0.8063905248474006, 0.5850486948450955]
@@ -64,7 +96,7 @@ end
 end
 
 
-@testset "HGO_4Fibers+MooneyRivlin3D" begin
+@testset "composition of HGO_4Fibers+MooneyRivlin3D" begin
   c1  =  [0.6639232500447778, 0.5532987701062146, 0.9912576142028674, 0.4951942011240962]
   c2  =  [0.800583033264982, 0.3141082734275339, 0.8063905248474006, 0.5850486948450955]
   M1  =  [ 0.36799630150742724, 0.8353476258002335, 0.57704047269419]