Add support for finite strain material models (#59)

* Add support for finite strain material models

* Add tests

* Add Newton as specific test-dependency with sources

Author: KnutAM

Project.toml

@@ -20,9 +20,13 @@ julia = "1.11"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+MechanicalMaterialModels = "b3282f9b-607f-4337-ab95-e5488ab5652c"
+Newton = "83aa5b51-0588-403c-85e4-434ec185aae7"
 
 [targets]
-test = ["Test", "Logging", "SparseArrays"]
+test = ["Test", "Logging", "SparseArrays", "MechanicalMaterialModels", "Newton"]
 
 [sources]
 MaterialModelsBase = {url = "https://github.com/KnutAM/MaterialModelsBase.jl.git"}
+MechanicalMaterialModels = {url = "https://github.com/KnutAM/MechanicalMaterialModels.jl.git"}
+Newton = {url = "https://github.com/KnutAM/Newton.jl.git"}

src/Utils/MaterialModelsBase.jl

@@ -16,6 +16,12 @@ where ``\\sigma`` is calculated with the `material_response` function from
 Note that `create_cell_state` is already implemented for `<:AbstractMaterial`. 
 """
 function FerriteAssembly.element_routine!(
+    Ke, re, state::Vector{<:MMB.AbstractMaterialState},
+    ae, material::MMB.AbstractMaterial, cellvalues::AbstractCellValues, buffer)
+    return mechanical_element_routine!(MMB.get_tensorbase(material), Ke, re, state, ae, material, cellvalues, buffer)
+end
+
+function mechanical_element_routine!(::Type{<:SymmetricTensor{2}},
     Ke, re, state::Vector{<:MMB.AbstractMaterialState},
     ae, material::MMB.AbstractMaterial, cellvalues::AbstractCellValues, buffer)
     cache = FerriteAssembly.get_user_cache(buffer)
@@ -40,6 +46,31 @@ function FerriteAssembly.element_routine!(
     end
 end
 
+function mechanical_element_routine!(::Type{<:Tensor{2}},
+    Ke, re, state::Vector{<:MMB.AbstractMaterialState},
+    ae, material::MMB.AbstractMaterial, cellvalues::AbstractCellValues, buffer)
+    cache = FerriteAssembly.get_user_cache(buffer)
+    Δt = FerriteAssembly.get_time_increment(buffer)
+    state_old = FerriteAssembly.get_old_state(buffer)
+    n_basefuncs = getnbasefunctions(cellvalues)
+    for q_point in 1:getnquadpoints(cellvalues)
+        # For each integration point, compute stress and material stiffness
+        H = function_gradient(cellvalues, q_point, ae) # Displacement gradient
+        F = H + one(H) # Deformation gradient
+        P, D, state[q_point] = MMB.material_response(material, F, state_old[q_point], Δt, cache)
+        dΩ = getdetJdV(cellvalues, q_point)
+        for i in 1:n_basefuncs
+            ∇δN = shape_gradient(cellvalues, q_point, i)
+            re[i] += (∇δN ⊡ P) * dΩ # add internal force to residual
+            ∇δN_D = ∇δN ⊡ D         # temporary value for speed
+            for j in 1:n_basefuncs
+                ∇N = shape_gradient(cellvalues, q_point, j)
+                Ke[i, j] += (∇δN_D ⊡ ∇N) * dΩ
+            end
+        end
+    end
+end
+
 """
     FerriteAssembly.element_residual!(
         re, state::Vector{<:MMB.AbstractMaterialState}, ae, 
@@ -49,6 +80,12 @@ The `element_residual!` implementation corresponding to the `element_routine!` i
 for a `MaterialModelsBase.AbstractMaterial`
 """
 function FerriteAssembly.element_residual!(
+    re, state::Vector{<:MMB.AbstractMaterialState},
+    ae, material::MMB.AbstractMaterial, cellvalues::AbstractCellValues, buffer)
+    return mechanical_element_residual!(MMB.get_tensorbase(material), re, state, ae, material, cellvalues, buffer)
+end
+
+function mechanical_element_residual!(::Type{<:SymmetricTensor{2}},
     re, state::Vector{<:MMB.AbstractMaterialState},
     ae, material::MMB.AbstractMaterial, cellvalues::AbstractCellValues, buffer)
     cache = FerriteAssembly.get_user_cache(buffer)
@@ -67,6 +104,26 @@ function FerriteAssembly.element_residual!(
     end
 end
 
+function mechanical_element_residual!(::Type{<:Tensor{2}},
+    re, state::Vector{<:MMB.AbstractMaterialState},
+    ae, material::MMB.AbstractMaterial, cellvalues::AbstractCellValues, buffer)
+    cache = FerriteAssembly.get_user_cache(buffer)
+    Δt = FerriteAssembly.get_time_increment(buffer)
+    state_old = FerriteAssembly.get_old_state(buffer)
+    n_basefuncs = getnbasefunctions(cellvalues)
+    for q_point in 1:getnquadpoints(cellvalues)
+        # For each integration point, compute stress and material stiffness
+        H = function_gradient(cellvalues, q_point, ae) # Displacement gradient
+        F = H + one(H) # Deformation gradient
+        P, _, state[q_point] = MMB.material_response(material, F, state_old[q_point], Δt, cache)
+        dΩ = getdetJdV(cellvalues, q_point)
+        for i in 1:n_basefuncs
+            ∇δN = shape_gradient(cellvalues, q_point, i)
+            re[i] += (∇δN ⊡ P) * dΩ # add internal force to residual
+        end
+    end
+end
+
 """
     FerriteAssembly.create_cell_state(m::MMB.AbstractMaterial, cv::AbstractCellValues, args...)
 

test/example_elements.jl

@@ -73,6 +73,25 @@ end
         test_equality(m, weak, Val(true))
         test_equality(m, mmb, Val(true))
     end
+    @testset "Hyperelasticity" begin
+        G, K = (1 + rand(), 1 + rand())
+        m_3d = MMM.CompressibleNeoHooke(; G, K)
+        m = MMB.ReducedStressState(MMB.PlaneStrain(), m_3d)
+        function weakform_neohooke(δu, ∇δu::Tensor{2,2}, u, ∇u::Tensor{2,2}, u_dot, ∇u_dot)
+            # 2D, consider as plane strain]
+            make3d(t::Tensor{2,2,T}) where {T} = Tensor{2,3,T}((i,j) -> i ≤ 2 && j ≤ 2 ? t[i,j] : zero(T))
+            return weakform_neohooke(δu, make3d(∇δu), u, make3d(∇u), u_dot, ∇u_dot)
+        end
+        function weakform_neohooke(δu, ∇δu::Tensor{2,3}, u, ∇u::Tensor{2,3}, u_dot, ∇u_dot)
+            F = one(∇u) + ∇u
+            C = tdot(F)
+            S = MMM.compute_stress(m_3d, C)
+            P = F ⋅ S
+            return ∇δu ⊡ P
+        end
+        weak = EE.WeakForm(weakform_neohooke)
+        test_equality(m, weak, Val(true))
+    end
     @testset "J2Plasticity" begin
         E = 1.0 + rand()
         H = 0.1 + 0.9*rand()

test/runtests.jl

@@ -5,6 +5,7 @@ using Test
 import FerriteAssembly as FA
 import FerriteAssembly.ExampleElements as EE
 import MaterialModelsBase as MMB
+import MechanicalMaterialModels as MMM
 using Logging
 
 include("replacements.jl")