Bond-associated Correspondence of Chen and Spencer

benchmark/benchmarks.jl

@@ -8,6 +8,7 @@ const SUITE = BenchmarkGroup()
 include(joinpath(@__DIR__, "mode_i.jl"))
 
 SUITE["mode_i"] = BenchmarkGroup()
-SUITE["mode_i"]["BBMaterial, 40"] = mode_i(BBMaterial(), 40)
-SUITE["mode_i"]["OSBMaterial, 40"] = mode_i(OSBMaterial(), 40)
-SUITE["mode_i"]["CCMaterial, 40"] = mode_i(NOSBMaterial(), 40)
+
+SUITE["mode_i"]["BBMaterial, 40"] = @benchmarkable submit(mode_i(BBMaterial(), 40))
+SUITE["mode_i"]["OSBMaterial, 40"] = @benchmarkable submit(mode_i(OSBMaterial(), 40))
+SUITE["mode_i"]["CMaterial, 40"] = @benchmarkable submit(mode_i(CMaterial(), 40))

benchmark/mode_i.jl

@@ -1,3 +1,4 @@
+using Peridynamics
 
 function mode_i(mat, npyz)
     l, Δx, δ, a = 1.0, 1/npyz, 3.015/npyz, 0.5
@@ -13,6 +14,7 @@ function mode_i(mat, npyz)
     velocity_bc!(t -> -30, body, :set_bottom, :y)
     velocity_bc!(t -> 30, body, :set_top, :y)
     vv = VelocityVerlet(steps=2000)
-    job = Job(body, vv; path=joinpath(@__DIR__, "results", "mode_i-BB-npyz$(npyz)"))
-    return @benchmarkable submit($job)
+    path = joinpath(@__DIR__, "results", "mode_i")
+    job = Job(body, vv; path)
+    return job
 end

docs/src/api_reference.md

@@ -12,14 +12,22 @@ Pages = ["api_reference.md"]
 ```@docs
 BBMaterial
 OSBMaterial
-NOSBMaterial
+CMaterial
+BACMaterial
 CKIMaterial
 ```
 
-## System related types
+## System or material related types
 ```@docs
 NoCorrection
 EnergySurfaceCorrection
+ZEMSilling
+LinearElastic
+NeoHooke
+MooneyRivlin
+SaintVenantKirchhoff
+linear_kernel
+cubic_b_spline_kernel
 ```
 
 ## Discretization

docs/src/literate/tutorial_wave_interface.jl

@@ -19,7 +19,7 @@ npyz = 4
 # With these parameters we now create a body, here using the non-ordinary state-based
 # correspondence formulation.
 pos, vol = uniform_box(lx, lyz, lyz, Δx)
-body = Body(NOSBMaterial(), pos, vol)
+body = Body(CMaterial(), pos, vol)
 # Again, failure is not allowed in the whole body.
 failure_permit!(body, false)
 

src/Peridynamics.jl

@@ -9,7 +9,13 @@ end
 import LibGit2, Dates
 
 # Material models
-export BBMaterial, OSBMaterial, NOSBMaterial, CKIMaterial
+export BBMaterial, OSBMaterial, CMaterial, BACMaterial, CKIMaterial
+
+# CMaterial related types
+export LinearElastic, NeoHooke, MooneyRivlin, SaintVenantKirchhoff, ZEMSilling
+
+# Kernels
+export linear_kernel, cubic_b_spline_kernel
 
 # Systems related types
 export NoCorrection, EnergySurfaceCorrection
@@ -50,6 +56,7 @@ abstract type AbstractTimeSolver end
 abstract type AbstractJob end
 abstract type AbstractJobOptions end
 abstract type AbstractSystem end
+abstract type AbstractBondSystem <: AbstractSystem end
 abstract type AbstractPredefinedCrack end
 abstract type AbstractBodyChunk{S<:AbstractSystem,T<:AbstractMaterial} end
 abstract type AbstractParameterHandler <: AbstractParameterSetup end
@@ -65,17 +72,25 @@ abstract type AbstractCorrection end
 abstract type AbstractStorage end
 abstract type AbstractCondition end
 abstract type AbstractBondSystemMaterial{Correction} <: AbstractMaterial end
+abstract type AbstractCorrespondenceMaterial{CM,ZEM} <: AbstractBondSystemMaterial{ZEM} end
+abstract type AbstractBondAssociatedSystemMaterial <: AbstractMaterial end
+abstract type AbstractConstitutiveModel end
+abstract type AbstractStressIntegration end
+abstract type AbstractZEMStabilization <: AbstractCorrection end
 abstract type AbstractInteractionSystemMaterial <: AbstractMaterial end
 
 include("auxiliary/function_arguments.jl")
 include("auxiliary/io.jl")
 include("auxiliary/logs.jl")
 include("auxiliary/mpi.jl")
 include("auxiliary/errors.jl")
+include("auxiliary/static_arrays.jl")
+include("auxiliary/nans.jl")
 
 include("physics/boundary_conditions.jl")
 include("physics/initial_conditions.jl")
 include("physics/material_parameters.jl")
+include("physics/stress.jl")
 include("physics/fracture.jl")
 include("physics/short_range_force_contact.jl")
 
@@ -86,8 +101,11 @@ include("discretization/body.jl")
 include("discretization/multibody_setup.jl")
 include("discretization/decomposition.jl")
 include("discretization/chunk_handler.jl")
+include("discretization/kernels.jl")
 include("discretization/bond_system.jl")
 include("discretization/bond_system_corrections.jl")
+include("discretization/zem_stabilization.jl")
+include("discretization/bond_associated_system.jl")
 include("discretization/interaction_system.jl")
 include("discretization/body_chunk.jl")
 
@@ -111,7 +129,9 @@ include("time_solvers/dynamic_relaxation.jl")
 include("physics/bond_based.jl")
 include("physics/continuum_kinematics_inspired.jl")
 include("physics/ordinary_state_based.jl")
+include("physics/constitutive_models.jl")
 include("physics/correspondence.jl")
+include("physics/ba_correspondence.jl")
 
 include("VtkReader/VtkReader.jl")
 using .VtkReader

src/auxiliary/nans.jl

@@ -0,0 +1,14 @@
+function containsnan(K::T) where {T<:AbstractArray}
+    @simd for i in eachindex(K)
+        isnan(K[i]) && return true
+    end
+    return false
+end
+
+function nancheck(chunk::AbstractBodyChunk, t)
+    if containsnan(chunk.storage.b_int)
+        msg = "NaN's found in field `b_int` at simulation time $(t)!\n"
+        error(msg)
+    end
+    return nothing
+end

src/auxiliary/precompile_workload.jl

@@ -5,12 +5,12 @@
         msg *= "Trigger package precompilation manually and then restart the mpirun!\n"
         error(msg)
     end
-    root = joinpath(@__DIR__, "temp_precompilation")
+    root = mktempdir()
     pos1, vol1 = uniform_box(1, 1, 1, 0.5; center=(0.5, 0.5, 0.5))
     pos2, vol2 = uniform_box(1, 1, 1, 0.5; center=(-0.5, 0.5, 0.5))
     path_bb = joinpath(root, "BB")
     path_osb = joinpath(root, "OSB")
-    path_nosb = joinpath(root, "NOSB")
+    path_cc = joinpath(root, "CC")
     path_ms = joinpath(root, "MS")
 
     @compile_workload begin
@@ -43,7 +43,7 @@
         velocity_bc!(t -> -1, b2, :set_a, :x)
         velocity_bc!(t -> 1, b2, :set_b, 1)
 
-        b3 = Body(NOSBMaterial(), pos1, vol1)
+        b3 = Body(CMaterial(), pos1, vol1)
         material!(b3; horizon=2, E=2.1e5, nu=0.25, rho=8e-6, Gc=2.7)
         point_set!(p -> p[1] ≤ 0.5, b3, :set_a)
         point_set!(x -> x > 0.5, b3, :set_b)
@@ -65,9 +65,7 @@
 
         submit(Job(b1, vv; path=path_bb, freq=1); quiet=true)
         submit(Job(b2, vv; path=path_osb, freq=1); quiet=true)
-        submit(Job(b3, vv; path=path_nosb, freq=1); quiet=true)
+        submit(Job(b3, vv; path=path_cc, freq=1); quiet=true)
         submit(Job(ms, vv; path=path_ms, freq=1); quiet=true)
     end
-
-    rm(root; recursive=true, force=true)
 end

src/auxiliary/static_arrays.jl

@@ -0,0 +1,44 @@
+
+@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])
+    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]
+    return nothing
+end
+
+@inline function get_vector(M::AbstractMatrix, i::Int)
+    V = SVector{3}(M[1, i], M[2, i], M[3, i])
+    return V
+end
+
+@inline function update_vector!(Mₙ::AbstractMatrix, i::Int, Vₙ₊₁::SVector{3})
+    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)
+    V = SVector{3}(M[1, j] - M[1, i], M[2, j] - M[2, i], M[3, j] - M[3, i])
+    return V
+end
+
+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)
+    Sinv = Diagonal{T,SVector{N,T}}(Sinvreg)
+    Minv = V * Sinv * U'
+    return Minv
+end

src/core/mpi_body_data_handler.jl

@@ -341,6 +341,14 @@ function recv_htl!(get_field_function::F, dh::MPIBodyDataHandler,
     return nothing
 end
 
+function calc_force_density!(dh::MPIBodyDataHandler, t, Δt)
+    (; chunk) = dh
+    @timeit_debug TO "exchange_loc_to_halo!" exchange_loc_to_halo!(dh)
+    @timeit_debug TO "calc_force_density!" calc_force_density!(chunk, t, Δt)
+    @timeit_debug TO "exchange_halo_to_loc!" exchange_halo_to_loc!(dh)
+    return nothing
+end
+
 function export_results(dh::MPIBodyDataHandler, options::AbstractJobOptions, n::Int,
                         t::Float64)
     options.export_allowed || return nothing

src/core/parameters.jl

@@ -48,6 +48,7 @@ function macrocheck_input_material(material)
     material isa Symbol && return nothing
     (material isa Expr && material.head === :.) && return nothing
     (material isa Expr && material.head === :escape) && return nothing
+    (material isa Expr && material.head === :curly) && return nothing
     return throw(ArgumentError("argument `$material` is not a valid material input!\n"))
 end
 

src/core/threads_body_data_handler.jl

@@ -157,6 +157,17 @@
     return nothing
 end
 
+function calc_force_density!(dh::ThreadsBodyDataHandler, t, Δt)
+    @threads :static for chunk_id in eachindex(dh.chunks)
+        exchange_loc_to_halo!(dh, chunk_id)
+        calc_force_density!(dh.chunks[chunk_id], t, Δt)
+    end
+    @threads :static for chunk_id in eachindex(dh.chunks)
+        exchange_halo_to_loc!(dh, chunk_id)
+    end
+    return nothing
+end
+
 function export_results(dh::ThreadsBodyDataHandler, options::AbstractJobOptions,
                         chunk_id::Int, timestep::Int, time::Float64)
     options.export_allowed || return nothing

src/core/threads_multibody_data_handler.jl

@@ -33,6 +33,17 @@
 @inline each_body_name(dh::ThreadsMultibodyDataHandler) = dh.body_names
 @inline each_body_idx(dh::ThreadsMultibodyDataHandler) = eachindex(dh.body_dhs)
 
+function calc_force_density!(dh::ThreadsMultibodyDataHandler, t, Δt)
+    for body_idx in each_body_idx(dh)
+        body_dh = get_body_dh(dh, body_idx)
+        @threads :static for chunk_id in eachindex(body_dh.chunks)
+            exchange_loc_to_halo!(body_dh, chunk_id)
+            calc_force_density!(body_dh.chunks[chunk_id], t, Δt)
+        end
+    end
+    return nothing
+end
+
 function update_caches!(dh::ThreadsMultibodyDataHandler)
     for body_idx in each_body_idx(dh)
         body_dh = get_body_dh(dh, body_idx)

src/discretization/body.jl

@@ -225,6 +225,10 @@ function pre_submission_check(body::Body; body_in_multibody_setup::Bool=false)
     return nothing
 end
 
+@inline function get_point_param(b::AbstractBody, i::Int)
+    return b.point_params[b.params_map[i]]
+end
+
 @inline function get_point_param(b::AbstractBody, key::Symbol, i::Int)
     return getfield(b.point_params[b.params_map[i]], key)
 end