Add tests for, and fix, finite strain stress iterations

Author: KnutAM

src/stressiterations.jl

@@ -234,11 +234,11 @@ const State1D = Union{UniaxialStrain, UniaxialStress}
 # Translate from reduced dim input to full tensors
 expand_tensordim(::AbstractStressState, a::Tensors.AllTensors{3}) = a
 #expand_tensordim(::AbstractStressState, v::Vec{dim,T}) where {dim,T} = Vec{3,T}(i->i>dim ? zero(T) : v[i])
-function expand_tensordim(::AbstractStressState, a::Tensor{2,dim,T}) where {dim,T} 
-    return Tensor{2,3}((i,j)-> (i<=dim && j<=dim) ? a[i,j] : zero(T))
+function expand_tensordim(::AbstractStressState, F::Tensor{2,dim,T}) where {dim,T} 
+    return Tensor{2,3}((i,j)-> (i<=dim && j<=dim) ? F[i,j] : (i == j ? one(T) : zero(T)))
 end
-function expand_tensordim(::AbstractStressState, a::SymmetricTensor{2,dim,T}) where {dim,T} 
-    return SymmetricTensor{2,3}((i,j)-> (i<=dim && j<=dim) ? a[i,j] : zero(T))
+function expand_tensordim(::AbstractStressState, ϵ::SymmetricTensor{2,dim,T}) where {dim,T} 
+    return SymmetricTensor{2,3}((i,j)-> (i<=dim && j<=dim) ? ϵ[i,j] : zero(T))
 end
 
 # Translate from full tensors to reduced output

src/vector_conversion.jl

@@ -4,10 +4,10 @@
 
 Return the type of the primary input (strain-like) and associated output (stress-like) 
 to the material model. The default is `SymmetricTensor{2, 3}` for small-strain material 
-models in 3D, but typically, it can be 
-* Small strain material model: `SymmetricTensor{2, dim}`
-* Large strain material model: `Tensor{2, dim}`
-* Traction-separation law: `Vec{dim}`
+models in 3D, but it can be 
+* Small strain material model: Small strain tensor, `ϵ::SymmetricTensor{2, dim}`
+* Large strain material model: Deformation gradient, `F::Tensor{2, dim}`
+* Traction-separation law: Displacement jump, `Δu::Vec{dim}`
 
 Where `dim = 3` is most common, but any value, 1, 2, or 3, is valid. 
 """

test/TestMaterials.jl

@@ -6,8 +6,7 @@ using Tensors, StaticArrays, ForwardDiff
 import MaterialModelsBase as MMB
 
 # Exported materials
-export LinearElastic, ViscoElastic
-# Todo: Also NeoHookean to test nonlinear geometric materials
+#export LinearElastic, ViscoElastic
 
 # LinearElastic material 
 struct LinearElastic{T} <: AbstractMaterial 

test/runtests.jl

@@ -4,8 +4,42 @@ using Tensors, StaticArrays
 
 const MMB = MaterialModelsBase
 
+# <temporary, move to TestMaterials>
+@kwdef struct StVenant{T} <: AbstractMaterial
+    G::T
+    K::T
+end
+function MMB.material_response(m::StVenant, F::Tensor{2,3}, old, args...)
+    function free_energy(defgrad)
+        E = (tdot(defgrad) - one(defgrad)) / 2
+        Edev = dev(E)
+        return m.G * Edev ⊡ Edev + m.K * tr(E)^2 / 2
+    end
+    ∂P∂F, P, _ = hessian(free_energy, F, :all)
+    return P, ∂P∂F, old
+end
+
+@kwdef struct NeoHooke{T} <: AbstractMaterial
+    G::T
+    K::T
+end
+function MMB.material_response(m::NeoHooke, F::Tensor{2,3}, old, args...)
+    function free_energy(defgrad)
+        C = tdot(defgrad)
+        detC = det(C)
+        ΨG = (m.G / 2) * (tr(C) / cbrt(detC) - 3)
+        ΨK = (m.K / 2) * (sqrt(detC) - 1)^2
+        return ΨG + ΨK
+    end
+    ∂P∂F, P, _ = hessian(free_energy, F, :all)
+    return P, ∂P∂F, old
+end
+
+MMB.get_tensorbase(::Union{StVenant, NeoHooke}) = Tensor{2,3}
+# </temporary>
+
 include("TestMaterials.jl")
-using .TestMaterials
+using .TestMaterials: LinearElastic, ViscoElastic
 include("utils4testing.jl")
 
 include("vector_conversion.jl")

test/stressiterations.jl

@@ -16,7 +16,11 @@ iter_mandel = ( ([2,3,4,5,6,7,8,9],[2,3,4,5,6]),
             ared = rand(TT{2,_getdim(stress_state)})
             afull = MMB.expand_tensordim(stress_state, ared)
             @test _getdim(afull) == 3    # Full dimensional tensor
-            @test norm(ared) ≈ norm(afull)    # Only add zeros
+            if TT == SymmetricTensor
+                @test norm(ared) ≈ norm(afull)    # Only add zeros
+            else # Add ones to diagonal
+                @test sqrt(norm(ared)^2 + (3 - _getdim(ared))) ≈ norm(afull)
+            end
             # Test invertability of conversions
             @test ared ≈ MMB.reduce_tensordim(stress_state, afull)
         end
@@ -137,6 +141,41 @@ end
     @test ϵfull_w == ϵfull
 end
 
+@testset "hyperelastic" begin
+    G = 80.e3           # Shear modulus (μ)
+    K = 160.e3          # Bulk modulus
+    E = 9*K*G/(3*K+G)   # Young's modulus
+    ν = E/2G - 1        # Poisson's ratio
+    λ = K - 2G/3        # Lame parameter
+
+    Δϵ = 1e-6   # Small strain to compare with linear case
+    rtol = 1e-5 # Relative tolerance to compare with linear case
+    @testset for m in (StVenant(;G, K), NeoHooke(;G, K))
+        old = initial_material_state(m)
+        # UniaxialStress
+        P, dPdF, state, Ffull = material_response(UniaxialStress(), m, Tensor{2,1}((1 + Δϵ,)), old, 0.0)
+        @test isapprox(P[1,1], E*Δϵ; rtol)
+        @test isapprox(dPdF[1,1,1,1], E; rtol)
+        @test Ffull[2,2] ≈ Ffull[3,3]
+        @test Ffull[2,2] ≈ 1-ν*Δϵ
+        
+        # UniaxialStrain
+        P, dPdF, state, Ffull = material_response(UniaxialStrain(), m, Tensor{2,1}((1 + Δϵ,)), old, 0.0)
+        @test isapprox(P[1,1], (2G+λ)*Δϵ; rtol)
+        @test isapprox(dPdF[1,1,1,1], (2G+λ); rtol)
+        P_full, dPdF_full, _ = material_response(m, Ffull, old, 0.0)
+        @test isapprox(P_full[2,2], λ*Δϵ; rtol)
+
+        # PlaneStress
+        ϵ = Δϵ * rand(SymmetricTensor{2,2})
+        F = one(Tensor{2,2}) + ϵ
+        P, dPdF, state, Ffull = material_response(PlaneStress(), m, F, old, 0.0)
+        Dvoigt = (E/(1-ν^2))*[1 ν 0; ν 1 0; 0 0 (1-ν)/2]
+        @test isapprox(tovoigt(symmetric(dPdF)), Dvoigt; rtol)
+        @test isapprox(P, fromvoigt(SymmetricTensor{4,2}, Dvoigt)⊡ϵ; rtol)
+    end
+end
+
 @testset "visco_elastic" begin
     # Test a nonlinear, state and rate dependent, material by running stress iterations.
     # Check that state variables and time dependence are handled correctly