Merge branch 'main' into fix-mpi

Author: miguelmaso

Project.toml

@@ -1,7 +1,7 @@
 name = "HyperFEM"
 uuid = "ddbf423a-3df4-409e-9685-3dc416cae46c"
+version = "0.0.2"
 authors = ["MultiSimo_Group"]
-version = "1.0.0-DEV"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
@@ -28,3 +28,25 @@ PartitionedArrays = "5a9dfac6-5c52-46f7-8278-5e2210713be9"
 
 [extensions]
 HyperFEMDistributed = ["GridapDistributed", "GridapPETSc", "PartitionedArrays"]
+
+[compat]
+AbstractTrees = "0.4.5"
+BlockArrays = "1.9.3"
+DrWatson = "2.19.1"
+ForwardDiff = "1.3.2"
+Gridap = "0.19.7"
+GridapDistributed = "0.4.10"
+GridapGmsh = "0.7.3"
+GridapPETSc = "0.5.6"
+GridapSolvers = "0.6.1"
+IterativeSolvers = "0.9.4"
+JSON = "1.4.0"
+PartitionedArrays = "0.3.5"
+ProfileView = "1.10.2"
+Roots = "2.2.10"
+StaticArrays = "1.9.16"
+TimerOutputs = "0.5.29"
+WriteVTK = "1.21.2"
+LinearAlgebra = "1.10.0"
+Profile = "1.10.0"
+julia = "1.10"

README.md

@@ -25,7 +25,7 @@ src="https://github.com/jmartfrut/HyperFEM/blob/main/docs/imgs/logo.png?raw=true
 ## Installation
 Open the Julia REPL, type `]` to enter package mode, and install as follows
 ```julia
-pkg> add https://github.com/MultiSimOLab/HyperFEM
+pkg> add HyperFEM
 ```
 
 ## Usage
@@ -99,7 +99,7 @@ U = MultiFieldFESpace([Uu, Uφ])
 km=Kinematics(Mechano,Solid)
 ke=Kinematics(Electro,Solid)
 
-F,_,_ = get_Kinematics(kM)
+F,_,_ = get_Kinematics(km)
 E     = get_Kinematics(ke)
 
 # residual and jacobian function of load factor
@@ -154,7 +154,7 @@ src="https://github.com/MultiSimOLab/HyperFEM/blob/main/docs/imgs/sims_.png?raw=
 
 In order to give credit to the HyperFEM contributors, we ask that you please reference the paper:
 
-C. Perez‐Garcia, R. Ortigosa, J. Martínez‐Frutos, and D. Garcia‐Gonzalez, **Topology and material optimization in ultra-soft magnetoactive structures: making advantage of residual anisotropies.** Adv. Mater. (2025): e18489. https://https://doi.org/10.1002/adma.202518489
+C. Perez‐Garcia, R. Ortigosa, J. Martínez‐Frutos, and D. Garcia‐Gonzalez, **Topology and material optimization in ultra-soft magnetoactive structures: making advantage of residual anisotropies.** Adv. Mater. (2025): e18489. https://doi.org/10.1002/adma.202518489
 
 
 along with the required citations for [Gridap](https://github.com/gridap/Gridap.jl).

src/ComputationalModels/BoundaryConditions.jl

@@ -67,7 +67,7 @@ function _get_bc_func(tags_::Vector{String}, values_, bc_timesteps)
     @inbounds for i in eachindex(tags_)
 
         if values_[i] === DirichletCoupling || typeof(values_[i]) <: DirichletCoupling
-            bc_func_[i] = (Λ) -> (x) -> x 
+            bc_func_[i] = (Λ) -> (x) -> x
         else
             if isnothing(bc_timesteps[i])
                 bc_func_[i] = values_[i]
@@ -97,6 +97,16 @@ struct DirichletBC{A} <: BoundaryCondition
         caches = (bc_values)
         new{typeof(caches)}(tags_, funcs_, bc_timesteps, caches)
     end
+
+
+    function DirichletBC(bc_tags::Vector{String}, bc_values)
+        @assert(length(bc_tags) == length(bc_values))
+        bc_timesteps = [Λ -> 1.0 for _ in 1:length(bc_tags)]
+        tags_, funcs_ = _get_bc_func(bc_tags, bc_values, bc_timesteps)
+        caches = (bc_values)
+        new{typeof(caches)}(tags_, funcs_, bc_timesteps, caches)
+    end
+
 end
 
 function updateBC!(m::DirichletBC, bc_values, bc_timesteps)
@@ -206,7 +216,7 @@ struct InterpolableBC{A,B,C} <: DirichletCoupling
         dim = length(U.space.fe_dof_basis.trian.model.grid.node_coordinates[1])
         mask = U.space.dirichlet_dof_tag .== dcmask[1]
         vals = U.dirichlet_values[mask]
-        Interface_coords_ = reshape(vals, dim,:)'
+        Interface_coords_ = reshape(vals, dim, :)'
         coords = VectorValue.(eachrow(Interface_coords_))
         v = evaluate(Interpolable(1.0), coords)
         bc_values = reduce(vcat, map(x -> get_array(x), v))
@@ -228,4 +238,3 @@ function InterpolableBC!(U::TrialFESpace, bc::DirichletBC, interface_tags::Strin
     return obj
 end
 
- 

src/ComputationalModels/ComputationalModels.jl

@@ -78,5 +78,5 @@ export reset!
 export interpolate_L2_tensor
 export interpolate_L2_vector
 export interpolate_L2_scalar
-
+export L2_Projection
 end

src/ComputationalModels/Drivers.jl

@@ -43,6 +43,7 @@ end
 function solve!(m::StaggeredModel;
     stepping=(nsteps=20, nsubsteps=1, maxbisec=15),
     presolver=(τ, ∆τ) -> nothing,
+    evolτ=t-> max(min(t, 1.0), 0.0),
     kargsolve)
 
     nsubsteps = stepping[:nsubsteps]
@@ -59,7 +60,10 @@ function solve!(m::StaggeredModel;
         println("           Staggered Step: $τ              ")
         println("*******************************************")
         presolver(τ, ∆τ)
-        stevol(Λ) = ∆τ * (Λ + τ)
+        τ⁻=evolτ(τ*∆τ)
+        τ⁺=evolτ((τ+1)*∆τ)
+        stevol(Λ) = τ⁻+ (τ⁺-τ⁻)*Λ
+        # stevol(Λ) = ∆τ * (Λ + τ)
         map(x -> updateBC!(x.dirichlet, x.dirichlet.caches, [stevol for _ in 1:length(x.dirichlet.caches)]), m.compmodels)
         for τ_inner in 1:nsubsteps
             map((x) -> TrialFESpace!(x.spaces[1], x.dirichlet, 1.0), m.compmodels)
@@ -141,7 +145,7 @@ function solve!(m::StaticNonlinearModel;
     x⁻ .= x
     α = 1.0
 
-    while Λ < 1.0
+    while !isapprox(Λ, 1.0; atol=1e-12)
         Λ += ∆Λ
         if Λ > 1.0
             ∆Λ = 1.0 - (Λ - ∆Λ)

src/Exports.jl

@@ -17,6 +17,9 @@ end
 @publish TensorAlgebra logreg
 @publish TensorAlgebra Box
 @publish TensorAlgebra Ellipsoid
+@publish TensorAlgebra I2
+@publish TensorAlgebra I3
+@publish TensorAlgebra I4
 @publish TensorAlgebra I9
 @publish TensorAlgebra Tensorize
 
@@ -25,6 +28,7 @@ end
 @publish PhysicalModels LinearElasticity3D
 @publish PhysicalModels LinearElasticity2D
 @publish PhysicalModels Yeoh3D
+@publish PhysicalModels Gent2D
 @publish PhysicalModels NeoHookean3D
 @publish PhysicalModels IncompressibleNeoHookean3D
 @publish PhysicalModels IncompressibleNeoHookean2D
@@ -43,6 +47,7 @@ end
 @publish PhysicalModels TransverseIsotropy3D
 @publish PhysicalModels TransverseIsotropy2D
 @publish PhysicalModels ThermalModel
+@publish PhysicalModels ThermalModel3rdLaw
 @publish PhysicalModels IdealDielectric
 @publish PhysicalModels Magnetic
 @publish PhysicalModels IdealMagnetic
@@ -56,7 +61,6 @@ end
 @publish PhysicalModels FlexoElectroModel
 @publish PhysicalModels ThermoElectroMech_Govindjee
 @publish PhysicalModels ThermoElectroMech_PINNs
-@publish PhysicalModels ThermoElectroMech_Bonet
 @publish PhysicalModels MagnetoMechModel
 @publish PhysicalModels ARAP2D
 @publish PhysicalModels ARAP2D_regularized
@@ -135,7 +139,8 @@ end
 @publish ComputationalModels  evaluate!
 @publish ComputationalModels  InterpolableBC
 @publish ComputationalModels  InterpolableBC!
-@publish ComputationalModels TrialFESpace! # Exporting internal function of Gridap
+@publish ComputationalModels  TrialFESpace! # Exporting internal function of Gridap
+@publish ComputationalModels  L2_Projection
 
 # Note: the files FaceLabeling, CartesianTags and Evolution functions should be moved to a module different than ComputationalModels
 @publish ComputationalModels add_tag_from_vertex_filter!

src/PhysicalModels/ElectroMechanicalModels.jl

@@ -12,8 +12,8 @@ struct ElectroMechModel{E<:Electro,M<:Mechano} <: ElectroMechano
   end
 
   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, Λ)
+    Ψ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)
     ∂Ψu(F, E) = ∂Ψm_u(F) + ∂Ψem_u(F, E)
     ∂Ψφ(F, E) = ∂Ψem_φ(F, E)
@@ -23,8 +23,8 @@ struct ElectroMechModel{E<:Electro,M<:Mechano} <: ElectroMechano
     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, Λ)
+    Ψ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)
@@ -35,7 +35,7 @@ struct ElectroMechModel{E<:Electro,M<:Mechano} <: ElectroMechano
   end
   function (obj::ElectroMechModel{<:Electro,<:ViscoElastic{<: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, Λ)
+    Ψem, ∂Ψem_u, ∂Ψem_φ, ∂Ψem_uu, ∂Ψem_φu, ∂Ψem_φφ = _getCoupling(obj.electro, obj.mechano)
     Ψ(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)
@@ -44,7 +44,7 @@ struct ElectroMechModel{E<:Electro,M<:Mechano} <: ElectroMechano
     ∂Ψφφ(F, E, Fn, A...) = ∂Ψem_φφ(F, E)
     return (Ψ, ∂Ψu, ∂Ψφ, ∂Ψuu, ∂Ψφu, ∂Ψφφ)
   end
-  function (obj::ElectroMechModel{<:Electro,<:ViscoElastic{<:AnisoElastic}})()
+  function (obj::ElectroMechModel{<:Electro,<:ViscoElastic{<: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, Fn, A...) = Ψm(F, n, Fn, A...) + Ψem(F, E)
@@ -57,8 +57,6 @@ struct ElectroMechModel{E<:Electro,M<:Mechano} <: ElectroMechano
   end
 end
 
-ViscoElectricModel = ElectroMechModel{<:Electro,<:ViscoElastic}
-
 function update_time_step!(obj::ElectroMechModel, Δt::Float64)
   update_time_step!(obj.electro, Δt)
   update_time_step!(obj.mechano, Δt)
@@ -72,8 +70,13 @@ function update_state!(obj::ElectroMechModel, state, F, E, args...)
   update_state!(obj.mechano, state, F, args...)
 end
 
+function Dissipation(obj::ElectroMechModel)
+  Dvis = Dissipation(obj.mechano)
+  D(F, E, X...) = Dvis(F, X...)
+end
+
 
-function _getCoupling(elec::Electro, mec::Mechano, Λ::Float64)
+function _getCoupling(elec::Electro, mec::Mechano, Λ::Float64=0.0)
   J(F) = det(F)
   H(F) = det(F) * inv(F)'
   # Energy #

src/PhysicalModels/MagneticModels.jl

@@ -162,13 +162,13 @@ end
 
 struct HardMagnetic2D <: Magneto
   μ::Float64
-  αr::Float64
+  αr::Ref{Float64}
   χe::Float64
   χr::Float64
   βmok::Float64
   βcoup::Float64
   function HardMagnetic2D(; μ0::Float64, αr::Float64, χe::Float64=0.0, χr::Float64=8.0, βmok::Float64=0.0, βcoup::Float64=0.0)
-    new(μ0, αr, χe, χr, βmok, βcoup)
+    new(μ0, Ref(αr), χe, χr, βmok, βcoup)
   end
 
   function (obj::HardMagnetic2D)(Λ::Float64=1.0)

src/PhysicalModels/MagnetoMechanicalModels.jl

@@ -164,7 +164,7 @@ function _getCoupling(mag::HardMagnetic2D, ::Mechano, Λ::Float64=1.0)
   #-------------------------------------------------------------------------------------
   Hℋ₀(F, ℋ₀) = H(F) * ℋ₀
   Hℋ₀Hℋ₀(F, ℋ₀) = Hℋ₀(F, ℋ₀) ⋅ Hℋ₀(F, ℋ₀)
-  ℋᵣ(N) = αr * Λ* N
+  ℋᵣ(N) = αr[] * Λ* N
   Fℋᵣ(F, N) = F * ℋᵣ(N)
   Ψcoup(F, N) = (μ * J(F)) * (Fℋᵣ(F, N) ⋅ Fℋᵣ(F, N) - ℋᵣ(N) ⋅ ℋᵣ(N))
   ∂Ψcoup_∂F(F, N) = 2 * (μ * J(F)) * (Fℋᵣ(F, N) ⊗ ℋᵣ(N))