Merge branch 'main' into energy_surface_correction

Author: kaipartmann

src/auxiliary/io.jl

@@ -102,10 +102,10 @@ end
 function check_export_fields(::Type{S}, fields::Vector{Symbol}) where {S}
     allowed_fields = point_data_fields(S)
     for f in fields
-        if !in(f, allowed_fields) && !custom_field(f)
+        if !in(f, allowed_fields) && !custom_field(S, f)
             msg = "unknown point data field `:$(f)` specified for export!\n"
             msg *= "If you intend to export a custom field, please define a method:"
-            msg *= "\n   `custom_field(::Val{:$(f)}) = true`\n"
+            msg *= "\n   `custom_field(::Type{MyStorage}, ::Val{:$(f)}) = true`\n"
             msg *= "Otherwise, see here all available point data fields of $S:\n"
             for allowed_name in allowed_fields
                 msg *= "  - $allowed_name\n"
@@ -116,10 +116,10 @@ function check_export_fields(::Type{S}, fields::Vector{Symbol}) where {S}
     return nothing
 end
 
-custom_field(field::Symbol) = custom_field(Val(field))
+custom_field(S::Type{<:AbstractStorage}, field::Symbol) = custom_field(S, Val(field))
 # this function can be specialized to indicate custom fields
 # if this is not done, the export_field function will error out!
-custom_field(::Val{field}) where {field} = false
+custom_field(::Type{<:AbstractStorage}, ::Val{field}) where {field} = false
 
 function get_vtk_filebase(body::AbstractBody, root::AbstractString)
     body_name = replace(string(get_name(body)), " " => "_")

src/physics/correspondence.jl

@@ -105,6 +105,8 @@ When specifying the `fields` keyword of [`Job`](@ref) for a [`Body`](@ref) with
 - `n_active_bonds::Vector{Int}`: Number of intact bonds of each point.
 - `cauchy_stress::Matrix{Float64}`: Cauchy stress tensor of each point.
 - `von_mises_stress::Vector{Float64}`: Von Mises stress of each point.
+- `hydrostatic_stress::Vector{Float64}`: Hydrostatic stress of each point.
+- `strain_energy_density::Vector{Float64}`: Strain energy density of each point.
 """
 struct CMaterial{CM,ZEM,K,DM} <: AbstractCorrespondenceMaterial{CM,ZEM}
     kernel::K
@@ -182,36 +184,48 @@ end
     @pointfield damage::Vector{Float64}
     bond_active::Vector{Bool}
     @pointfield n_active_bonds::Vector{Int}
+    @pointfield defgrad::Matrix{Float64}
     @pointfield cauchy_stress::Matrix{Float64}
     @pointfield von_mises_stress::Vector{Float64}
+    @pointfield strain_energy_density::Vector{Float64}
 end
 
 function init_field(::CMaterial, ::AbstractTimeSolver, system::BondSystem, ::Val{:b_int})
     return zeros(3, get_n_points(system))
 end
 
-function init_field(::CMaterial, ::AbstractTimeSolver, system::BondSystem,
-                    ::Val{:cauchy_stress})
+function init_field(::AbstractCorrespondenceMaterial, ::AbstractTimeSolver,
+                    system::BondSystem, ::Val{:defgrad})
     return zeros(9, get_n_loc_points(system))
 end
 
-function init_field(::CMaterial, ::AbstractTimeSolver, system::BondSystem,
-                    ::Val{:von_mises_stress})
+function init_field(::AbstractCorrespondenceMaterial, ::AbstractTimeSolver,
+                    system::BondSystem, ::Val{:cauchy_stress})
+    return zeros(9, get_n_loc_points(system))
+end
+
+function init_field(::AbstractCorrespondenceMaterial, ::AbstractTimeSolver,
+                    system::BondSystem, ::Val{:von_mises_stress})
+    return zeros(get_n_loc_points(system))
+end
+
+function init_field(::AbstractCorrespondenceMaterial, ::AbstractTimeSolver,
+                    system::BondSystem, ::Val{:strain_energy_density})
     return zeros(get_n_loc_points(system))
 end
 
 function force_density_point!(storage::AbstractStorage, system::AbstractSystem,
                               mat::AbstractCorrespondenceMaterial,
                               params::AbstractPointParameters, t, Δt, i)
-    defgrad_res = calc_deformation_gradient(storage, system, mat, params, i)
+    defgrad_res = calc_deformation_gradient!(storage, system, mat, params, i)
     too_much_damage!(storage, system, mat, defgrad_res, i) && return nothing
     PKinv = calc_first_piola_kirchhoff!(storage, mat, params, defgrad_res, Δt, i)
     c_force_density!(storage, system, mat, params, mat.zem, PKinv, defgrad_res, i)
     return nothing
 end
 
-function calc_deformation_gradient(storage::CStorage, system::BondSystem, ::CMaterial,
-                                   ::CPointParameters, i)
+function calc_deformation_gradient!(storage::CStorage, system::BondSystem, ::CMaterial,
+                                    ::CPointParameters, i)
     (; bonds, volume) = system
     (; bond_active) = storage
     K = zero(SMatrix{3,3,Float64,9})
@@ -231,6 +245,7 @@ function calc_deformation_gradient(storage::CStorage, system::BondSystem, ::CMat
     end
     Kinv = inv(K)
     F = _F * Kinv
+    Peridynamics.update_tensor!(storage.defgrad, i, F)
     return (; F, Kinv, ω0)
 end
 
@@ -317,11 +332,35 @@ function too_much_damage!(storage::AbstractStorage, system::AbstractSystem,
 end
 
 # calculate the von Mises stress from the Cauchy stress tensor just when exporting
-function export_field(::Val{:von_mises_stress}, ::CMaterial, system::BondSystem,
-                      storage::AbstractStorage, ::AbstractParameterSetup, t)
+function export_field(::Val{:von_mises_stress}, ::CMaterial,
+                      system::BondSystem, storage::AbstractStorage, ::AbstractParameterSetup, t)
     for i in each_point_idx(system)
         σ = get_tensor(storage.cauchy_stress, i)
         storage.von_mises_stress[i] = von_mises_stress(σ)
     end
     return storage.von_mises_stress
 end
+
+# calculate the von hydrostatic stress from the Cauchy stress tensor just when exporting,
+# use the `von_mises_stress` field to store the hydrostatic stress
+function export_field(::Val{:hydrostatic_stress}, ::CMaterial,
+                      system::BondSystem, storage::CStorage, ::AbstractParameterSetup, t)
+    for i in each_point_idx(system)
+        σ = get_tensor(storage.cauchy_stress, i)
+        storage.von_mises_stress[i] = 1/3 * (σ[1,1] + σ[2,2] + σ[3,3])
+    end
+    return storage.von_mises_stress
+end
+custom_field(::Type{CStorage}, ::Val{:hydrostatic_stress}) = true
+
+# calculate the strain energy density from the deformation gradient just when exporting
+function export_field(::Val{:strain_energy_density}, mat::CMaterial, system::BondSystem,
+                      storage::AbstractStorage, paramsetup::AbstractParameterSetup, t)
+    model = mat.constitutive_model
+    for i in each_point_idx(system)
+        params = get_params(paramsetup, i)
+        F = get_tensor(storage.defgrad, i)
+        storage.strain_energy_density[i] = strain_energy_density(model, storage, params, F)
+    end
+    return storage.strain_energy_density
+end

src/physics/correspondence_rotated.jl

@@ -57,7 +57,10 @@ end
     bond_active::Vector{Bool}
     @pointfield n_active_bonds::Vector{Int}
     @pointfield unrotated_stress::Matrix{Float64}
+    @pointfield defgrad::Matrix{Float64}
+    @pointfield cauchy_stress::Matrix{Float64}
     @pointfield von_mises_stress::Vector{Float64}
+    @pointfield strain_energy_density::Vector{Float64}
     @pointfield left_stretch::Matrix{Float64}
     @pointfield rotation::Matrix{Float64}
     zem_stiffness_rotated::MArray{NTuple{4,3},Float64,4,81}
@@ -77,11 +80,6 @@ function init_field(::CRMaterial, ::AbstractTimeSolver, system::BondSystem,
     return zeros(9, get_n_loc_points(system))
 end
 
-function init_field(::CRMaterial, ::AbstractTimeSolver, system::BondSystem,
-                    ::Val{:von_mises_stress})
-    return zeros(get_n_loc_points(system))
-end
-
 function init_field(::CRMaterial, ::AbstractTimeSolver, system::BondSystem,
                     ::Val{:left_stretch})
     V = zeros(9, get_n_loc_points(system))
@@ -101,8 +99,8 @@ function init_field(::CRMaterial, ::AbstractTimeSolver, system::BondSystem,
     return zero(MArray{NTuple{4,3},Float64,4,81})
 end
 
-function calc_deformation_gradient(storage::CRStorage, system::BondSystem, ::CRMaterial,
-                                   ::CPointParameters, i)
+function calc_deformation_gradient!(storage::CRStorage, system::BondSystem, ::CRMaterial,
+                                    ::CPointParameters, i)
     (; bonds, volume) = system
     (; bond_active) = storage
     K = zero(SMatrix{3,3,Float64,9})
@@ -126,6 +124,7 @@ function calc_deformation_gradient(storage::CRStorage, system::BondSystem, ::CRM
     Kinv = inv(K)
     F = _F * Kinv
     Ḟ = _Ḟ * Kinv
+    Peridynamics.update_tensor!(storage.defgrad, i, F)
     return (; F, Ḟ, Kinv, ω0)
 end
 
@@ -158,6 +157,17 @@ function calc_zem_stiffness_tensor!(storage::CRStorage, system::BondSystem, mat:
     return C_1
 end
 
+# calculate the von Mises stress from the Cauchy stress tensor just when exporting
+function export_field(::Val{:cauchy_stress}, ::CRMaterial, system::BondSystem,
+                      storage::AbstractStorage, ::AbstractParameterSetup, t)
+    for i in each_point_idx(system)
+        σ = get_tensor(storage.unrotated_stress, i)
+        T = rotate_stress(storage::AbstractStorage, σ, i)
+        Peridynamics.update_tensor!(storage.cauchy_stress, i, T)
+    end
+    return storage.cauchy_stress
+end
+
 # calculate the von Mises stress from the Cauchy stress tensor just when exporting
 function export_field(::Val{:von_mises_stress}, ::CRMaterial, system::BondSystem,
                       storage::AbstractStorage, ::AbstractParameterSetup, t)
@@ -168,3 +178,28 @@ function export_field(::Val{:von_mises_stress}, ::CRMaterial, system::BondSystem
     end
     return storage.von_mises_stress
 end
+
+# calculate the von hydrostatic stress from the Cauchy stress tensor just when exporting,
+# use the `von_mises_stress` field to store the hydrostatic stress
+function export_field(::Val{:hydrostatic_stress}, ::CRMaterial,
+                      system::BondSystem, storage::CRStorage, ::AbstractParameterSetup, t)
+    for i in each_point_idx(system)
+        σ = get_tensor(storage.unrotated_stress, i)
+        T = rotate_stress(storage::AbstractStorage, σ, i)
+        storage.von_mises_stress[i] = 1/3 * (T[1,1] + T[2,2] + T[3,3])
+    end
+    return storage.von_mises_stress
+end
+custom_field(::Type{CRStorage}, ::Val{:hydrostatic_stress}) = true
+
+# calculate the strain energy density from the deformation gradient just when exporting
+function export_field(::Val{:strain_energy_density}, mat::CRMaterial, system::BondSystem,
+                      storage::AbstractStorage, paramsetup::AbstractParameterSetup, t)
+    model = mat.constitutive_model
+    for i in each_point_idx(system)
+        params = get_params(paramsetup, i)
+        F = get_tensor(storage.defgrad, i)
+        storage.strain_energy_density[i] = strain_energy_density(model, storage, params, F)
+    end
+    return storage.strain_energy_density
+end

src/physics/rk_correspondence.jl

@@ -113,8 +113,10 @@ When specifying the `fields` keyword of [`Job`](@ref) for a [`Body`](@ref) with
 - `b_ext::Matrix{Float64}`: External force density of each point
 - `damage::Vector{Float64}`: Damage of each point
 - `n_active_bonds::Vector{Int}`: Number of intact bonds of each point
-- `stress::Matrix{Float64}`: Stress tensor of each point
+- `cauchy_stress::Matrix{Float64}`: Cauchy stress tensor of each point
 - `von_mises_stress::Vector{Float64}`: Von Mises stress of each point
+- `hydrostatic_stress::Vector{Float64}`: Hydrostatic stress of each point.
+- `strain_energy_density::Vector{Float64}`: Strain energy density of each point.
 """
 struct RKCMaterial{CM,K,DM} <: AbstractRKCMaterial{CM,NoCorrection}
     kernel::K
@@ -172,12 +174,13 @@ end
     bond_active::Vector{Bool}
     @pointfield n_active_bonds::Vector{Int}
     @pointfield update_gradients::Vector{Bool}
-    @pointfield stress::Matrix{Float64}
+    @pointfield cauchy_stress::Matrix{Float64}
     @pointfield von_mises_stress::Vector{Float64}
+    @pointfield strain_energy_density::Vector{Float64}
     @lthfield defgrad::Matrix{Float64}
     @lthfield weighted_volume::Vector{Float64}
     gradient_weight::Matrix{Float64}
-    first_piola_kirchhoff::Matrix{Float64}
+    bond_first_piola_kirchhoff::Matrix{Float64}
 end
 
 function init_field(::AbstractRKCMaterial, ::AbstractTimeSolver, system::BondSystem,
@@ -191,7 +194,7 @@ function init_field(::AbstractRKCMaterial, ::AbstractTimeSolver, system::BondSys
 end
 
 function init_field(::AbstractRKCMaterial, ::AbstractTimeSolver, system::BondSystem,
-                    ::Val{:stress})
+                    ::Val{:cauchy_stress})
     return zeros(9, get_n_loc_points(system))
 end
 
@@ -200,6 +203,11 @@ function init_field(::AbstractRKCMaterial, ::AbstractTimeSolver, system::BondSys
     return zeros(get_n_loc_points(system))
 end
 
+function init_field(::AbstractRKCMaterial, ::AbstractTimeSolver, system::BondSystem,
+                    ::Val{:strain_energy_density})
+    return zeros(get_n_loc_points(system))
+end
+
 function init_field(::AbstractRKCMaterial, ::AbstractTimeSolver, system::BondSystem,
                     ::Val{:defgrad})
     return zeros(9, get_n_points(system))
@@ -216,7 +224,7 @@ function init_field(::AbstractRKCMaterial, ::AbstractTimeSolver, system::BondSys
 end
 
 function init_field(::AbstractRKCMaterial, ::AbstractTimeSolver, system::BondSystem,
-                    ::Val{:first_piola_kirchhoff})
+                    ::Val{:bond_first_piola_kirchhoff})
     return zeros(9, get_n_bonds(system))
 end
 
@@ -499,15 +507,15 @@ function rkc_force_density!(storage::AbstractStorage, system::AbstractBondSystem
                             mat::AbstractRKCMaterial, params::AbstractPointParameters,
                             ∑P, t, Δt, i)
     (; bonds, volume) = system
-    (; bond_active, gradient_weight, first_piola_kirchhoff, weighted_volume,
+    (; bond_active, gradient_weight, bond_first_piola_kirchhoff, weighted_volume,
        b_int) = storage
     wi = weighted_volume[i]
     for bond_id in each_bond_idx(system, i)
         if bond_active[bond_id]
             bond = bonds[bond_id]
             j, L = bond.neighbor, bond.length
             ΔXij = get_vector_diff(system.position, i, j)
-            Pij = get_tensor(first_piola_kirchhoff, bond_id)
+            Pij = get_tensor(bond_first_piola_kirchhoff, bond_id)
             Φij = get_vector(gradient_weight, bond_id)
             ϕ = 1 / wi
             ω̃ij = kernel(system, bond_id) * ϕ
@@ -522,7 +530,7 @@ end
 function calc_first_piola_kirchhoff!(storage::RKCStorage, mat::RKCMaterial,
                                      params::StandardPointParameters, F, bond_id)
     P = first_piola_kirchhoff(mat.constitutive_model, storage, params, F)
-    update_tensor!(storage.first_piola_kirchhoff, bond_id, P)
+    update_tensor!(storage.bond_first_piola_kirchhoff, bond_id, P)
     return P
 end
 
@@ -542,3 +550,77 @@ function bond_avg(Fi, Fj, ΔXij, Δxij, L)
     Fij = Favg + Fcor
     return Fij
 end
+
+function cauchy_stress_point!(storage::AbstractStorage, system::BondSystem,
+                              ::AbstractRKCMaterial, ::StandardPointParameters, i)
+    (; bonds) = system
+    (; bond_active, defgrad, bond_first_piola_kirchhoff, n_active_bonds) = storage
+    Fi = get_tensor(defgrad, i)
+    σi = zero(SMatrix{3,3,Float64,9})
+    for bond_id in each_bond_idx(system, i)
+        if bond_active[bond_id]
+            bond = bonds[bond_id]
+            j, L = bond.neighbor, bond.length
+            ΔXij = get_vector_diff(system.position, i, j)
+            Δxij = get_vector_diff(storage.position, i, j)
+            Fj = get_tensor(defgrad, j)
+            Fij = bond_avg(Fi, Fj, ΔXij, Δxij, L)
+            Pij = get_tensor(bond_first_piola_kirchhoff, bond_id)
+            σij = cauchy_stress(Pij, Fij)
+            σi += σij
+        end
+    end
+    update_tensor!(storage.cauchy_stress, i, σi ./ n_active_bonds[i])
+    return nothing
+end
+
+# calculate the Cauchy stress tensor just when exporting
+function export_field(::Val{:cauchy_stress}, mat::AbstractRKCMaterial, system::BondSystem,
+                      storage::AbstractStorage, paramsetup::AbstractParameterSetup, t)
+    for i in each_point_idx(system)
+        params = get_params(paramsetup, i)
+        cauchy_stress_point!(storage, system, mat, params, i)
+    end
+    return storage.cauchy_stress
+end
+
+# calculate the von Mises stress from the Cauchy stress tensor just when exporting
+function export_field(::Val{:von_mises_stress}, mat::AbstractRKCMaterial,
+                      system::BondSystem, storage::AbstractStorage,
+                      paramsetup::AbstractParameterSetup, t)
+    for i in each_point_idx(system)
+        params = get_params(paramsetup, i)
+        cauchy_stress_point!(storage, system, mat, params, i)
+        σ = get_tensor(storage.cauchy_stress, i)
+        storage.von_mises_stress[i] = von_mises_stress(σ)
+    end
+    return storage.von_mises_stress
+end
+
+# calculate the von hydrostatic stress from the Cauchy stress tensor just when exporting,
+# use the `von_mises_stress` field to store the hydrostatic stress
+function export_field(::Val{:hydrostatic_stress}, mat::AbstractRKCMaterial,
+                      system::BondSystem, storage::AbstractStorage,
+                      paramsetup::AbstractParameterSetup, t)
+    for i in each_point_idx(system)
+        params = get_params(paramsetup, i)
+        cauchy_stress_point!(storage, system, mat, params, i)
+        σ = get_tensor(storage.cauchy_stress, i)
+        storage.von_mises_stress[i] = 1/3 * (σ[1,1] + σ[2,2] + σ[3,3])
+    end
+    return storage.von_mises_stress
+end
+custom_field(::Type{RKCStorage}, ::Val{:hydrostatic_stress}) = true
+
+# calculate the strain energy density from the deformation gradient just when exporting
+function export_field(::Val{:strain_energy_density}, mat::AbstractRKCMaterial,
+                      system::BondSystem, storage::AbstractStorage,
+                      paramsetup::AbstractParameterSetup, t)
+    model = mat.constitutive_model
+    for i in each_point_idx(system)
+        params = get_params(paramsetup, i)
+        F = get_tensor(storage.defgrad, i)
+        storage.strain_energy_density[i] = strain_energy_density(model, storage, params, F)
+    end
+    return storage.strain_energy_density
+end