diff --git a/Project.toml b/Project.toml index 346b6869..13559fb5 100644 --- a/Project.toml +++ b/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"} diff --git a/src/Utils/MaterialModelsBase.jl b/src/Utils/MaterialModelsBase.jl index b8fb8648..689c701b 100644 --- a/src/Utils/MaterialModelsBase.jl +++ b/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...) diff --git a/test/example_elements.jl b/test/example_elements.jl index 9b742bc6..70694327 100644 --- a/test/example_elements.jl +++ b/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() diff --git a/test/runtests.jl b/test/runtests.jl index 1ba3cf97..32035017 100644 --- a/test/runtests.jl +++ b/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")