Merge branch 'main' into mpi_test_fixes

Author: kaipartmann

Project.toml

@@ -1,7 +1,7 @@
 name = "Peridynamics"
 uuid = "4dc47793-80f3-4232-b30e-ca78ca9d621b"
 authors = ["Kai Partmann"]
-version = "0.4.0"
+version = "0.4.1"
 
 [deps]
 AbaqusReader = "bc6b9049-e460-56d6-94b4-a597b2c0390d"

docs/src/public_api_reference.md

@@ -11,6 +11,7 @@ Pages = ["public_api_reference.md"]
 ## Material models
 ```@docs
 BBMaterial
+DHBBMaterial
 OSBMaterial
 CMaterial
 CRMaterial
@@ -24,6 +25,7 @@ CriticalStretch
 NoCorrection
 EnergySurfaceCorrection
 ZEMSilling
+ZEMWan
 LinearElastic
 NeoHooke
 MooneyRivlin

src/Peridynamics.jl

@@ -9,10 +9,11 @@ end
 import LibGit2, Dates
 
 # Material models
-export BBMaterial, OSBMaterial, CMaterial, CRMaterial, BACMaterial, CKIMaterial
+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
@@ -75,6 +76,7 @@ abstract type AbstractCorrection end
 abstract type AbstractStorage end
 abstract type AbstractCondition end
 abstract type AbstractBondSystemMaterial{Correction} <: AbstractMaterial end
+abstract type AbstractBondBasedMaterial{CM} <: AbstractBondSystemMaterial{CM} end
 abstract type AbstractCorrespondenceMaterial{CM,ZEM} <: AbstractBondSystemMaterial{ZEM} end
 abstract type AbstractBondAssociatedSystemMaterial <: AbstractMaterial end
 abstract type AbstractConstitutiveModel end
@@ -132,6 +134,7 @@ include("time_solvers/velocity_verlet.jl")
 include("time_solvers/dynamic_relaxation.jl")
 
 include("physics/bond_based.jl")
+include("physics/dh_bond_based.jl")
 include("physics/continuum_kinematics_inspired.jl")
 include("physics/ordinary_state_based.jl")
 include("physics/constitutive_models.jl")

src/core/storages.jl

@@ -57,7 +57,7 @@ Example definition of the storage for the bond-based material:
     bond_active::Vector{Bool}
     @pointfield n_active_bonds::Vector{Int}
 end
-
+````
 """
 macro storage end
 

src/discretization/body.jl

@@ -7,11 +7,15 @@ Construct a `Body` for a peridynamics simulation.
 # Arguments
 - `material::AbstractMaterial`: The material which is defined for the whole body.
     Available material models:
-    - [`BBMaterial`](@ref): Bond-based peridynamics
-    - [`OSBMaterial`](@ref): Ordinary state-based peridynamics
-    - [`CMaterial`](@ref): Correspondence formulation
-    - [`BACMaterial`](@ref): Bond-associated correspondence formulation of Chen and Spencer
-    - [`CKIMaterial`](@ref): Continuum-kinematics-inspired peridynamics
+    - [`BBMaterial`](@ref): Bond-based peridynamics.
+    - [`DHBBMaterial`](@ref): Dual-horizon bond-based peridynamics.
+    - [`OSBMaterial`](@ref): Ordinary state-based peridynamics, also called linear
+        peridynamic solid (LPS).
+    - [`CMaterial`](@ref): Correspondence formulation.
+    - [`CRMaterial`](@ref): Correspondence formulation stress rotation for objectivity
+        enforcement.
+    - [`BACMaterial`](@ref): Bond-associated correspondence formulation of Chen and Spencer.
+    - [`CKIMaterial`](@ref): Continuum-kinematics-inspired peridynamics.
 - `position::AbstractMatrix`: A `3×n` matrix with the point position of the `n` points.
 - `volume::AbstractVector`: A vector with the volume of each point.
 - `inp_file::AbstractString`: An Abaqus input file containing meshes, imported with

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/bond_based.jl

@@ -70,7 +70,7 @@ When specifying the `fields` keyword of [`Job`](@ref) for a [`Body`](@ref) with
 - `damage::Vector{Float64}`: Damage of each point.
 - `n_active_bonds::Vector{Int}`: Number of intact bonds of each point.
 """
-struct BBMaterial{Correction,DM} <: AbstractBondSystemMaterial{Correction}
+struct BBMaterial{Correction,DM} <: AbstractBondBasedMaterial{Correction}
     dmgmodel::DM
     function BBMaterial{C}(dmgmodel::DM) where {C,DM}
         new{C,DM}(dmgmodel)
@@ -83,6 +83,14 @@ end
 BBMaterial(; kwargs...) = BBMaterial{NoCorrection}(; kwargs...)
 
 function StandardPointParameters(mat::BBMaterial, p::Dict{Symbol,Any})
+    (; δ, rho, E, nu, G, K, λ, μ) = get_required_point_parameters_bb(mat, p)
+    (; Gc, εc) = get_frac_params(mat.dmgmodel, p, δ, K)
+    bc = 18 * K / (π * δ^4) # bond constant
+    return StandardPointParameters(δ, rho, E, nu, G, K, λ, μ, Gc, εc, bc)
+end
+
+function get_required_point_parameters_bb(mat::AbstractBondBasedMaterial,
+                                          p::Dict{Symbol,Any})
     par = get_given_elastic_params(p)
     (; E, nu, G, K, λ, μ) = par
     if isfinite(nu) && !isapprox(nu, 0.25)
@@ -100,9 +108,7 @@ function StandardPointParameters(mat::BBMaterial, p::Dict{Symbol,Any})
         msg *= "Please define either only one or two fitting elastic parameters!\n"
         throw(ArgumentError(msg))
     end
-    (; Gc, εc) = get_frac_params(mat.dmgmodel, p, δ, K)
-    bc = 18 * K / (π * δ^4) # bond constant
-    return StandardPointParameters(δ, rho, E, nu, G, K, λ, μ, Gc, εc, bc)
+    return (; δ, rho, E, nu, G, K, λ, μ)
 end
 
 @params BBMaterial StandardPointParameters
@@ -130,8 +136,8 @@ function init_field(::BBMaterial, ::AbstractTimeSolver, system::BondSystem,
 end
 
 # Customized calc_failure to save the bond stretch ε for force density calculation
-function calc_failure!(storage::BBStorage, system::BondSystem,
-                       ::BBMaterial, ::CriticalStretch,
+function calc_failure!(storage::AbstractStorage, system::BondSystem,
+                       ::AbstractBondBasedMaterial, ::CriticalStretch,
                        paramsetup::AbstractParameterSetup, i)
     (; εc) = get_params(paramsetup, i)
     (; position, n_active_bonds, bond_active, bond_stretch) = storage

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