Merge pull request #441 from ONSAS/mvanzulli/437

In place linear stress tensor (#437)

src/Entities/Tetrahedrons.jl

@@ -55,12 +55,17 @@ struct TetrahedronCache{T,ST<:Symmetric{T}} <: AbstractElementCache
     fint::Vector{T}
     "Stiffness matrix."
     Ks::ST
+    "Cosserat stress."
+    S::ST
+    "Constitutive driver."
+    ∂S∂E::Matrix{T}
     "Piola stress."
-    σ::ST
+    P::ST
     "Cauchy-Green strain."
     ε::ST
     "Deformation gradient."
     F::Matrix{T}
+    "U material derivative."
     H::Matrix{T}
     "Reference coordinates."
     X::Matrix{T}
@@ -74,10 +79,16 @@ struct TetrahedronCache{T,ST<:Symmetric{T}} <: AbstractElementCache
     aux_geometric_Ks::Matrix{T}
     "Lagrange Green Strain"
     E::ST
+    "Aux eye matrix"
+    I₃₃::Matrix{T}
+    "Aux ones matrix"
+    ones₃₃::Matrix{T}
     function TetrahedronCache()
         fint = zeros(12)
         Ks = Symmetric(zeros(12, 12))
-        σ = Symmetric(zeros(3, 3))
+        S = Symmetric(zeros(3, 3))
+        ∂S∂E = zeros(6, 6)
+        P = Symmetric(zeros(3, 3))
         ε = Symmetric(zeros(3, 3))
         F = zeros(3, 3)
         H = zeros(3, 3)
@@ -87,8 +98,10 @@ struct TetrahedronCache{T,ST<:Symmetric{T}} <: AbstractElementCache
         B = zeros(6, 12)
         aux_geometric_Ks = zeros(4, 4)
         E = Symmetric(zeros(3, 3))
-        new{Float64,Symmetric{Float64}}(fint, Ks, σ, ε, F, H, X, J,
-                                        funder, B, aux_geometric_Ks, E)
+        I₃₃ = eye(3)
+        ones₃₃ = ones(3, 3)
+        new{Float64,Symmetric{Float64}}(fint, Ks, S, ∂S∂E, P, ε, F, H, X, J,
+                                        funder, B, aux_geometric_Ks, E, I₃₃, ones₃₃)
     end
 end
 
@@ -167,7 +180,7 @@ end
 and a an element displacement vector `u_e`. This function modifies the cache to avoid memory allocations."
 function internal_forces(m::AbstractHyperElasticMaterial, t::Tetrahedron, u_e::AbstractVector,
                          cache::TetrahedronCache)
-    (; fint, Ks, σ, ε, F, H, X, J, funder, B, aux_geometric_Ks, E) = cache
+    (; fint, Ks, P, ε, F, H, X, J, funder, B, aux_geometric_Ks, E) = cache
 
     # Kinematics
     U = reshape(u_e, 3, 4)
@@ -196,12 +209,12 @@ function internal_forces(m::AbstractHyperElasticMaterial, t::Tetrahedron, u_e::A
     Ks .= Km + Ks
 
     # Piola stress
-    σ .= Symmetric(F * 𝕊)
+    P .= Symmetric(F * 𝕊)
 
     # Right hand Cauchy strain tensor
     ε .= Symmetric(F' * F)
 
-    fint, Ks, σ, ε
+    fint, Ks, P, ε
 end
 
 "
@@ -220,7 +233,7 @@ A 4-tuple containing:
 "
 function internal_forces(m::IsotropicLinearElastic, t::Tetrahedron, u_e::AbstractVector,
                          cache::TetrahedronCache)
-    (; fint, Ks, σ, ε, F, H, X, J, funder, B) = cache
+    (; fint, Ks, S, ∂S∂E, ε, F, H, X, J, funder, B, I₃₃, ones₃₃) = cache
 
     # Kinematics
     ∂X∂ζ = _shape_functions_derivatives(t)
@@ -231,19 +244,19 @@ function internal_forces(m::IsotropicLinearElastic, t::Tetrahedron, u_e::Abstrac
     U = reshape(u_e, 3, 4)
     H .= U * funder'
     ε .= Symmetric(0.5 * (H + H'))
-    F .= eye(3)
+    F .= I₃₃
     _B_mat!(B, funder, F)
 
-    # Stresses
-    σaux, ∂𝕊∂𝔼 = cauchy_stress(m, ε)
-    σ .= Symmetric(σaux)
+    # Stresses (due to stresses are all the same for linear elastic materials cosserat
+    # is used as cache)
+    stress!(S, ∂S∂E, m, ε; cache_ones=ones₃₃, cache_eye=I₃₃)
 
     # Stiffness matrix
-    Ks .= Symmetric(B' * ∂𝕊∂𝔼 * B * vol)
+    Ks .= Symmetric(B' * ∂S∂E * B * vol)
 
     fint .= Ks * u_e
 
-    fint, Ks, σ, ε
+    fint, Ks, S, ε
 end
 
 "Shape function derivatives."

src/Materials/IsotropicLinearElasticMaterial.jl

@@ -7,7 +7,7 @@ using ..LinearElasticMaterials
 using ..Utils
 
 @reexport import ..LinearElasticMaterials: lame_parameters, shear_modulus, poisson_ratio,
-                                           elasticity_modulus, bulk_modulus, cauchy_stress
+                                           elasticity_modulus, bulk_modulus, stress!
 
 export IsotropicLinearElastic
 
@@ -69,13 +69,16 @@ end
 
 "Return the cauchy stress tensor `σ` and the constitutive driver `∂σ∂ϵ`
 considering a `IsotropicLinearElastic` material `m`."
-function cauchy_stress(m::IsotropicLinearElastic, ϵ::AbstractMatrix)
+function stress!(σ::AbstractMatrix{<:Real}, ∂σ∂ϵ::Matrix{<:Real},
+                 m::IsotropicLinearElastic{<:Real}, ϵ::AbstractMatrix{<:Real};
+                 cache_eye::AbstractMatrix{<:Real}=eye(3),
+                 cache_ones::Matrix{<:Real}=ones(3, 3))
     λ, G = lame_parameters(m)
-    σ = Symmetric(λ * tr(ϵ) * eye(3) + 2 * G * ϵ)
 
-    ∂σ∂ϵ = zeros(6, 6)
-    ∂σ∂ϵ[1:3, 1:3] = λ * ones(3, 3) + 2 * G * eye(3)
-    ∂σ∂ϵ[4:6, 4:6] = G * eye(3)
+    σ .= Symmetric(λ * tr(ϵ) * cache_eye + 2 * G * ϵ)
+
+    ∂σ∂ϵ[1:3, 1:3] .= λ * cache_ones + 2 * G * cache_eye
+    ∂σ∂ϵ[4:6, 4:6] .= G * cache_eye
 
     σ, ∂σ∂ϵ
 end

src/Materials/LinearElasticMaterials.jl

@@ -4,7 +4,7 @@ module LinearElasticMaterials
 using ..Materials
 
 export AbstractLinearElasticMaterial, lame_parameters, elasticity_modulus, shear_modulus,
-       bulk_modulus, poisson_ratio, cauchy_stress
+       bulk_modulus, poisson_ratio, stress!
 
 """ Abstract supertype for all elastic material models.
 
@@ -38,8 +38,10 @@ function elasticity_modulus(m::AbstractLinearElasticMaterial) end
 "Return the bulk modulus `K` for an `AbstractLinearElasticMaterial` material `m`."
 function bulk_modulus(m::AbstractLinearElasticMaterial) end
 
-"Return the cauchy stress tensor `σ` and the constitutive driver `∂σ∂ϵ`
+"Return the stress tensor `σ` and the constitutive driver `∂σ∂ϵ`
 considering a `IsotropicLinearElastic` material `m`."
-function cauchy_stress(m::AbstractLinearElasticMaterial, ϵ::AbstractMatrix) end
+function stress!(σ::AbstractMatrix, ∂σ∂ϵ::AbstractMatrix,
+                 m::AbstractLinearElasticMaterial,
+                 ϵ::AbstractMatrix) end
 
 end

test/materials/materials.jl

@@ -66,7 +66,10 @@ mat_label = "steel"
                             σ_vogit[6] σ_vogit[2] σ_vogit[4]
                             σ_vogit[5] σ_vogit[4] σ_vogit[3]])
 
-    σ, ∂σ∂ϵ = cauchy_stress(linear_steel, ϵ)
+    σ = Symmetric(zeros(3, 3))
+    ∂σ∂ϵ = zeros(6, 6)
+
+    σ, ∂σ∂ϵ = stress!(σ, ∂σ∂ϵ, linear_steel, ϵ)
 
     @test σ ≈ σ_expected rtol = RTOL
 end