Merge pull request #219 from kaipartmann/zem-wan

Added the ZEM algorithm of Wan et al. (2019).

Author: kaipartmann

docs/src/public_api_reference.md

@@ -25,6 +25,7 @@ CriticalStretch
 NoCorrection
 EnergySurfaceCorrection
 ZEMSilling
+ZEMWan
 LinearElastic
 NeoHooke
 MooneyRivlin

src/Peridynamics.jl

@@ -13,7 +13,7 @@ export BBMaterial, DHBBMaterial, OSBMaterial, CMaterial, CRMaterial, BACMaterial
        CKIMaterial
 
 # CMaterial related types
-export LinearElastic, NeoHooke, MooneyRivlin, SaintVenantKirchhoff, ZEMSilling
+export LinearElastic, NeoHooke, MooneyRivlin, SaintVenantKirchhoff, ZEMSilling, ZEMWan
 
 # Kernels
 export linear_kernel, cubic_b_spline_kernel

src/discretization/zem_stabilization.jl

@@ -1,6 +1,6 @@
 function get_correction(mat::AbstractBondSystemMaterial{<:AbstractZEMStabilization},
                         ::Int, ::Int, ::Int)
-    return mat.zem_stabilization
+    return mat.zem
 end
 
 """
@@ -19,3 +19,35 @@ struct ZEMSilling <: AbstractZEMStabilization
         return new(Cs)
     end
 end
+
+"""
+    ZEMWan()
+
+Zero-energy mode stabilization algorithm of Wan et al. (2019), which is an improvement to
+Silling's algorithm that does not require a stabilization parameter.
+See also [`CMaterial`](@ref) on how to use this stabilization algorithm.
+"""
+struct ZEMWan <: AbstractZEMStabilization end
+
+function calc_zem_stiffness_tensor(C, Kinv)
+    C_1 = zero(SMatrix{3,3,Float64,9})
+    for i in 1:3, j in 1:3
+        sum_value = 0.0
+        for k in 1:3, l in 1:3
+            @inbounds sum_value += C[i, j, k, l] * Kinv[k, l]
+        end
+        C_1 = setindex(C_1, sum_value, i, j)
+    end
+    return C_1
+end
+
+function calc_rotated_zem_stiffness_tensor!(C_rotated, C, Kinv, R)
+    for m in 1:3, n in 1:3, o in 1:3, p in 1:3
+        sum_value = 0.0
+        for i in 1:3, j in 1:3, k in 1:3, l in 1:3
+            @inbounds sum_value += C[i, j, k, l] * R[i, m] * R[j, n] * R[k, o] * R[l, p]
+        end
+        @inbounds C_rotated[m, n, o, p] = sum_value
+    end
+    return calc_zem_stiffness_tensor(C_rotated, Kinv)
+end

src/physics/constitutive_models.jl

@@ -20,21 +20,43 @@ function first_piola_kirchhoff(::LinearElastic, storage::AbstractStorage,
                                params::AbstractPointParameters, F::SMatrix{3,3,T,9}) where T
     E = 0.5 .* (F' * F - I)
     Evoigt = SVector{6,Float64}(E[1,1], E[2,2], E[3,3], 2 * E[2,3], 2 * E[3,1], 2 * E[1,2])
-    Cvoigt = get_hooke_matrix(params.nu, params.λ, params.μ)
+    Cvoigt = get_hooke_matrix_voigt(params.nu, params.λ, params.μ)
     Pvoigt = Cvoigt * Evoigt
     P = SMatrix{3,3,Float64,9}(Pvoigt[1], Pvoigt[6], Pvoigt[5],
                                Pvoigt[6], Pvoigt[2], Pvoigt[4],
                                Pvoigt[5], Pvoigt[4], Pvoigt[3])
     return P
 end
 
-function get_hooke_matrix(nu, λ, μ)
+function get_hooke_matrix_voigt(nu, λ, μ)
     a = (1 - nu) * λ / nu
-    CVoigt = SMatrix{6,6,Float64,36}(a, λ, λ, 0, 0, 0, λ, a, λ, 0, 0, 0, λ, λ, a, 0, 0, 0,
-                                     0, 0, 0, μ, 0, 0, 0, 0, 0, 0, μ, 0, 0, 0, 0, 0, 0, μ)
-    return CVoigt
+    Cvoigt = SMatrix{6,6,Float64,36}(
+        a, λ, λ, 0, 0, 0,
+        λ, a, λ, 0, 0, 0,
+        λ, λ, a, 0, 0, 0,
+        0, 0, 0, μ, 0, 0,
+        0, 0, 0, 0, μ, 0,
+        0, 0, 0, 0, 0, μ
+    )
+    return Cvoigt
 end
 
+function get_hooke_matrix(nu, λ, μ)
+    Cvoigt = get_hooke_matrix_voigt(nu, λ, μ)
+    C = zero(MArray{NTuple{4,3},Float64,4,81})
+    voigt_map = @SVector [(1,1), (2,2), (3,3), (2,3), (1,3), (1,2)]
+    for I in 1:6, J in 1:6
+        i1, i2 = voigt_map[I]
+        j1, j2 = voigt_map[J]
+        C[i1, i2, j1, j2] = Cvoigt[I, J]
+        C[i2, i1, j1, j2] = Cvoigt[I, J]
+        C[i1, i2, j2, j1] = Cvoigt[I, J]
+        C[i2, i1, j2, j1] = Cvoigt[I, J]
+    end
+    return C
+end
+
+
 @doc raw"""
     NeoHooke
 

src/physics/correspondence.jl

@@ -17,9 +17,11 @@ consistent (correspondence) formulation of non-ordinary state-based peridynamics
     - [`NeoHooke`](@ref)
     - [`MooneyRivlin`](@ref)
     - [`SaintVenantKirchhoff`](@ref)
-- `zem::AbstractZEMStabilization`: Zero-energy mode stabilization. The
-    stabilization algorithm of Silling (2017) is used as default. \\
-    (default: [`ZEMSilling`](@ref))
+- `zem::AbstractZEMStabilization`: Algorithm of zero-energy mode stabilization. \\
+    (default: [`ZEMSilling`](@ref)) \\
+    The following algorithms can be used:
+    - [`ZEMSilling`](@ref)
+    - [`ZEMWan`](@ref)
 - `dmgmodel::AbstractDamageModel`: Damage model defining the damage behavior. \\
     (default: [`CriticalStretch`](@ref))
 - `maxdmg::Float64`: Maximum value of damage a point is allowed to obtain. If this value is
@@ -106,7 +108,7 @@ When specifying the `fields` keyword of [`Job`](@ref) for a [`Body`](@ref) with
 struct CMaterial{CM,ZEM,K,DM} <: AbstractCorrespondenceMaterial{CM,ZEM}
     kernel::K
     constitutive_model::CM
-    zem_stabilization::ZEM
+    zem::ZEM
     dmgmodel::DM
     maxdmg::Float64
     function CMaterial(kernel::K, cm::CM, zem::ZEM, dmgmodel::DM,
@@ -132,15 +134,38 @@ function log_material_property(::Val{:constitutive_model}, mat; indentation)
     return msg_qty("constitutive model", mat.constitutive_model; indentation)
 end
 
-function log_material_property(::Val{:zem_stabilization}, mat; indentation)
-    return msg_qty("zero-energy mode stabilization", mat.zem_stabilization; indentation)
+function log_material_property(::Val{:zem}, mat; indentation)
+    return msg_qty("zero-energy mode stabilization", mat.zem; indentation)
 end
 
 function log_material_property(::Val{:maxdmg}, mat; indentation)
     return msg_qty("maximum damage", mat.maxdmg; indentation)
 end
 
-@params CMaterial StandardPointParameters
+struct CPointParameters <: AbstractPointParameters
+    δ::Float64
+    rho::Float64
+    E::Float64
+    nu::Float64
+    G::Float64
+    K::Float64
+    λ::Float64
+    μ::Float64
+    Gc::Float64
+    εc::Float64
+    bc::Float64
+    C::MArray{NTuple{4,3},Float64,4,81}
+end
+
+function CPointParameters(mat::AbstractMaterial, p::Dict{Symbol,Any})
+    (; δ, rho, E, nu, G, K, λ, μ) = get_required_point_parameters(mat, p)
+    (; Gc, εc) = get_frac_params(mat.dmgmodel, p, δ, K)
+    bc = 18 * K / (π * δ^4) # bond constant
+    C = get_hooke_matrix(nu, λ, μ)
+    return CPointParameters(δ, rho, E, nu, G, K, λ, μ, Gc, εc, bc, C)
+end
+
+@params CMaterial CPointParameters
 
 @storage CMaterial struct CStorage
     @lthfield position::Matrix{Float64}
@@ -187,13 +212,12 @@ function force_density_point!(storage::AbstractStorage, system::AbstractSystem,
     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)
-    zem = mat.zem_stabilization
-    c_force_density!(storage, system, mat, params, zem, PKinv, defgrad_res, 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,
-                                   ::StandardPointParameters, i)
+                                   ::CPointParameters, i)
     (; bonds, volume) = system
     (; bond_active) = storage
     K = zero(SMatrix{3,3,Float64,9})
@@ -217,7 +241,7 @@ function calc_deformation_gradient(storage::CStorage, system::BondSystem, ::CMat
 end
 
 function calc_first_piola_kirchhoff!(storage::CStorage, mat::CMaterial,
-                                     params::StandardPointParameters, defgrad_res, Δt, i)
+                                     params::CPointParameters, defgrad_res, Δt, i)
     (; F, Kinv) = defgrad_res
     P = first_piola_kirchhoff(mat.constitutive_model, storage, params, F)
     PKinv = P * Kinv
@@ -252,6 +276,40 @@ function c_force_density!(storage::AbstractStorage, system::AbstractSystem,
     return nothing
 end
 
+function c_force_density!(storage::AbstractStorage, system::AbstractSystem,
+                          mat::AbstractCorrespondenceMaterial,
+                          params::AbstractPointParameters, zem::ZEMWan, PKinv, defgrad_res,
+                          i)
+    (; bonds, volume) = system
+    (; bond_active) = storage
+    (; F) = defgrad_res
+    C_1 = calc_zem_stiffness_tensor!(storage, system, mat, params, zem, defgrad_res, i)
+    for bond_id in each_bond_idx(system, i)
+        bond = bonds[bond_id]
+        j = bond.neighbor
+        ΔXij = get_vector_diff(system.position, i, j)
+        Δxij = get_vector_diff(storage.position, i, j)
+
+        # improved stabilization from this article:
+        # https://doi.org/10.1007/s10409-019-00873-y
+        ωij = kernel(system, bond_id) * bond_active[bond_id]
+        tzem = ωij * C_1 * (Δxij .- F * ΔXij)
+
+        # update of force density
+        tij = ωij * PKinv * ΔXij + tzem
+        update_add_vector!(storage.b_int, i, tij .* volume[j])
+        update_add_vector!(storage.b_int, j, -tij .* volume[i])
+    end
+    return nothing
+end
+
+function calc_zem_stiffness_tensor!(storage::CStorage, system::BondSystem, mat::CMaterial,
+                                    params::CPointParameters, zem::ZEMWan, defgrad_res, i)
+    (; Kinv) = defgrad_res
+    C_1 = calc_zem_stiffness_tensor(params.C, Kinv)
+    return C_1
+end
+
 function too_much_damage!(storage::AbstractStorage, system::AbstractSystem,
                           mat::AbstractCorrespondenceMaterial, defgrad_res, i)
     (; F) = defgrad_res

src/physics/correspondence_rotated.jl

@@ -14,9 +14,11 @@ consistent (correspondence) formulation of non-ordinary state-based peridynamics
     (default: `LinearElastic()`) \\
     Only the following model can be used:
     - [`LinearElastic`](@ref)
-- `zem::AbstractZEMStabilization`: Zero-energy mode stabilization. The
-    stabilization algorithm of Silling (2017) is used as default. \\
-    (default: [`ZEMSilling`](@ref))
+- `zem::AbstractZEMStabilization`: Algorithm of zero-energy mode stabilization. \\
+    (default: [`ZEMSilling`](@ref)) \\
+    The following algorithms can be used:
+    - [`ZEMSilling`](@ref)
+    - [`ZEMWan`](@ref)
 - `dmgmodel::AbstractDamageModel`: Damage model defining the damage behavior. \\
     (default: [`CriticalStretch`](@ref))
 - `maxdmg::Float64`: Maximum value of damage a point is allowed to obtain. If this value is
@@ -103,7 +105,7 @@ When specifying the `fields` keyword of [`Job`](@ref) for a [`Body`](@ref) with
 struct CRMaterial{CM,ZEM,K,DM} <: AbstractCorrespondenceMaterial{CM,ZEM}
     kernel::K
     constitutive_model::CM
-    zem_stabilization::ZEM
+    zem::ZEM
     dmgmodel::DM
     maxdmg::Float64
     function CRMaterial(kernel::K, cm::CM, zem::ZEM, dmgmodel::DM,
@@ -129,7 +131,7 @@ function Base.show(io::IO, @nospecialize(mat::CRMaterial))
     return nothing
 end
 
-@params CRMaterial StandardPointParameters
+@params CRMaterial CPointParameters
 
 @storage CRMaterial struct CRStorage
     @lthfield position::Matrix{Float64}
@@ -149,6 +151,7 @@ end
     @pointfield von_mises_stress::Vector{Float64}
     @pointfield left_stretch::Matrix{Float64}
     @pointfield rotation::Matrix{Float64}
+    zem_stiffness_rotated::MArray{NTuple{4,3},Float64,4,81}
 end
 
 function init_field(::CRMaterial, ::AbstractTimeSolver, system::BondSystem,
@@ -183,8 +186,13 @@ function init_field(::CRMaterial, ::AbstractTimeSolver, system::BondSystem,
     return R
 end
 
+function init_field(::CRMaterial, ::AbstractTimeSolver, system::BondSystem,
+                    ::Val{:zem_stiffness_rotated})
+    return zero(MArray{NTuple{4,3},Float64,4,81})
+end
+
 function calc_deformation_gradient(storage::CRStorage, system::BondSystem, ::CRMaterial,
-                                   ::StandardPointParameters, i)
+                                   ::CPointParameters, i)
     (; bonds, volume) = system
     (; bond_active) = storage
     K = zero(SMatrix{3,3,Float64,9})
@@ -212,7 +220,7 @@ function calc_deformation_gradient(storage::CRStorage, system::BondSystem, ::CRM
 end
 
 function calc_first_piola_kirchhoff!(storage::CRStorage, mat::CRMaterial,
-                                     params::StandardPointParameters, defgrad_res, Δt, i)
+                                     params::CPointParameters, defgrad_res, Δt, i)
     (; F, Ḟ, Kinv) = defgrad_res
     D = init_stress_rotation!(storage, F, Ḟ, Δt, i)
     if iszero(D)
@@ -231,3 +239,12 @@ function calc_first_piola_kirchhoff!(storage::CRStorage, mat::CRMaterial,
     PKinv = P * Kinv
     return PKinv
 end
+
+function calc_zem_stiffness_tensor!(storage::CRStorage, system::BondSystem, mat::CRMaterial,
+                                    params::CPointParameters, zem::ZEMWan, defgrad_res, i)
+    (; Kinv) = defgrad_res
+    (; rotation, zem_stiffness_rotated) = storage
+    R = get_tensor(rotation, i)
+    C_1 = calc_rotated_zem_stiffness_tensor!(zem_stiffness_rotated, params.C, Kinv, R)
+    return C_1
+end

test/discretization/test_bond_system.jl

@@ -223,6 +223,6 @@ end
     @test contains(msg, "CriticalStretch")
     msg = Peridynamics.log_material_property(Val(:kernel), mat; indentation)
     @test contains(msg, "linear_kernel")
-    msg = Peridynamics.log_material_property(Val(:zem_stabilization), mat; indentation)
+    msg = Peridynamics.log_material_property(Val(:zem), mat; indentation)
     @test contains(msg, "ZEMSilling")
 end