Merge branch 'main' into thermo-visco

Author: miguelmaso

src/ComputationalModels/BoundaryConditions.jl

@@ -10,6 +10,8 @@ struct MultiFieldBC <: BoundaryCondition
     BoundaryCondition::Vector{BoundaryCondition}
 end
 
+getindex(bc::MultiFieldBC, i) = bc.BoundaryCondition[i]
+
 include("EvolutionFunctions.jl")
 include("CartesianTags.jl")
 

src/ComputationalModels/ComputationalModels.jl

@@ -34,6 +34,7 @@ using WriteVTK
 using GridapGmsh 
 using GridapGmsh: GmshDiscreteModel
 
+import Base.getindex
 
 include("BoundaryConditions.jl")
 export DirichletBC

src/ComputationalModels/FESpaces.jl

@@ -44,30 +44,23 @@ end
 
 function Gridap.FESpaces.TrialFESpace!(space::MultiFieldFESpace, bc::MultiFieldBC, Λ::Float64)
   @inbounds for (i, space) in enumerate(space.spaces)
-    TrialFESpace!(space, bc.BoundaryCondition[i], Λ)
+    TrialFESpace!(space, bc[i], Λ)
   end
 end
 
 function Gridap.FESpaces.TrialFESpace!(space::MultiFieldFESpace, bc::MultiFieldBC, Λ::Float64, ΔΛ::Float64)
   @inbounds for (i, space) in enumerate(space.spaces)
-    TrialFESpace!(space, bc.BoundaryCondition[i], Λ, ΔΛ)
+    TrialFESpace!(space, bc[i], Λ, ΔΛ)
   end
 end
 
-function Gridap.FESpaces.TrialFESpace(space::SingleFieldFESpace, bc::DirichletBC)
-    TrialFESpace(space, bc, 0.0)
-end
-
-function Gridap.FESpaces.TrialFESpace(space::MultiFieldFESpace, bc::MultiFieldBC)
-    TrialFESpace(space, bc, 0.0)
-end
-
 
-function Gridap.FESpaces.TrialFESpace(space::SingleFieldFESpace, ::NothingBC, Λ::Float64)
+function Gridap.FESpaces.TrialFESpace(space::SingleFieldFESpace, ::NothingBC, Λ::Float64=0.0)
   space
 end
 
-function Gridap.FESpaces.TrialFESpace(space::SingleFieldFESpace, bc::DirichletBC, Λ::Float64)
+
+function Gridap.FESpaces.TrialFESpace(space::SingleFieldFESpace, bc::DirichletBC, Λ::Float64=0.0)
     trialspace= TrialFESpace(space, map(f -> f(Λ), bc.values))
    @inbounds for i in eachindex(bc.tags)
     if bc.caches[i] isa InterpolableBC
@@ -78,13 +71,13 @@ function Gridap.FESpaces.TrialFESpace(space::SingleFieldFESpace, bc::DirichletBC
  return trialspace
 end
 
-function Gridap.FESpaces.TrialFESpace(space::MultiFieldFESpace, bc::MultiFieldBC, Λ::Float64)
+
+function Gridap.FESpaces.TrialFESpace(space::MultiFieldFESpace, bc::MultiFieldBC, Λ::Float64=0.0)
   U_ = Vector{Union{TrialFESpace,UnconstrainedFESpace}}(undef, length(space))
   @inbounds for (i, space) in enumerate(space.spaces)
-    U_[i] = TrialFESpace(space, bc.BoundaryCondition[i], Λ)
+    U_[i] = TrialFESpace(space, bc[i], Λ)
   end
   return MultiFieldFESpace(U_)
-
 end
 
 

src/PhysicalModels/MagnetoMechanicalModels.jl

@@ -15,7 +15,7 @@ struct MagnetoMechModel{MG<:Magneto,MC<:Mechano} <: MagnetoMechano
     new{MG,MC}(magneto, mechano)
   end
 
-  function (obj::MagnetoMechModel)(Λ::Float64=1.0)
+  function (obj::MagnetoMechModel{<:Magneto,<:IsoElastic})(Λ::Float64=1.0)
     Ψm, ∂Ψm_u, ∂Ψm_uu = obj.mechano(Λ)
     Ψmm, ∂Ψmm_u, ∂Ψmm_φ, ∂Ψmm_uu, ∂Ψmm_φu, ∂Ψmm_φφ = _getCoupling(obj.magneto, obj.mechano, Λ)
 
@@ -28,9 +28,26 @@ struct MagnetoMechModel{MG<:Magneto,MC<:Mechano} <: MagnetoMechano
 
     return (Ψ, ∂Ψu, ∂Ψφ, ∂Ψuu, ∂Ψφu, ∂Ψφφ)
   end
+
+    function (obj::MagnetoMechModel{<:Magneto,<:AnisoElastic})(Λ::Float64=1.0)
+    Ψm, ∂Ψm_u, ∂Ψm_uu = obj.mechano(Λ)
+    Ψmm, ∂Ψmm_u, ∂Ψmm_φ, ∂Ψmm_uu, ∂Ψmm_φu, ∂Ψmm_φφ = _getCoupling(obj.magneto, obj.mechano, Λ)
+
+    Ψ(F, ℋ₀, N) = Ψm(F,N) + Ψmm(F, ℋ₀, N)
+    ∂Ψu(F, ℋ₀, N) = ∂Ψm_u(F,N) + ∂Ψmm_u(F, ℋ₀, N)
+    ∂Ψφ(F, ℋ₀, N) = ∂Ψmm_φ(F, ℋ₀, N)
+    ∂Ψuu(F, ℋ₀, N) = ∂Ψm_uu(F,N) + ∂Ψmm_uu(F, ℋ₀, N)
+    ∂Ψφu(F, ℋ₀, N) = ∂Ψmm_φu(F, ℋ₀, N)
+    ∂Ψφφ(F, ℋ₀, N) = ∂Ψmm_φφ(F, ℋ₀, N)
+
+    return (Ψ, ∂Ψu, ∂Ψφ, ∂Ψuu, ∂Ψφu, ∂Ψφφ)
+  end
+
+
 end
 
 
+
 function _getCoupling(mag::Union{IdealMagnetic,IdealMagnetic2D}, ::Mechano, Λ::Float64)
   Ψmm, ∂Ψmm_∂u, ∂Ψmm_∂φ, ∂Ψmm_∂uu, ∂Ψmm_∂φu, ∂Ψmm_∂φφ = mag(Λ)
 

test/TestConstitutiveModels/PhysicalModelTests.jl

@@ -644,6 +644,47 @@ end
 end
 
 
+@testset "HardMagnetic_SoftMaterial3D_aniso" begin
+
+  ∇u = TensorValue(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0) * 1e-3
+  ∇φ = VectorValue(1.0, 2.0, 3.0)
+  N = VectorValue(0.0, 0.0, 1.0)
+
+  modelMRiso = MooneyRivlin3D(λ=3.0, μ1=1.0, μ2=2.0)
+  modelMRaniso = MooneyRivlin3D(λ=3.0, μ1=1.0, μ2=2.0)
+  modelMR=modelMRiso+modelMRaniso
+  modelID = HardMagnetic(μ0=1.2566e-6, αr=40e-3, χe=0.0, χr=8.0;  βmok=1.0, βcoup=1.0)
+  modelmagneto=modelMR+modelID
+  Ψ, ∂Ψu, ∂Ψφ, ∂Ψuu, ∂Ψφu, ∂Ψφφ = modelmagneto()
+  K=Kinematics(Mechano,Solid)
+  F, _, _ = get_Kinematics(K)
+  Km=Kinematics(Magneto,Solid)
+  H0 = get_Kinematics(Km)
+
+
+  # ∂Ψu_(F) =TensorValue(ForwardDiff.gradient(x -> Ψ(x, get_array( H0(∇φ)),get_array(N) ), get_array(F)))
+  # ∂Ψφ_(H) =VectorValue(ForwardDiff.gradient(x -> Ψ(get_array( F(∇u)), x,get_array(N) ), get_array(H)))
+  # ∂Ψuu_(F) =TensorValue(ForwardDiff.hessian(x -> Ψ(x, get_array( H0(∇φ)),get_array(N) ), get_array(F)))
+  # ∂Ψφu_(H) =TensorValue(ForwardDiff.jacobian(x -> ∂Ψφ(x, get_array( H0(∇φ)),get_array(N) ), get_array(F(∇u))))
+  # ∂Ψφφ_(H) =TensorValue(ForwardDiff.jacobian(x -> ∂Ψφ(get_array( F(∇u)), x,get_array(N) ), get_array(H)))
+
+  # norm(∂Ψu_(F(∇u))) 
+  # norm(∂Ψφ_(H0(∇φ))) 
+  # norm(∂Ψuu_(F(∇u))) 
+  # norm(∂Ψφu_(H0(∇φ))) 
+  # norm(∂Ψφφ_(H0(∇φ))) 
+
+
+  @test Ψ(F(∇u), H0(∇φ), N) == 0.003187725760804994 
+  @test norm(∂Ψu(F(∇u), H0(∇φ), N)) == 0.4966662732306754
+  @test norm(∂Ψφ(F(∇u), H0(∇φ), N)) == 4.660348298920368e-6 
+  @test norm(∂Ψuu(F(∇u), H0(∇φ), N)) ==60.735729041745294 
+  @test norm(∂Ψφu(F(∇u), H0(∇φ), N)) == 1.2369035467980284e-5 
+  @test norm(∂Ψφφ(F(∇u), H0(∇φ), N)) ==  2.1878750641250348e-6
+end
+
+
+
 
 @testset "HardMagnetic_SoftMaterial3D" begin