Cleanup + docs

Project.toml

@@ -8,7 +8,6 @@ AbaqusReader = "bc6b9049-e460-56d6-94b4-a597b2c0390d"
 Base64 = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
 CodecZlib = "944b1d66-785c-5afd-91f1-9de20f533193"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
-Infiltrator = "5903a43b-9cc3-4c30-8d17-598619ec4e9b"
 LibGit2 = "76f85450-5226-5b5a-8eaa-529ad045b433"
 LightXML = "9c8b4983-aa76-5018-a973-4c85ecc9e179"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -27,7 +26,6 @@ AbaqusReader = "0.2"
 Base64 = "1.8"
 CodecZlib = "0.7"
 Dates = "1.8"
-Infiltrator = "1.8.3"
 LibGit2 = "1.8"
 LightXML = "0.9"
 LinearAlgebra = "1.8"

src/Peridynamics.jl

@@ -6,13 +6,10 @@ using Base.Threads, Printf, LinearAlgebra, StaticArrays, PointNeighbors, Progres
     using ThreadPinning
 end
 
-using Infiltrator
-
 import LibGit2, Dates
 
 # Material models
-export BBMaterial, OSBMaterial, CMaterial, CRMaterial, BACMaterial, RKCMaterial,
-       NSCMaterial, NSCRMaterial, CKIMaterial
+export BBMaterial, OSBMaterial, CMaterial, BACMaterial, CKIMaterial
 
 # CMaterial related types
 export NeoHookeNonlinear, SaintVenantKirchhoff, ZEMSilling
@@ -131,11 +128,7 @@ include("physics/continuum_kinematics_inspired.jl")
 include("physics/ordinary_state_based.jl")
 include("physics/constitutive_models.jl")
 include("physics/correspondence.jl")
-include("physics/correspondence_rotated.jl")
-include("physics/ns_correspondence.jl")
-include("physics/ns_correspondence_rotated.jl")
 include("physics/ba_correspondence.jl")
-include("physics/rk_correspondence.jl")
 
 include("VtkReader/VtkReader.jl")
 using .VtkReader

src/discretization/bond_associated_system.jl

@@ -105,34 +105,6 @@ function calc_hood_volumes!(chunk::AbstractBodyChunk{<:BondAssociatedSystem})
     return nothing
 end
 
-# function calc_hood_volumes!(chunk::AbstractBodyChunk{<:BondAssociatedSystem})
-#     (; mat, paramsetup, system, storage) = chunk
-#     (; volume, bonds, hood_volume, ba_hood_volume) = system
-#     (; bond_active) = storage
-
-#     for i in each_point_idx(chunk)
-#         _hood_volume = 0.0 # volume[i]
-#         params = get_params(paramsetup, i)
-#         for bond_idx in each_bond_idx(system, i)
-#             bond = bonds[bond_idx]
-#             j, L = bond.neighbor, bond.length
-#             ωij = influence_function(mat, params, L) * bond_active[bond_idx]
-#             _hood_volume += ωij * volume[j]
-#             _ba_hood_volume = 0.0
-#             for i_bond_idx in each_intersecting_bond_idx(system, i, bond_idx)
-#                 i_bond = bonds[i_bond_idx]
-#                 jj, LL = i_bond.neighbor, i_bond.length
-#                 ωijj = influence_function(mat, params, LL) * bond_active[i_bond_idx]
-#                 _ba_hood_volume += ωijj * volume[jj]
-#             end
-#             ba_hood_volume[bond_idx] = _ba_hood_volume
-#         end
-#         hood_volume[i] = _hood_volume
-#     end
-
-#     return nothing
-# end
-
 @inline get_hood_volume(chunk::AbstractBodyChunk) = chunk.system.hood_volume
 
 function initialize!(dh::AbstractThreadsBodyDataHandler{<:BondAssociatedSystem},

src/discretization/zem_stabilization.jl

@@ -4,9 +4,14 @@ function get_correction(mat::AbstractBondSystemMaterial{<:AbstractZEMStabilizati
 end
 
 """
-    ZEMSilling(; Cs=0.0)
+    ZEMSilling(; Cs)
 
-TODO
+Zero-energy mode stabilization algorithm of Silling (2017). This is necessary for the
+correspondence formulation to stabilize the zero-energy modes.
+See also [`CMaterial`](@ref) on how to use this stabilization algorithm.
+
+# Keywords
+- `Cs::Real`: Stabilization factor. (default: `100.0`)
 """
 struct ZEMSilling <: AbstractZEMStabilization
     Cs::Float64

src/physics/ba_correspondence.jl

@@ -1,47 +1,46 @@
 """
-    CMaterial(; maxdmg, zem)
+    BACMaterial(; kernel, model, maxdmg)
 
-A material type used to assign the material of a [`Body`](@ref) with the local continuum
-consistent (correspondence) formulation of non-ordinary state-based peridynamics.
+A material type used to assign the material of a [`Body`](@ref) with the bond-associated
+correspondence formulation of Chen and Spencer (2019).
 
 # Keywords
+- `kernel::Function`: Kernel function used for weighting the interactions between points.
+    (default: `linear_kernel`)
+- `model::AbstractConstitutiveModel`: Constitutive model defining the material behavior.
+    (default: `LinearElastic()`)
 - `maxdmg::Float64`: Maximum value of damage a point is allowed to obtain. If this value is
     exceeded, all bonds of that point are broken because the deformation gradient would then
     possibly contain `NaN` values.
-    (default: `0.95`)
-- `zem::AbstractZEMStabilization`: zero-energy mode stabilization. The
-    stabilization algorithm of Silling (2017) is used as default.
-    (default: `ZEMSilling()`)
-
-!!! note "Stability of fracture simulations"
-    This formulation is known to be not suitable for fracture simultations without
-    stabilization of the zero-energy modes. Therefore be careful when doing fracture
-    simulations and try out different parameters for `maxdmg` and `zem`.
+    (default: `0.85`)
 
 # Examples
 
 ```julia-repl
-julia> mat = CMaterial()
+julia> mat = BACMaterial()
 CMaterial(maxdmg=0.95, zem_fac=ZEMSilling())
 ```
 
 ---
 
 ```julia
-CMaterial
+BACMaterial{CM,K}
 ```
 
-Material type for the local continuum consistent (correspondence) formulation of
-non-ordinary state-based peridynamics.
+Material type for the bond-associated correspondence formulation of Chen and Spencer (2019).
+
+# Type Parameters
+- `CM`: A constitutive model type. See the constructor docs for more informations.
+- `K`: A kernel function type. See the constructor docs for more informations.
 
 # Fields
+- `kernel::Function`: Kernel function used for weighting the interactions between points.
+- `model::AbstractConstitutiveModel`: Constitutive model defining the material behavior.
 - `maxdmg::Float64`: Maximum value of damage a point is allowed to obtain. See the
     constructor docs for more informations.
-- `zem_fac::Float64`: Correction factor used for zero-energy mode stabilization. See the
-    constructor docs for more informations.
 
 # Allowed material parameters
-When using [`material!`](@ref) on a [`Body`](@ref) with `CMaterial`, then the following
+When using [`material!`](@ref) on a [`Body`](@ref) with `BACMaterial`, then the following
 parameters are allowed:
 - `horizon::Float64`: Radius of point interactions
 - `rho::Float64`: Density
@@ -52,7 +51,7 @@ parameters are allowed:
 
 # Allowed export fields
 When specifying the `fields` keyword of [`Job`](@ref) for a [`Body`](@ref) with
-`CMaterial`, the following fields are allowed:
+`BACMaterial`, the following fields are allowed:
 - `position::Matrix{Float64}`: Position of each point
 - `displacement::Matrix{Float64}`: Displacement of each point
 - `velocity::Matrix{Float64}`: Velocity of each point
@@ -79,7 +78,7 @@ function Base.show(io::IO, @nospecialize(mat::BACMaterial))
 end
 
 function BACMaterial(; kernel::Function=linear_kernel,
-                     model::AbstractConstitutiveModel=NeoHookeNonlinear(),
+                     model::AbstractConstitutiveModel=LinearElastic(),
                      maxdmg::Real=0.85)
     return BACMaterial(kernel, model, maxdmg)
 end

src/physics/constitutive_models.jl

@@ -2,9 +2,6 @@ struct LinearElastic <: AbstractConstitutiveModel end
 
 function first_piola_kirchhoff(::LinearElastic, storage::AbstractStorage,
                                params::AbstractPointParameters, F::SMatrix{3,3,T,9}) where T
-    # J = det(F)
-    # J < eps() && return zero(SMatrix{3,3,T,9})
-    # isnan(J) && return zero(SMatrix{3,3,T,9})
     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.μ)
@@ -22,12 +19,11 @@ function get_hooke_matrix(nu, λ, μ)
     return CVoigt
 end
 
-struct NeoHookeViera <: AbstractConstitutiveModel end
+struct NeoHooke <: AbstractConstitutiveModel end
 
-function first_piola_kirchhoff(::NeoHookeViera, storage::AbstractStorage,
+function first_piola_kirchhoff(::NeoHooke, storage::AbstractStorage,
                                params::AbstractPointParameters, F::SMatrix{3,3,T,9}) where T
     J = det(F)
-    # J < eps() && return zero(SMatrix{3,3,T,9})
     Cinv = inv(F' * F)
     S = params.μ * (I - Cinv) + params.λ * log(J) * Cinv
     P = F * S
@@ -39,8 +35,8 @@ struct NeoHookeNonlinear <: AbstractConstitutiveModel end
 function first_piola_kirchhoff(::NeoHookeNonlinear, storage::AbstractStorage,
                                params::AbstractPointParameters, F::SMatrix{3,3,T,9}) where T
     J = det(F)
-    # J < eps() && return zero(SMatrix{3,3,T,9})
-    # isnan(J) && return zero(SMatrix{3,3,T,9})
+    J < eps() && return zero(SMatrix{3,3,T,9})
+    isnan(J) && return zero(SMatrix{3,3,T,9})
     C = F' * F
     Cinv = inv(C)
     S = params.G .* (I - 1 / 3 .* tr(C) .* Cinv) .* J^(-2 / 3) .+
@@ -53,9 +49,6 @@ struct SaintVenantKirchhoff <: AbstractConstitutiveModel end
 
 function first_piola_kirchhoff(::SaintVenantKirchhoff, storage::AbstractStorage,
                                params::AbstractPointParameters, F::SMatrix{3,3,T,9}) where T
-    # J = det(F)
-    # J < eps() && return zero(SMatrix{3,3,T,9})
-    # isnan(J) && return zero(SMatrix{3,3,T,9})
     E = 0.5 .* (F' * F - I)
     S = params.λ * tr(E) * I + 2 * params.μ * E
     P = F * S

src/physics/correspondence.jl

@@ -8,7 +8,7 @@ consistent (correspondence) formulation of non-ordinary state-based peridynamics
 - `kernel::Function`: Kernel function used for weighting the interactions between points.
     (default: `linear_kernel`)
 - `model::AbstractConstitutiveModel`: Constitutive model defining the material behavior.
-    (default: `NeoHookeNonlinear()`)
+    (default: `LinearElastic()`)
 - `zem::AbstractZEMStabilization`: Zero-energy mode stabilization. The
     stabilization algorithm of Silling (2017) is used as default.
     (default: `ZEMSilling()`)
@@ -32,12 +32,18 @@ CMaterial(maxdmg=0.85, zem=ZEMSilling())
 ---
 
 ```julia
-CMaterial
+CMaterial{CM,ZEM,K}
 ```
 
 Material type for the local continuum consistent (correspondence) formulation of
 non-ordinary state-based peridynamics.
 
+# Type Parameters
+- `CM`: A constitutive model type. See the constructor docs for more informations.
+- `ZEM`: A zero-energy mode stabilization type. See the constructor docs for more
+         informations.
+- `K`: A kernel function type. See the constructor docs for more informations.
+
 # Fields
 - `kernel::Function`: Kernel function used for weighting the interactions between points.
     See the constructor docs for more informations.
@@ -88,7 +94,7 @@ function Base.show(io::IO, @nospecialize(mat::CMaterial))
 end
 
 function CMaterial(; kernel::Function=linear_kernel,
-                    model::AbstractConstitutiveModel=NeoHookeNonlinear(),
+                    model::AbstractConstitutiveModel=LinearElastic(),
                     zem::AbstractZEMStabilization=ZEMSilling(), maxdmg::Real=0.85)
     return CMaterial(kernel, model, zem, maxdmg)
 end
@@ -158,7 +164,6 @@ end
 function force_density_point!(storage::AbstractStorage, system::AbstractSystem,
                               mat::AbstractCorrespondenceMaterial,
                               params::AbstractPointParameters, t, Δt, i)
-    calc_failure_point!(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)
@@ -167,18 +172,6 @@ function force_density_point!(storage::AbstractStorage, system::AbstractSystem,
     return nothing
 end
 
-function calc_failure_point!(storage, system, mat, params, i)
-    for bond_id in each_bond_idx(system, i)
-        bond = system.bonds[bond_id]
-        j, L = bond.neighbor, bond.length
-        Δxij = get_coordinates_diff(storage, i, j)
-        l = norm(Δxij)
-        ε = (l - L) / L
-        stretch_based_failure!(storage, system, bond, params, ε, i, bond_id)
-    end
-    return nothing
-end
-
 function calc_deformation_gradient(storage::CStorage, system::BondSystem, ::CMaterial,
                                    ::CPointParameters, i)
     (; bonds, volume) = system
@@ -226,9 +219,9 @@ function c_force_density!(storage::AbstractStorage, system::AbstractSystem,
         j, L = bond.neighbor, bond.length
         ΔXij = get_coordinates_diff(system, i, j)
         Δxij = get_coordinates_diff(storage, i, j)
-        # l = norm(Δxij)
-        # ε = (l - L) / L
-        # stretch_based_failure!(storage, system, bond, params, ε, i, bond_id)
+        l = norm(Δxij)
+        ε = (l - L) / L
+        stretch_based_failure!(storage, system, bond, params, ε, i, bond_id)
 
         # stabilization
         ωij = kernel(system, bond_id) * bond_active[bond_id]
@@ -245,13 +238,7 @@ end
 function too_much_damage!(storage::AbstractStorage, system::AbstractSystem,
                           mat::AbstractCorrespondenceMaterial, defgrad_res, i)
     (; F) = defgrad_res
-    # @infiltrate storage.n_active_bonds[i] ≤ 3
     # if storage.n_active_bonds[i] ≤ 3 || storage.damage[i] > mat.maxdmg || containsnan(F)
-    #     # kill all bonds of this point
-    #     storage.bond_active[each_bond_idx(system, i)] .= false
-    #     storage.n_active_bonds[i] = 0
-    #     return true
-    # end
     if storage.damage[i] > mat.maxdmg || containsnan(F)
         # kill all bonds of this point
         storage.bond_active[each_bond_idx(system, i)] .= false