diff --git a/src/PhysicalModels/ElectricalModels.jl b/src/PhysicalModels/ElectricalModels.jl
index c4b5c20..c8f401b 100644
--- a/src/PhysicalModels/ElectricalModels.jl
+++ b/src/PhysicalModels/ElectricalModels.jl
@@ -4,8 +4,8 @@
 # ===================
 
 struct IdealDielectric <: Electro
-  ε::Float64
-  function IdealDielectric(; ε::Float64)
+  ε::Real
+  function IdealDielectric(; ε::Real)
     new(ε)
   end
 end
diff --git a/src/PhysicalModels/ElectroMechanicalModels.jl b/src/PhysicalModels/ElectroMechanicalModels.jl
index 930df47..ae2ff22 100644
--- a/src/PhysicalModels/ElectroMechanicalModels.jl
+++ b/src/PhysicalModels/ElectroMechanicalModels.jl
@@ -119,7 +119,7 @@ end
 
 struct FlexoElectroModel{EM<:ElectroMechano} <: FlexoElectro{EM}
   electromechano::EM
-  κ::Float64
+  κ::Real
 
   function FlexoElectroModel(electro::E, mechano::M; κ=1.0) where {E<:Electro,M<:Mechano}
     physmodel = ElectroMechModel(electro, mechano)
diff --git a/src/PhysicalModels/MagneticModels.jl b/src/PhysicalModels/MagneticModels.jl
index f6a316f..e7a87ac 100644
--- a/src/PhysicalModels/MagneticModels.jl
+++ b/src/PhysicalModels/MagneticModels.jl
@@ -4,17 +4,20 @@
 # ===================
 
 
-struct Magnetic <: Magneto
-  μ::Float64
-  αr::Ref{Float64}
-  χe::Float64
+struct Magnetic{T<:Real} <: Magneto
+  μ::Real
+  αr::Base.RefValue{T}
+  χe::Real
   
-function Magnetic(; μ0::Float64, αr::Float64, χe::Float64=0.0)
-  new(μ0, Ref(αr), χe)
+  function Magnetic(; μ0::Real, αr::Real, χe::Real=0.0)
+    T = typeof(αr)
+    new{T}(μ0, Ref{T}(αr), χe)
+  end
 end
+
 function (obj::Magnetic)(Λ::Float64=1.0)
-  μ, αr, χe = obj.μ, obj.αr, obj.χe
-  ℋᵣ(N) = αr[] * Λ * N
+  (; μ, χe) = obj
+  ℋᵣ(N) = obj.αr * Λ * N
   # Energy #
   Ψmm(ℋ₀, N) = (-μ / 2.0) * ((ℋ₀ + ℋᵣ(N)) ⋅ (ℋ₀ + ℋᵣ(N))) * (1 + χe)
   ∂Ψmm_∂φ(ℋ₀, N) = (-μ) * (ℋ₀ + ℋᵣ(N)) * (1 + χe)
@@ -23,18 +26,15 @@ function (obj::Magnetic)(Λ::Float64=1.0)
 end
 
 
-end
-
-
 struct IdealMagnetic <: Magneto
-  μ::Float64
-  χe::Float64
-  function IdealMagnetic(; μ0::Float64, χe::Float64=0.0)
+  μ::Real
+  χe::Real
+  function IdealMagnetic(; μ0::Real, χe::Real=0.0)
     new(μ0, χe)
   end
   function (obj::IdealMagnetic)(Λ::Float64=1.0)
 
-    μ, χe = obj.μ, obj.χe
+    (; μ, χe) = obj
     J(F) = det(F)
     H(F) = J(F) * inv(F)'
 
@@ -72,15 +72,15 @@ end
 
 
 struct IdealMagnetic2D <: Magneto
-  μ::Float64
-  χe::Float64
-  function IdealMagnetic2D(; μ0::Float64, χe::Float64=0.0)
+  μ::Real
+  χe::Real
+  function IdealMagnetic2D(; μ0::Real, χe::Real=0.0)
     new(μ0, χe)
   end
 
   function (obj::IdealMagnetic2D)(Λ::Float64=1.0)
 
-    μ, χe = obj.μ, obj.χe
+    (; μ, χe) = obj
     J(F) = det(F)
     H(F) = J(F) * inv(F)'
 
@@ -113,20 +113,21 @@ end
 
 
 struct HardMagnetic <: Magneto
-  μ::Float64
-  αr::Float64
-  χe::Float64
-  χr::Float64
-  χt::Float64
-  βmok::Float64
-  βcoup::Float64
-  function HardMagnetic(; μ0::Float64, αr::Float64, χe::Float64=0.0, χr::Float64=8.0, χt::Union{Float64,Nothing}=nothing, βmok::Float64=0.0, βcoup::Float64=0.0)
-    χt_val = isnothing(χt) ? χe : χt
-    new(μ0, αr, χe, χr, χt_val, βmok, βcoup)
+  μ::Real
+  αr::Real
+  χe::Real
+  χr::Real
+  χt::Real
+  βmok::Real
+  βcoup::Real
+
+  function HardMagnetic(; μ0::Real, αr::Real, χe::Real=0.0, χr::Real=8.0, χt::Union{Real,Nothing}=nothing, βmok::Real=0.0, βcoup::Real=0.0)
+    χt = isnothing(χt) ? χe : χt
+    new(μ0, αr, χe, χr, χt, βmok, βcoup)
   end
 
   function (obj::HardMagnetic)(Λ::Float64=1.0)
-    μ, χe = obj.μ, obj.χe
+    (; μ, χe) = obj
     J(F) = det(F)
     H(F) = J(F) * inv(F)'
 
@@ -162,22 +163,24 @@ struct HardMagnetic <: Magneto
 end
 
 
-struct HardMagnetic2D <: Magneto
-  μ::Float64
-  αr::Ref{Float64}
-  χe::Float64
-  χr::Float64
-  χt::Float64
-  βmok::Float64
-  βcoup::Float64
-  function HardMagnetic2D(; μ0::Float64, αr::Float64, χe::Float64=0.0, χr::Float64=8.0, χt::Union{Float64,Nothing}=nothing, βmok::Float64=0.0, βcoup::Float64=0.0)
-    χt_val = isnothing(χt) ? χe : χt
-    new(μ0, Ref(αr), χe, χr, χt_val, βmok, βcoup)
+struct HardMagnetic2D{T<:Real} <: Magneto
+  μ::Real
+  αr::Base.RefValue{T}
+  χe::Real
+  χr::Real
+  χt::Real
+  βmok::Real
+  βcoup::Real
+
+  function HardMagnetic2D(; μ0::Real, αr::Real, χe::Real=0.0, χr::Real=8.0, χt::Union{Real,Nothing}=nothing, βmok::Real=0.0, βcoup::Real=0.0)
+    χt = isnothing(χt) ? χe : χt
+    T = typeof(αr)
+    new{T}(μ0, Ref{T}(αr), χe, χr, χt, βmok, βcoup)
   end
 
   function (obj::HardMagnetic2D)(Λ::Float64=1.0)
 
-    μ, χe = obj.μ, obj.χe
+    (; μ, χe) = obj
     J(F) = det(F)
     H(F) = J(F) * inv(F)'
 
@@ -208,3 +211,21 @@ struct HardMagnetic2D <: Magneto
     return (Ψmm, ∂Ψmm_∂u, ∂Ψmm_∂φ, ∂Ψmm_∂uu, ∂Ψmm_∂φu, ∂Ψmm_∂φφ)
   end
 end
+
+
+function Base.getproperty(obj::Union{Magnetic,HardMagnetic2D}, prop::Symbol)
+  if prop === :αr
+    return getfield(obj, :αr)[]
+  else
+    return getfield(obj, prop)
+  end
+end
+
+
+function Base.setproperty!(obj::Union{Magnetic,HardMagnetic2D}, prop::Symbol, val)
+  if prop === :αr
+    return getfield(obj, :αr)[] = val
+  else
+    return setfield!(obj, prop, val)
+  end
+end
diff --git a/src/PhysicalModels/MagnetoMechanicalModels.jl b/src/PhysicalModels/MagnetoMechanicalModels.jl
index 9cd0250..c892d3d 100644
--- a/src/PhysicalModels/MagnetoMechanicalModels.jl
+++ b/src/PhysicalModels/MagnetoMechanicalModels.jl
@@ -68,12 +68,10 @@ function _getCoupling(mag::HardMagnetic, ::Mechano, Λ::Float64=1.0)
   # Hard magnetics in ultra-soft magnetorheological elastomers enhance fracture toughness and 
   # delay crack propagation, Journal of the Mechanics and Physics of Solids,
 
-
-  μ, αr, χr, χt, βcoup, βmok = mag.μ, mag.αr, mag.χr, mag.χt, mag.βcoup, mag.βmok
+  (; μ, χr, χt, βcoup, βmok) = mag
   J(F) = det(F)
   H(F) = det(F) * inv(F)'
  
-
   #-------------------------------------------------------------------------------------
   # FIRST TERM
   #-------------------------------------------------------------------------------------
@@ -85,7 +83,7 @@ function _getCoupling(mag::HardMagnetic, ::Mechano, Λ::Float64=1.0)
   Hℋ₀(F, ℋ₀) = H(F) * ℋ₀
   Hℋ₀Hℋ₀(F, ℋ₀) = Hℋ₀(F, ℋ₀) ⋅ Hℋ₀(F, ℋ₀)
 
-  ℋᵣ(N) = αr * Λ* N
+  ℋᵣ(N) = mag.α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))
@@ -150,7 +148,7 @@ function _getCoupling(mag::HardMagnetic2D, ::Mechano, Λ::Float64=1.0)
   # Hard magnetics in ultra-soft magnetorheological elastomers enhance fracture toughness and 
   # delay crack propagation, Journal of the Mechanics and Physics of Solids,
 
-  μ, αr, χr, χt, βcoup, βmok = mag.μ, mag.αr, mag.χr, mag.χt, mag.βcoup, mag.βmok
+  (; μ, χr, χt, βcoup, βmok) = mag
   J(F) = det(F)
   H(F) = det(F) * inv(F)'
  
@@ -164,7 +162,7 @@ function _getCoupling(mag::HardMagnetic2D, ::Mechano, Λ::Float64=1.0)
   #-------------------------------------------------------------------------------------
   Hℋ₀(F, ℋ₀) = H(F) * ℋ₀
   Hℋ₀Hℋ₀(F, ℋ₀) = Hℋ₀(F, ℋ₀) ⋅ Hℋ₀(F, ℋ₀)
-  ℋᵣ(N) = αr[] * Λ* N
+  ℋᵣ(N) = mag.α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))
diff --git a/src/PhysicalModels/MechanicalModels.jl b/src/PhysicalModels/MechanicalModels.jl
index dbf36dc..0bc9ec5 100644
--- a/src/PhysicalModels/MechanicalModels.jl
+++ b/src/PhysicalModels/MechanicalModels.jl
@@ -5,8 +5,8 @@
 # ============================================
 
 struct VolumetricEnergy <: Volumetric
-  λ::Float64
-  function VolumetricEnergy(; λ::Float64)
+  λ::Real
+  function VolumetricEnergy(; λ::Real)
     new(λ)
   end
 end
@@ -35,13 +35,13 @@ end
 
 struct HessianRegularization <: Mechano
   mechano::Mechano
-  δ::Float64
+  δ::Real
 
-  function HessianRegularization(mechano::Mechano; δ::Float64=1.0e-6)
+  function HessianRegularization(mechano::Mechano; δ::Real=1.0e-6)
     new(mechano, δ)
   end
 
-  function HessianRegularization(; mechano::Mechano, δ::Float64=1.0e-6)
+  function HessianRegularization(; mechano::Mechano, δ::Real=1.0e-6)
     new(mechano, δ)
   end
 
@@ -62,14 +62,14 @@ end
 
 struct Hessian∇JRegularization <: Mechano
   mechano::Mechano
-  δ::Float64
-  κ::Float64
+  δ::Real
+  κ::Real
 
-  function Hessian∇JRegularization(mechano::Mechano; δ::Float64=1.0e-6, κ::Float64=1.0)
+  function Hessian∇JRegularization(mechano::Mechano; δ::Real=1.0e-6, κ::Real=1.0)
     new(mechano, δ, κ)
   end
 
-  function Hessian∇JRegularization(; mechano::Mechano, δ::Float64=1.0e-6, κ::Float64=1.0)
+  function Hessian∇JRegularization(; mechano::Mechano, δ::Real=1.0e-6, κ::Real=1.0)
     new(mechano, δ, κ)
   end
 
@@ -97,11 +97,11 @@ end
 # Energy interpolations
 # ======================
 struct EnergyInterpolationScheme <: IsoElastic
-  p::Float64
+  p::Real
   model1::IsoElastic
   model2::IsoElastic
 
-  function EnergyInterpolationScheme(model1::IsoElastic, model2::IsoElastic; p::Float64=3.0)
+  function EnergyInterpolationScheme(model1::IsoElastic, model2::IsoElastic; p::Real=3.0)
     new(p, model1, model2)
   end
 
@@ -191,11 +191,11 @@ end
 # ===================
 
 struct Yeoh3D <: IsoElastic
-  λ::Float64
-  C10::Float64
-  C20::Float64
-  C30::Float64
-  function Yeoh3D(; λ::Float64, C10::Float64, C20::Float64, C30::Float64)
+  λ::Real
+  C10::Real
+  C20::Real
+  C30::Real
+  function Yeoh3D(; λ::Real, C10::Real, C20::Real, C30::Real)
     new(λ, C10, C20, C30)
   end
 
@@ -228,10 +228,10 @@ struct Yeoh3D <: IsoElastic
 end
 
 struct LinearElasticity2D <: IsoElastic
-  λ::Float64
-  μ::Float64
-  ρ::Float64
-  function LinearElasticity2D(; λ::Float64, μ::Float64, ρ::Float64=0.0)
+  λ::Real
+  μ::Real
+  ρ::Real
+  function LinearElasticity2D(; λ::Real, μ::Real, ρ::Real=0.0)
     new(λ, μ, ρ)
   end
 
@@ -247,10 +247,10 @@ end
 
 
 mutable struct LinearElasticity3D <: IsoElastic
-  λ::Float64
-  μ::Float64
-  ρ::Float64
-  function LinearElasticity3D(; λ::Float64, μ::Float64, ρ::Float64=0.0)
+  λ::Real
+  μ::Real
+  ρ::Real
+  function LinearElasticity3D(; λ::Real, μ::Real, ρ::Real=0.0)
     new(λ, μ, ρ)
   end
 
@@ -266,10 +266,10 @@ end
 
 
 struct NeoHookean3D <: IsoElastic
-  λ::Float64
-  μ::Float64
-  ρ::Float64
-  function NeoHookean3D(; λ::Float64, μ::Float64, ρ::Float64=0.0)
+  λ::Real
+  μ::Real
+  ρ::Real
+  function NeoHookean3D(; λ::Real, μ::Real, ρ::Real=0.0)
     new(λ, μ, ρ)
   end
 
@@ -292,13 +292,13 @@ end
 
 
 struct Gent2D <: IsoElastic
-  λ::Float64
-  μ::Float64
-  Jm::Float64
-  γ::Float64
-  ρ::Float64
+  λ::Real
+  μ::Real
+  Jm::Real
+  γ::Real
+  ρ::Real
 
-  function Gent2D(; λ::Float64, μ::Float64, Jm::Float64, γ::Float64, ρ::Float64=0.0)
+  function Gent2D(; λ::Real, μ::Real, Jm::Real, γ::Real, ρ::Real=0.0)
     new(λ, μ, Jm, γ, ρ)
   end
 
@@ -332,11 +332,11 @@ end
 
 
 struct MooneyRivlin3D <: IsoElastic
-  λ::Float64
-  μ1::Float64
-  μ2::Float64
-  ρ::Float64
-  function MooneyRivlin3D(; λ::Float64, μ1::Float64, μ2::Float64, ρ::Float64=0.0)
+  λ::Real
+  μ1::Real
+  μ2::Real
+  ρ::Real
+  function MooneyRivlin3D(; λ::Real, μ1::Real, μ2::Real, ρ::Real=0.0)
     new(λ, μ1, μ2, ρ)
   end
 
@@ -363,12 +363,12 @@ end
 
 
 struct MooneyRivlin2D <: IsoElastic
-  λ::Float64
-  μ1::Float64
-  μ2::Float64
-  ρ::Float64
+  λ::Real
+  μ1::Real
+  μ2::Real
+  ρ::Real
 
-  function MooneyRivlin2D(; λ::Float64, μ1::Float64, μ2::Float64, ρ::Float64=0.0)
+  function MooneyRivlin2D(; λ::Real, μ1::Real, μ2::Real, ρ::Real=0.0)
     new(λ, μ1, μ2, ρ)
   end
 
@@ -389,13 +389,13 @@ end
 
 
 struct NonlinearMooneyRivlin3D <: IsoElastic
-  λ::Float64
-  μ1::Float64
-  μ2::Float64
-  α1::Float64
-  α2::Float64
-  ρ::Float64
-  function NonlinearMooneyRivlin3D(; λ::Float64, μ1::Float64, μ2::Float64, α1::Float64, α2::Float64, ρ::Float64=0.0)
+  λ::Real
+  μ1::Real
+  μ2::Real
+  α1::Real
+  α2::Real
+  ρ::Real
+  function NonlinearMooneyRivlin3D(; λ::Real, μ1::Real, μ2::Real, α1::Real, α2::Real, ρ::Real=0.0)
     new(λ, μ1, μ2, α1, α2, ρ)
   end
 
@@ -424,13 +424,13 @@ struct NonlinearMooneyRivlin3D <: IsoElastic
 end
 
 struct NonlinearMooneyRivlin2D <: IsoElastic
-  λ::Float64
-  μ1::Float64
-  μ2::Float64
-  α1::Float64
-  α2::Float64
-  ρ::Float64
-  function NonlinearMooneyRivlin2D(; λ::Float64, μ1::Float64, μ2::Float64, α1::Float64, α2::Float64, ρ::Float64=0.0)
+  λ::Real
+  μ1::Real
+  μ2::Real
+  α1::Real
+  α2::Real
+  ρ::Real
+  function NonlinearMooneyRivlin2D(; λ::Real, μ1::Real, μ2::Real, α1::Real, α2::Real, ρ::Real=0.0)
     new(λ, μ1, μ2, α1, α2, ρ)
   end
 
@@ -461,14 +461,14 @@ end
 
 
 struct NonlinearMooneyRivlin2D_CV <: IsoElastic
-  λ::Float64
-  μ1::Float64
-  μ2::Float64
-  α1::Float64
-  α2::Float64
-  γ::Float64
-  ρ::Float64
-  function NonlinearMooneyRivlin2D_CV(; λ::Float64, μ1::Float64, μ2::Float64, α1::Float64, α2::Float64, γ::Float64, ρ::Float64=0.0)
+  λ::Real
+  μ1::Real
+  μ2::Real
+  α1::Real
+  α2::Real
+  γ::Real
+  ρ::Real
+  function NonlinearMooneyRivlin2D_CV(; λ::Real, μ1::Real, μ2::Real, α1::Real, α2::Real, γ::Real, ρ::Real=0.0)
     new(λ, μ1, μ2, α1, α2, γ, ρ)
   end
 
@@ -496,14 +496,14 @@ struct NonlinearMooneyRivlin2D_CV <: IsoElastic
 end
 
 struct NonlinearMooneyRivlin_CV <: IsoElastic
-  λ::Float64
-  μ1::Float64
-  μ2::Float64
-  α1::Float64
-  α2::Float64
-  γ::Float64
-  ρ::Float64
-  function NonlinearMooneyRivlin_CV(; λ::Float64, μ1::Float64, μ2::Float64, α1::Float64, α2::Float64, γ::Float64, ρ::Float64=0.0)
+  λ::Real
+  μ1::Real
+  μ2::Real
+  α1::Real
+  α2::Real
+  γ::Real
+  ρ::Real
+  function NonlinearMooneyRivlin_CV(; λ::Real, μ1::Real, μ2::Real, α1::Real, α2::Real, γ::Real, ρ::Real=0.0)
     new(λ, μ1, μ2, α1, α2, γ, ρ)
   end
 
@@ -536,12 +536,12 @@ end
 
 
 struct NonlinearNeoHookean_CV <: IsoElastic
-  λ::Float64
-  μ::Float64
-  α::Float64
-  γ::Float64
-  ρ::Float64
-  function NonlinearNeoHookean_CV(; λ::Float64, μ::Float64, α::Float64, γ::Float64, ρ::Float64=0.0)
+  λ::Real
+  μ::Real
+  α::Real
+  γ::Real
+  ρ::Real
+  function NonlinearNeoHookean_CV(; λ::Real, μ::Real, α::Real, γ::Real, ρ::Real=0.0)
     new(λ, μ, α, γ, ρ)
   end
 
@@ -568,12 +568,12 @@ end
 
 
 struct NonlinearIncompressibleMooneyRivlin2D_CV <: IsoElastic
-  λ::Float64
-  μ::Float64
-  α::Float64
-  γ::Float64
-  ρ::Float64
-  function NonlinearIncompressibleMooneyRivlin2D_CV(; λ::Float64, μ::Float64, α::Float64, γ::Float64, ρ::Float64=0.0)
+  λ::Real
+  μ::Real
+  α::Real
+  γ::Real
+  ρ::Real
+  function NonlinearIncompressibleMooneyRivlin2D_CV(; λ::Real, μ::Real, α::Real, γ::Real, ρ::Real=0.0)
     new(λ, μ, α, γ, ρ)
   end
 
@@ -615,9 +615,9 @@ end
 
 
 struct EightChain <: IsoElastic
-  μ::Float64
-  N::Float64
-  function EightChain(; μ::Float64, N::Float64)
+  μ::Real
+  N::Real
+  function EightChain(; μ::Real, N::Real)
     new(μ, N)
   end
 
@@ -680,11 +680,11 @@ end
 
 
 struct TransverseIsotropy3D <: AnisoElastic
-  μ::Float64
-  α1::Float64
-  α2::Float64
-  ρ::Float64
-  function TransverseIsotropy3D(; μ::Float64, α1::Float64, α2::Float64, ρ::Float64=0.0)
+  μ::Real
+  α1::Real
+  α2::Real
+  ρ::Real
+  function TransverseIsotropy3D(; μ::Real, α1::Real, α2::Real, ρ::Real=0.0)
     new(μ, α1, α2, ρ)
   end
 
@@ -715,11 +715,11 @@ end
 
 
 struct TransverseIsotropy2D <: AnisoElastic
-  μ::Float64
-  α1::Float64
-  α2::Float64
-  ρ::Float64
-  function TransverseIsotropy2D(; μ::Float64, α1::Float64, α2::Float64, ρ::Float64=0.0)
+  μ::Real
+  α1::Real
+  α2::Real
+  ρ::Real
+  function TransverseIsotropy2D(; μ::Real, α1::Real, α2::Real, ρ::Real=0.0)
     new(μ, α1, α2, ρ)
   end
 
@@ -756,12 +756,12 @@ struct TransverseIsotropy2D <: AnisoElastic
 end
 
 
-struct HGO_4Fibers <: AnisoElastic
-  c1::Vector{Float64}
-  c2::Vector{Float64}
-  function HGO_4Fibers(; c1::Vector{Float64}, c2::Vector{Float64})
+struct HGO_4Fibers{T<:Real} <: AnisoElastic
+  c1::Vector{T}
+  c2::Vector{T}
+  function HGO_4Fibers(; c1::Vector{T}, c2::Vector{T}) where {T<:Real}
     @assert length(c1) == length(c2) == 4
-    new(c1, c2)
+    new{T}(c1, c2)
   end
 
   function (obj::HGO_4Fibers)(Λ::Float64=1.0; Threshold=0.01)
@@ -807,9 +807,9 @@ end
 
 
 struct HGO_1Fiber <: AnisoElastic
-  c1::Float64
-  c2::Float64
-  function HGO_1Fiber(; c1::Float64, c2::Float64)
+  c1::Real
+  c2::Real
+  function HGO_1Fiber(; c1::Real, c2::Real)
     new(c1, c2)
   end
 
@@ -837,11 +837,11 @@ end
 
 
 struct IncompressibleNeoHookean3D <: IsoElastic
-  λ::Float64
-  μ::Float64
-  ρ::Float64
-  δ::Float64
-  function IncompressibleNeoHookean3D(; λ::Float64, μ::Float64, ρ::Float64=0.0, δ::Float64=0.1)
+  λ::Real
+  μ::Real
+  ρ::Real
+  δ::Real
+  function IncompressibleNeoHookean3D(; λ::Real, μ::Real, ρ::Real=0.0, δ::Real=0.1)
     new(λ, μ, ρ, δ)
   end
 
@@ -897,11 +897,11 @@ function SecondPiola(obj::IncompressibleNeoHookean3D, Λ::Float64=1.0)
 end
 
 struct IncompressibleNeoHookean2D <: IsoElastic
-  λ::Float64
-  μ::Float64
-  ρ::Float64
-  δ::Float64
-  function IncompressibleNeoHookean2D(; λ::Float64, μ::Float64, ρ::Float64=0.0, δ::Float64=0.1)
+  λ::Real
+  μ::Real
+  ρ::Real
+  δ::Real
+  function IncompressibleNeoHookean2D(; λ::Real, μ::Real, ρ::Real=0.0, δ::Real=0.1)
     new(λ, μ, ρ, δ)
   end
 
@@ -938,11 +938,11 @@ struct IncompressibleNeoHookean2D <: IsoElastic
 end
 
 struct IncompressibleNeoHookean2D_CV <: IsoElastic
-  λ::Float64
-  μ::Float64
-  γ::Float64
-  ρ::Float64
-  function IncompressibleNeoHookean2D_CV(; λ::Float64, μ::Float64, γ::Float64, ρ::Float64=0.0)
+  λ::Real
+  μ::Real
+  γ::Real
+  ρ::Real
+  function IncompressibleNeoHookean2D_CV(; λ::Real, μ::Real, γ::Real, ρ::Real=0.0)
     new(λ, μ, γ, ρ)
   end
 
@@ -970,10 +970,10 @@ end
 
 
 struct ARAP2D_regularized <: IsoElastic
-  μ::Float64
-  ρ::Float64
-  δ::Float64
-  function ARAP2D_regularized(; μ::Float64, ρ::Float64=0.0, δ::Float64=0.1)
+  μ::Real
+  ρ::Real
+  δ::Real
+  function ARAP2D_regularized(; μ::Real, ρ::Real=0.0, δ::Real=0.1)
     new(μ, ρ, δ)
   end
 
@@ -1011,9 +1011,9 @@ end
 
 
 struct ARAP2D <: IsoElastic
-  μ::Float64
-  ρ::Float64
-  function ARAP2D(; μ::Float64, ρ::Float64=0.0)
+  μ::Real
+  ρ::Real
+  function ARAP2D(; μ::Real, ρ::Real=0.0)
     new(μ, ρ)
   end
 
@@ -1038,7 +1038,7 @@ end
 
 
 struct IsochoricNeoHookean3D <: IsoElastic
-  μ::Float64
+  μ::Real
 end
 
 function IsochoricNeoHookean3D(; μ::Real)
@@ -1083,12 +1083,12 @@ end
 
 
 struct IncompressibleNeoHookean3D_2dP <: Mechano
-  μ::Float64
-  τ::Float64
-  Δt::Float64
-  ρ::Float64
+  μ::Real
+  τ::Real
+  Δt::Real
+  ρ::Real
 
-  function IncompressibleNeoHookean3D_2dP(; μ::Float64, τ::Float64, Δt::Float64, ρ::Float64=0.0)
+  function IncompressibleNeoHookean3D_2dP(; μ::Real, τ::Real, Δt::Real, ρ::Real=0.0)
     new(μ, τ, Δt, ρ)
   end
 
diff --git a/src/PhysicalModels/PINNs.jl b/src/PhysicalModels/PINNs.jl
index 54b076f..c0a3aea 100644
--- a/src/PhysicalModels/PINNs.jl
+++ b/src/PhysicalModels/PINNs.jl
@@ -35,9 +35,9 @@ struct ThermoElectroMech_PINNs{A,B,C,D} <: ThermoElectroMechano{Thermo,Electro,M
   ϵ::C
   β::D
   nLayer::Int64
-  κ::Float64
+  κ::Real
 
-  function ThermoElectroMech_PINNs(; W, b, ϵ, β, nLayer::Int64, κ::Float64=1e-4)
+  function ThermoElectroMech_PINNs(; W, b, ϵ, β, nLayer::Int64, κ::Real=1e-4)
     A, B, C, D = typeof(W), typeof(b), typeof(ϵ), typeof(β)
     new{A,B,C,D}(W, b, ϵ, β, nLayer, κ)
   end
diff --git a/src/PhysicalModels/ThermalModels.jl b/src/PhysicalModels/ThermalModels.jl
index 0a005ac..ba6fb7e 100644
--- a/src/PhysicalModels/ThermalModels.jl
+++ b/src/PhysicalModels/ThermalModels.jl
@@ -4,11 +4,11 @@
 # ===================
 
 struct ThermalModel <: Thermo
-  Cv::Float64
-  θr::Float64
-  α::Float64
-  κ::Float64
-  function ThermalModel(; Cv::Float64, θr::Float64, α::Float64, κ::Float64=10.0)
+  Cv::Real
+  θr::Real
+  α::Real
+  κ::Real
+  function ThermalModel(; Cv::Real, θr::Real, α::Real, κ::Real=10.0)
     new(Cv, θr, α, κ)
   end
 end
@@ -26,8 +26,8 @@ end
 # ===================
 
 struct EntropicElasticityLaw <: ThermalLaw
-  θr::Float64
-  γ::Float64
+  θr::Real
+  γ::Real
   EntropicElasticityLaw(; θr, γ) = new(θr, γ)
 end
 
@@ -40,9 +40,9 @@ function (law::EntropicElasticityLaw)()
 end
 
 struct NonlinearMeltingLaw <: ThermalLaw
-  θr::Float64
-  θM::Float64
-  γ::Float64
+  θr::Real
+  θM::Real
+  γ::Real
   NonlinearMeltingLaw(; θr, θM, γ) = new(θr, θM, γ)
 end
 
@@ -55,10 +55,10 @@ function (law::NonlinearMeltingLaw)()
 end
 
 struct NonlinearSofteningLaw <: ThermalLaw
-  θr::Float64
-  θt::Float64
-  γ::Float64
-  δ::Float64
+  θr::Real
+  θt::Real
+  γ::Real
+  δ::Real
   NonlinearSofteningLaw(; θr, θt, γ, δ=0) = new(θr, θt, γ, δ)
 end
 
@@ -73,8 +73,8 @@ function (law::NonlinearSofteningLaw)()
 end
 
 struct TrigonometricLaw <: ThermalLaw
-  θr::Float64
-  θM::Float64
+  θr::Real
+  θM::Real
 end
 
 function (law::TrigonometricLaw)()
@@ -89,10 +89,10 @@ function (law::TrigonometricLaw)()
 end
 
 struct PolynomialLaw <: ThermalLaw
-  θr::Float64
-  a::Float64
-  b::Float64
-  c::Float64
+  θr::Real
+  a::Real
+  b::Real
+  c::Real
 end
 
 function (law::PolynomialLaw)()
diff --git a/src/PhysicalModels/ThermoElectroMechanicalModels.jl b/src/PhysicalModels/ThermoElectroMechanicalModels.jl
index 1fff381..f49bcf6 100644
--- a/src/PhysicalModels/ThermoElectroMechanicalModels.jl
+++ b/src/PhysicalModels/ThermoElectroMechanicalModels.jl
@@ -61,13 +61,13 @@ struct ThermoElectroMech_Govindjee{T<:Thermo,E<:Electro,M<:Mechano} <: ThermoEle
   dfdθ::Function
   gθ::Function
   dgdθ::Function
-  β::Float64
+  β::Real
 
-  function ThermoElectroMech_Govindjee(thermo::T, electro::E, mechano::M; fθ::Function, dfdθ::Function, gθ::Function, dgdθ::Function, β::Float64=0.0) where {T<:Thermo,E<:Electro,M<:Mechano}
+  function ThermoElectroMech_Govindjee(thermo::T, electro::E, mechano::M; fθ::Function, dfdθ::Function, gθ::Function, dgdθ::Function, β::Real=0.0) where {T<:Thermo,E<:Electro,M<:Mechano}
     new{T,E,M}(thermo, electro, mechano, fθ, dfdθ, gθ, dgdθ, β)
   end
 
-  function ThermoElectroMech_Govindjee(; thermo::T, electro::E, mechano::M, fθ::Function, dfdθ::Function, gθ::Function, dgdθ::Function, β::Float64=0.0) where {T<:Thermo,E<:Electro,M<:Mechano}
+  function ThermoElectroMech_Govindjee(; thermo::T, electro::E, mechano::M, fθ::Function, dfdθ::Function, gθ::Function, dgdθ::Function, β::Real=0.0) where {T<:Thermo,E<:Electro,M<:Mechano}
     new{T,E,M}(thermo, electro, mechano, fθ, dfdθ, gθ, dgdθ, β)
   end
 
diff --git a/src/PhysicalModels/ThermoMechanicalModels.jl b/src/PhysicalModels/ThermoMechanicalModels.jl
index 5ad7af1..337cc8f 100644
--- a/src/PhysicalModels/ThermoMechanicalModels.jl
+++ b/src/PhysicalModels/ThermoMechanicalModels.jl
@@ -151,16 +151,16 @@ end
 struct ThermoMech_EntropicPolyconvex{T<:Thermo,M<:Mechano} <: ThermoMechano{T,M}
   thermo::T
   mechano::M
-  β::Float64
+  β::Real
   G::Function
   ϕ::Function
   s::Function
 
-  function ThermoMech_EntropicPolyconvex(thermo::T, mechano::M; β::Float64, G::Function, ϕ::Function, s::Function) where {T <: Thermo, M <: Mechano}
+  function ThermoMech_EntropicPolyconvex(thermo::T, mechano::M; β::Real, G::Function, ϕ::Function, s::Function) where {T <: Thermo, M <: Mechano}
     new{T,M}(thermo, mechano, β, G, ϕ, s)
   end
 
-  function ThermoMech_EntropicPolyconvex(; thermo::T, mechano::M, β::Float64, G::Function, ϕ::Function, s::Function) where {T <: Thermo, M <: Mechano}
+  function ThermoMech_EntropicPolyconvex(; thermo::T, mechano::M, β::Real, G::Function, ϕ::Function, s::Function) where {T <: Thermo, M <: Mechano}
     new{T,M}(thermo, mechano, β, G, ϕ, s)
   end
 
diff --git a/src/PhysicalModels/ViscousModels.jl b/src/PhysicalModels/ViscousModels.jl
index 4113aa9..5b10b07 100644
--- a/src/PhysicalModels/ViscousModels.jl
+++ b/src/PhysicalModels/ViscousModels.jl
@@ -8,9 +8,9 @@ using ..TensorAlgebra
 
 struct ViscousIncompressible <: Visco
   elasto::Elasto
-  τ::Float64
-  Δt::Ref{Float64}
-  function ViscousIncompressible(elasto; τ::Float64)
+  τ::Real
+  Δt::Ref{Real}
+  function ViscousIncompressible(elasto; τ::Real)
     new(elasto, τ, 0)
   end
   function (obj::ViscousIncompressible)()
@@ -96,7 +96,7 @@ end
 struct GeneralizedMaxwell{E<:Elasto} <: ViscoElastic{E}
   longterm::E
   branches::NVisco{N} where N
-  Δt::Ref{Float64}
+  Δt::Ref{Real}
   function GeneralizedMaxwell(longTerm::E, branches::Vararg{Visco}) where {E<:Elasto}
     new{E}(longTerm,NVisco(branches),0)
   end
diff --git a/test/TestConstitutiveModels/PhysicalModelTests.jl b/test/TestConstitutiveModels/PhysicalModelTests.jl
index 3ef7a49..ef0a6cb 100644
--- a/test/TestConstitutiveModels/PhysicalModelTests.jl
+++ b/test/TestConstitutiveModels/PhysicalModelTests.jl
@@ -787,6 +787,19 @@ end
   @test norm(∂Ψφφ(F(∇u), H0(∇φ), N)) == 2.1878750641250348e-6
 end
 
+@testset "Magnetic_αr_mutability" begin
+  model = Magnetic(μ0=1.2566e-6, αr=4e-2, χe=0.0)
+  @test model.αr == 4e-2
+  model.αr = 1e-2
+  @test model.αr == 1e-2
+  Ψ, _... = model()
+  H0 = VectorValue(rand(3))
+  N  = VectorValue(rand(3))
+  Ψ1 = Ψ(H0, N)
+  model.αr *= 2
+  Ψ2 = Ψ(H0, N)
+  @test Ψ2 < Ψ1
+end
 
 @testset "Magnetic3D" begin
   #  Memory estimate: 0 bytes, allocs estimate: 0.