Merge pull request #220 from kaipartmann/rk

Reproducing kernel peridynamics formulation

Author: kaipartmann

docs/src/private_api_reference.md

@@ -25,6 +25,7 @@ Peridynamics.apply_precrack!
 Peridynamics.point_sets_intersect
 Peridynamics.velocity_databc!
 Peridynamics.forcedensity_databc!
+Peridynamics.invreg
 
 Peridynamics.InterfaceError
 Peridynamics.HaloExchange
@@ -55,4 +56,5 @@ Peridynamics.PosDepSingleDimIC
 Peridynamics.ShortRangeForceContact
 
 Peridynamics.@storage
+Peridynamics.@autoinfiltrate
 ```

docs/src/public_api_reference.md

@@ -17,6 +17,8 @@ CMaterial
 CRMaterial
 BACMaterial
 CKIMaterial
+RKCMaterial
+RKCRMaterial
 ```
 
 ## System or material related types
@@ -30,8 +32,10 @@ LinearElastic
 NeoHooke
 MooneyRivlin
 SaintVenantKirchhoff
+const_one_kernel
 linear_kernel
 cubic_b_spline_kernel
+cubic_b_spline_kernel_norm
 ```
 
 ## Discretization

src/Peridynamics.jl

@@ -10,13 +10,13 @@ import LibGit2, Dates
 
 # Material models
 export BBMaterial, DHBBMaterial, OSBMaterial, CMaterial, CRMaterial, BACMaterial,
-       CKIMaterial
+       CKIMaterial, RKCMaterial, RKCRMaterial
 
 # CMaterial related types
 export LinearElastic, NeoHooke, MooneyRivlin, SaintVenantKirchhoff, ZEMSilling, ZEMWan
 
 # Kernels
-export linear_kernel, cubic_b_spline_kernel
+export const_one_kernel, linear_kernel, cubic_b_spline_kernel, cubic_b_spline_kernel_norm
 
 # Damage models
 export CriticalStretch
@@ -78,6 +78,7 @@ 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 AbstractRKCMaterial{CM,C} <: AbstractBondSystemMaterial{C} end
 abstract type AbstractBondAssociatedSystemMaterial <: AbstractMaterial end
 abstract type AbstractConstitutiveModel end
 abstract type AbstractStressIntegration end
@@ -93,6 +94,7 @@ include("auxiliary/mpi.jl")
 include("auxiliary/errors.jl")
 include("auxiliary/static_arrays.jl")
 include("auxiliary/nans.jl")
+include("auxiliary/autoinfiltrate.jl")
 
 include("physics/boundary_conditions.jl")
 include("physics/initial_conditions.jl")
@@ -140,6 +142,8 @@ include("physics/ordinary_state_based.jl")
 include("physics/constitutive_models.jl")
 include("physics/correspondence.jl")
 include("physics/correspondence_rotated.jl")
+include("physics/rk_correspondence.jl")
+include("physics/rk_correspondence_rotated.jl")
 include("physics/ba_correspondence.jl")
 
 include("VtkReader/VtkReader.jl")

src/auxiliary/autoinfiltrate.jl

@@ -0,0 +1,35 @@
+"""
+    @autoinfiltrate
+    @autoinfiltrate condition::Bool
+
+$(internal_api_warning())
+
+Invoke the `@infiltrate` macro of the package Infiltrator.jl to create a breakpoint for
+ad-hoc interactive debugging in the REPL. If the optional argument `condition` is given, the
+breakpoint is only enabled if `condition` evaluates to `true`.
+
+As opposed to using `Infiltrator.@infiltrate` directly, this macro does not require
+Infiltrator.jl to be added as a dependency to Peridynamics.jl. As a bonus, the macro will
+also attempt to load the Infiltrator module if it has not yet been loaded manually.
+
+Note: For this macro to work, the Infiltrator.jl package needs to be installed in your
+current Julia environment stack.
+
+See also: [Infiltrator.jl](https://github.com/JuliaDebug/Infiltrator.jl)
+"""
+macro autoinfiltrate(condition=true)
+    pkgid = Base.PkgId(Base.UUID("5903a43b-9cc3-4c30-8d17-598619ec4e9b"), "Infiltrator")
+    if !haskey(Base.loaded_modules, pkgid)
+        try
+            Base.eval(Main, :(using Infiltrator))
+        catch
+        end
+    end
+    i = get(Base.loaded_modules, pkgid, nothing)
+    lnn = LineNumberNode(__source__.line, __source__.file)
+    error_msg = "Cannot load Infiltrator.jl!"
+    error_msg *= "Make sure it is included in your environment stack."
+    isnothing(i) && return Expr(:macrocall, Symbol("@warn"), lnn, error_msg)
+    return Expr(:macrocall, Expr(:., i, QuoteNode(Symbol("@infiltrate"))), lnn,
+                esc(condition))
+end

src/auxiliary/nans.jl

@@ -5,9 +5,10 @@ function containsnan(K::T) where {T<:AbstractArray}
     return false
 end
 
-function nancheck(chunk::AbstractBodyChunk, t)
+function nancheck(chunk::AbstractBodyChunk, t, Δt)
     if containsnan(chunk.storage.b_int)
-        msg = "NaN's found in field `b_int` at simulation time $(t)!\n"
+        n = (t > 0 && Δt > 0) ? Int(t ÷ Δt) : 0
+        msg = "NaN's found in field `b_int` at simulation time $(t), step $(n)!\n"
         error(msg)
     end
     return nothing

src/auxiliary/static_arrays.jl

@@ -1,45 +1,64 @@
 
-@inline function get_tensor(T::AbstractMatrix, i::Int)
-    tensor = SMatrix{3,3}(T[1, i], T[2, i], T[3, i], T[4, i], T[5, i], T[6, i], T[7, i],
-                          T[8, i], T[9, i])
+@inline function get_tensor(Mₙ::AbstractMatrix{T}, i::Int) where {T}
+    tensor = SMatrix{3,3,T,9}(Mₙ[1, i], Mₙ[2, i], Mₙ[3, i], Mₙ[4, i], Mₙ[5, i], Mₙ[6, i],
+                              Mₙ[7, i], Mₙ[8, i], Mₙ[9, i])
     return tensor
 end
 
-@inline function update_tensor!(Tₙ::AbstractMatrix, i::Int, Tₙ₊₁::SMatrix{3,3})
-    Tₙ[1, i] = Tₙ₊₁[1, 1]
-    Tₙ[2, i] = Tₙ₊₁[1, 2]
-    Tₙ[3, i] = Tₙ₊₁[1, 3]
-    Tₙ[4, i] = Tₙ₊₁[2, 1]
-    Tₙ[5, i] = Tₙ₊₁[2, 2]
-    Tₙ[6, i] = Tₙ₊₁[2, 3]
-    Tₙ[7, i] = Tₙ₊₁[3, 1]
-    Tₙ[8, i] = Tₙ₊₁[3, 2]
-    Tₙ[9, i] = Tₙ₊₁[3, 3]
+@inline function update_tensor!(Mₙ::AbstractMatrix{T}, i::Int,
+                                Tₙ₊₁::StaticMatrix{3,3,T}) where {T}
+    Mₙ[1, i] = Tₙ₊₁[1]
+    Mₙ[2, i] = Tₙ₊₁[2]
+    Mₙ[3, i] = Tₙ₊₁[3]
+    Mₙ[4, i] = Tₙ₊₁[4]
+    Mₙ[5, i] = Tₙ₊₁[5]
+    Mₙ[6, i] = Tₙ₊₁[6]
+    Mₙ[7, i] = Tₙ₊₁[7]
+    Mₙ[8, i] = Tₙ₊₁[8]
+    Mₙ[9, i] = Tₙ₊₁[9]
     return nothing
 end
 
-@inline function get_vector(M::AbstractMatrix, i::Int)
-    return SVector{3}(M[1, i], M[2, i], M[3, i])
+@inline function get_vector(M::AbstractMatrix{T}, i::Int) where {T}
+    return SVector{3,T}(M[1, i], M[2, i], M[3, i])
 end
 
-@inline function update_vector!(Mₙ::AbstractMatrix, i::Int, Vₙ₊₁::SVector{3})
+@inline function update_vector!(Mₙ::AbstractMatrix{T}, i::Int,
+                                Vₙ₊₁::StaticVector{3,T}) where {T}
     Mₙ[1, i] = Vₙ₊₁[1]
     Mₙ[2, i] = Vₙ₊₁[2]
     Mₙ[3, i] = Vₙ₊₁[3]
     return nothing
 end
 
-@inline function update_add_vector!(Mₙ::AbstractMatrix, i::Int, Vₙ₊₁::SVector{3})
+@inline function update_add_vector!(Mₙ::AbstractMatrix{T}, i::Int,
+                                    Vₙ₊₁::StaticVector{3,T}) where {T}
     Mₙ[1, i] += Vₙ₊₁[1]
     Mₙ[2, i] += Vₙ₊₁[2]
     Mₙ[3, i] += Vₙ₊₁[3]
     return nothing
 end
 
-@inline function get_vector_diff(M::AbstractMatrix, i::Int, j::Int)
-    return SVector{3}(M[1, j] - M[1, i], M[2, j] - M[2, i], M[3, j] - M[3, i])
+@inline function get_vector_diff(M::AbstractMatrix{T}, i::Int, j::Int) where {T}
+    return SVector{3,T}(M[1, j] - M[1, i], M[2, j] - M[2, i], M[3, j] - M[3, i])
 end
 
+"""
+    invreg(M::StaticMatrix{N,N,T}, threshold::Real=eps()) where {N,T}
+
+Computes the pseudo-inverse of a square static matrix `M` using singular value
+decomposition (SVD) with regularization. The regularization is applied to the
+singular values, where singular values below the specified `threshold` are set
+to zero in the pseudo-inverse. This helps to avoid numerical instability and
+ill-conditioning in the inversion process.
+
+# Arguments
+- `M::StaticMatrix{N,N,T}`: The square `N`×`N` static matrix with element type `T` to be
+    inverted.
+- `threshold::Real=eps()`: The threshold value for regularization, should be a positive
+    real number between 0 and 1. Singular values below this threshold are set to zero in the
+    pseudo-inverse.
+"""
 function invreg(M::StaticMatrix{N,N,T}, threshold::Real=eps()) where {N,T}
     U, S, V = svd(M)
     Sinvreg = SVector{N,T}((s > threshold ? 1/s : 0) for s in S)

src/discretization/bond_system.jl

@@ -116,7 +116,7 @@ function find_bonds!(bonds::Vector{Bond}, nhs::PointNeighbors.GridNeighborhoodSe
 end
 
 @inline function check_point_duplicates(L::Float64, i::Int, j::Int)
-    if L == 0
+    if L < eps()
         msg = "point duplicate found!\n"
         msg *= "Point #$(i) has a duplicate #$(j) which will lead to `NaN`s!\n"
         error(msg)
@@ -232,18 +232,18 @@ end
 function break_bonds!(storage::AbstractStorage, system::AbstractBondSystem,
                       set_a::Vector{Int}, set_b::Vector{Int})
     storage.n_active_bonds .= 0
-    for point_id in each_point_idx(system)
-        for bond_id in each_bond_idx(system, point_id)
+    for i in each_point_idx(system)
+        for bond_id in each_bond_idx(system, i)
             bond = system.bonds[bond_id]
             neighbor_id = bond.neighbor
-            point_in_a = in(point_id, set_a)
-            point_in_b = in(point_id, set_b)
+            point_in_a = in(i, set_a)
+            point_in_b = in(i, set_b)
             neigh_in_a = in(neighbor_id, set_a)
             neigh_in_b = in(neighbor_id, set_b)
             if (point_in_a && neigh_in_b) || (point_in_b && neigh_in_a)
                 storage.bond_active[bond_id] = false
             end
-            storage.n_active_bonds[point_id] += storage.bond_active[bond_id]
+            storage.n_active_bonds[i] += storage.bond_active[bond_id]
         end
     end
     return nothing
@@ -264,12 +264,22 @@ function calc_force_density!(chunk::AbstractBodyChunk{<:AbstractBondSystem}, t,
     (; dmgmodel) = mat
     storage.b_int .= 0
     storage.n_active_bonds .= 0
-    for point_id in each_point_idx(chunk)
-        calc_failure!(storage, system, mat, dmgmodel, paramsetup, point_id)
-        calc_damage!(storage, system, mat, dmgmodel, paramsetup, point_id)
-        force_density_point!(storage, system, mat, paramsetup, t, Δt, point_id)
+    for i in each_point_idx(chunk)
+        calc_failure!(storage, system, mat, dmgmodel, paramsetup, i)
+        calc_damage!(storage, system, mat, dmgmodel, paramsetup, i)
+        force_density_point!(storage, system, mat, paramsetup, t, Δt, i)
     end
-    nancheck(chunk, t)
+    nancheck(chunk, t, Δt)
+    return nothing
+end
+
+# a placeholder function for all force density calculations with multiple parameters that
+# are not specifially handled by the material
+function force_density_point!(storage::AbstractStorage, system::AbstractBondSystem,
+                              mat::AbstractBondSystemMaterial,
+                              paramhandler::AbstractParameterHandler, t, Δt, i)
+    params = get_params(paramhandler, i)
+    force_density_point!(storage, system, mat, params, t, Δt, i)
     return nothing
 end
 
@@ -393,7 +403,9 @@ function log_material_property(::Val{:kernel}, mat::AbstractBondSystemMaterial;
     return msg
 end
 
-function log_material(mat::M; indentation::Int=2) where {M<:AbstractCorrespondenceMaterial}
+function log_material(mat::M;
+                      indentation::Int=2) where {M<:Union{AbstractCorrespondenceMaterial,
+                                                          AbstractRKCMaterial}}
     msg = msg_qty("material type", nameof(M); indentation)
     for prop in fieldnames(M)
         msg *= log_material_property(Val(prop), mat; indentation)

src/discretization/interaction_system.jl

@@ -427,7 +427,7 @@ function calc_force_density!(chunk::AbstractBodyChunk{<:InteractionSystem}, t, 
         calc_damage!(storage, system, mat, dmgmodel, paramsetup, point_id)
         force_density_point!(storage, system, mat, paramsetup, t, Δt, point_id)
     end
-    nancheck(chunk, t)
+    nancheck(chunk, t, Δt)
     return nothing
 end
 

src/discretization/kernels.jl

@@ -1,3 +1,14 @@
+@doc raw"""
+    const_one_kernel(δ, L)
+
+A kernel function ``\omega`` (also called influence function) used for weighting the
+bonds in a family. The kernel function is simply defined as a constant value 1:
+```math
+\omega = 1
+```
+"""
+@inline const_one_kernel(δ, L) = 1
+
 @doc raw"""
     linear_kernel(δ, L)
 
@@ -39,3 +50,35 @@ with the horizon ``\delta`` and the initial bond length ``L``.
         return 0
     end
 end
+
+@doc raw"""
+    cubic_b_spline_kernel_norm(δ, L)
+
+A cubic B-spline kernel function ``\omega`` used for weighting the bonds in a family.
+The kernel function is defined as
+```math
+\begin{aligned}
+\xi &= \frac{L}{\delta} \; , \\
+\omega &= \frac{8}{\pi \, \delta^3} \cdot \left\{
+    \begin{array}{ll}
+        1 - 6 \xi^2 + 6 \xi^3 & \quad \text{if} \; 0 < \xi \leq 0.5 \; , \\[3pt]
+        2 (1 - \xi)^3 & \quad \text{if} \; 0.5 < \xi \leq 1 \; , \\[3pt]
+        0 & \quad \text{else} \; ,
+    \end{array}
+    \right.
+\end{aligned}
+```
+with the horizon ``\delta`` and the initial bond length ``L``. This kernel is properly
+normalized to satisfy the condition ``\int_{\mathcal{H}(X)} \omega(\Delta X) dV' = 1``.
+"""
+@inline function cubic_b_spline_kernel_norm(δ, L)
+    ξ = L / δ
+    C = 8.0/(π * δ^3)
+    if 0 < ξ ≤ 0.5
+        return C * (1 - 6 * ξ^2 + 6 * ξ^3)
+    elseif 0.5 < ξ ≤ 1
+        return C * 2 * (1 - ξ)^3
+    else
+        return 0
+    end
+end

src/physics/ba_correspondence.jl

@@ -189,14 +189,6 @@ function init_field(::BACMaterial, ::AbstractTimeSolver, system::BondAssociatedS
     return zeros(get_n_loc_points(system))
 end
 
-function force_density_point!(storage::BACStorage, system::BondAssociatedSystem,
-                              mat::BACMaterial, paramhandler::AbstractParameterHandler,
-                              t, Δt, i)
-    params = get_params(paramhandler, i)
-    force_density_point!(storage, system, mat, params, t, Δt, i)
-    return nothing
-end
-
 function force_density_point!(storage::BACStorage, system::BondAssociatedSystem,
                               mat::BACMaterial, params::BACPointParameters, t, Δt, i)
     for bond_idx in each_bond_idx(system, i)