Refactoring, docstrings & test fixes

Author: kaipartmann

src/Peridynamics.jl

@@ -79,7 +79,7 @@ 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 AbstractBondAssociatedSystemMaterial <: AbstractBondSystemMaterial{Nothing} end
 abstract type AbstractConstitutiveModel end
 abstract type AbstractStressIntegration end
 abstract type AbstractZEMStabilization <: AbstractCorrection end

src/discretization/interaction_system.jl

@@ -424,7 +424,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, Δt)
+    nancheck(storage, t, Δt)
     return nothing
 end
 

src/physics/bond_based.jl

@@ -69,6 +69,7 @@ When specifying the `fields` keyword of [`Job`](@ref) for a [`Body`](@ref) with
 - `b_ext::Matrix{Float64}`: External force density of each point.
 - `damage::Vector{Float64}`: Damage of each point.
 - `n_active_bonds::Vector{Int}`: Number of intact bonds of each point.
+- `strain_energy_density::Vector{Float64}`: Strain energy density of each point.
 """
 struct BBMaterial{Correction,DM} <: AbstractBondBasedMaterial{Correction}
     dmgmodel::DM
@@ -200,8 +201,7 @@ function force_density_point!(storage::BBStorage, system::BondSystem, ::BBMateri
 end
 
 function strain_energy_density_point!(storage::AbstractStorage, system::BondSystem,
-                                      ::AbstractBondBasedMaterial,
-                                      paramsetup::AbstractParameterSetup, i)
+                                      ::BBMaterial, paramsetup::AbstractParameterSetup, i)
     (; bond_active, bond_length, strain_energy_density) = storage
     (; bonds, correction, volume) = system
     params_i = get_params(paramsetup, i)

src/physics/correspondence.jl

@@ -103,7 +103,7 @@ When specifying the `fields` keyword of [`Job`](@ref) for a [`Body`](@ref) with
 - `b_ext::Matrix{Float64}`: External force density of each point.
 - `damage::Vector{Float64}`: Damage of each point.
 - `n_active_bonds::Vector{Int}`: Number of intact bonds of each point.
-- `stress::Matrix{Float64}`: Stress tensor of each point.
+- `cauchy_stress::Matrix{Float64}`: Cauchy stress tensor of each point.
 - `von_mises_stress::Vector{Float64}`: Von Mises stress of each point.
 """
 struct CMaterial{CM,ZEM,K,DM} <: AbstractCorrespondenceMaterial{CM,ZEM}

src/physics/dh_bond_based.jl

@@ -69,6 +69,7 @@ When specifying the `fields` keyword of [`Job`](@ref) for a [`Body`](@ref) with
 - `b_ext::Matrix{Float64}`: External force density of each point.
 - `damage::Vector{Float64}`: Damage of each point.
 - `n_active_bonds::Vector{Int}`: Number of intact bonds of each point.
+- `strain_energy_density::Vector{Float64}`: Strain energy density of each point.
 """
 struct DHBBMaterial{Correction,DM} <: AbstractBondBasedMaterial{Correction}
     dmgmodel::DM
@@ -146,3 +147,23 @@ function force_density_point!(storage::DHBBStorage, system::BondSystem, ::DHBBMa
     end
     return nothing
 end
+
+function strain_energy_density_point!(storage::AbstractStorage, system::BondSystem,
+                                      ::DHBBMaterial, paramsetup::AbstractParameterSetup, i)
+    (; bond_active, bond_length, strain_energy_density) = storage
+    (; bonds, correction, volume) = system
+    params_i = get_params(paramsetup, i)
+    Ψ = 0.0
+    for bond_id in each_bond_idx(system, i)
+        bond = bonds[bond_id]
+        j, L = bond.neighbor, bond.length
+        l = bond_length[bond_id]
+        ε = (l - L) / L
+        params_j = get_params(paramsetup, j)
+        ωij = bond_active[bond_id] * surface_correction_factor(correction, bond_id)
+        bc = (params_i.bc + params_j.bc) / 2
+        Ψ += 0.5 * ωij * bc * ε * ε * L * volume[j] # added factor 2 here!
+    end
+    strain_energy_density[i] = Ψ
+    return nothing
+end

src/physics/ordinary_state_based.jl

@@ -77,6 +77,7 @@ When specifying the `fields` keyword of [`Job`](@ref) for a [`Body`](@ref) with
 - `b_ext::Matrix{Float64}`: External force density of each point.
 - `damage::Vector{Float64}`: Damage of each point.
 - `n_active_bonds::Vector{Int}`: Number of intact bonds of each point.
+- `strain_energy_density::Vector{Float64}`: Strain energy density of each point.
 """
 struct OSBMaterial{Correction,K,DM} <: AbstractBondSystemMaterial{Correction}
     kernel::K
@@ -226,33 +227,39 @@ function calc_dilatation(storage::OSBStorage, system::BondSystem, mat::OSBMateri
     return dil
 end
 
-function export_field(::Val{:strain_energy_density}, mat::OSBMaterial, system::BondSystem,
-                      storage::AbstractStorage, paramsetup::AbstractParameterSetup, t)
+function strain_energy_density_point!(storage::AbstractStorage, system::BondSystem,
+                                      mat::OSBMaterial, paramsetup::AbstractParameterSetup,
+                                      i)
     (; bond_active, bond_length, strain_energy_density) = storage
     (; bonds, correction, volume) = system
-    strain_energy_density .= 0.0
+    params_i = get_params(paramsetup, i)
+    wvol = calc_weighted_volume(storage, system, mat, params_i, i)
+    iszero(wvol) && return nothing
+    dil = calc_dilatation(storage, system, mat, params_i, wvol, i)
+    Ψvol = 0.5 * params_i.K * dil^2
+    Ψdev = 0.0
+    for bond_id in each_bond_idx(system, i)
+        bond = bonds[bond_id]
+        j, L = bond.neighbor, bond.length
+        l = bond_length[bond_id]
+        e = l - L
+        edev = e - 1/3 * dil * L
+        ωij = kernel(system, bond_id) * bond_active[bond_id]
+        β = surface_correction_factor(correction, bond_id)
+        params_j = get_params(paramsetup, j)
+        G = (params_i.G + params_j.G) / 2
+        cdev = 15.0 * G / (2 * wvol)
+        Ψdev += cdev * ωij * β * edev * edev * volume[j]
+    end
+    Ψ = Ψvol + Ψdev
+    strain_energy_density[i] = Ψ
+    return nothing
+end
+
+function export_field(::Val{:strain_energy_density}, mat::OSBMaterial, system::BondSystem,
+                      storage::AbstractStorage, paramsetup::AbstractParameterSetup, t)
     for i in each_point_idx(system)
-        params_i = get_params(paramsetup, i)
-        wvol = calc_weighted_volume(storage, system, mat, params_i, i)
-        iszero(wvol) && continue
-        dil = calc_dilatation(storage, system, mat, params_i, wvol, i)
-        Ψvol = 0.5 * params_i.K * dil^2
-        Ψdev = 0.0
-        for bond_id in each_bond_idx(system, i)
-            bond = bonds[bond_id]
-            j, L = bond.neighbor, bond.length
-            l = bond_length[bond_id]
-            e = l - L
-            edev = e - 1/3 * dil * L
-            ωij = kernel(system, bond_id) * bond_active[bond_id]
-            β = surface_correction_factor(correction, bond_id)
-            params_j = get_params(paramsetup, j)
-            G = (params_i.G + params_j.G) / 2
-            cdev = 15.0 * G / (2 * wvol)
-            Ψdev += cdev * ωij * β * edev * edev * volume[j]
-        end
-        Ψ = Ψvol + Ψdev
-        strain_energy_density[i] = Ψ
+        strain_energy_density_point!(storage, system, mat, paramsetup, i)
     end
-    return strain_energy_density
+    return storage.strain_energy_density
 end

src/physics/rk_correspondence.jl

@@ -447,13 +447,14 @@ function rkc_defgrad!(storage::AbstractStorage, system::AbstractBondSystem,
     return nothing
 end
 
-function calc_force_density!(chunk::BodyChunk{<:BondSystem,<:AbstractRKCMaterial}, t, Δt)
-    (; system, mat, paramsetup, storage) = chunk
-    storage.b_int .= 0
-    for i in each_point_idx(chunk)
+function calc_force_density!(storage::AbstractStorage, system::AbstractBondSystem,
+                             mat::AbstractRKCMaterial, paramsetup::AbstractParameterSetup,
+                             t, Δt)
+    storage.b_int .= 0.0
+    for i in each_point_idx(system)
         force_density_point!(storage, system, mat, paramsetup, t, Δt, i)
     end
-    nancheck(chunk, t, Δt)
+    nancheck(storage, t, Δt)
     return nothing
 end
 

test/auxiliary/test_nans.jl

@@ -17,11 +17,11 @@ end
     body = Body(BBMaterial(), ref_position, volume)
     material!(body, horizon=δ, rho=1, E=E)
     dh = Peridynamics.threads_data_handler(body, VelocityVerlet(steps=1), 1)
-    chunk = dh.chunks[1]
-    (; b_int) = chunk.storage
+    (; storage) = dh.chunks[1]
+    (; b_int) = storage
 
-    @test Peridynamics.nancheck(chunk, 0.0, 0.0) === nothing
+    @test Peridynamics.nancheck(storage, 0.0, 0.0) === nothing
 
     b_int[3, end] = NaN
-    @test_throws ErrorException Peridynamics.nancheck(chunk, 0.0, 0.0)
+    @test_throws ErrorException Peridynamics.nancheck(storage, 0.0, 0.0)
 end

test/physics/test_strain_energy_density.jl

@@ -36,7 +36,6 @@
         @printf("  mean error:      %7.2f %%\n", (Ψ̂_pd - Ψ_a) / Ψ_a * 100)
         @printf("  min rel. error:  %7.2f %%\n", ΔΨ_min * 100)
         @printf("  max rel. error:  %7.2f %%\n", ΔΨ_max * 100)
-        println("-"^50)
 
         @test ΔΨ_min > tols[1]
         @test ΔΨ_max < tols[2]
@@ -80,6 +79,42 @@ end
     test_stendens(body, ts, F_a, Ψ_a, tols; testcase)
 end
 
+@testitem "Strain energy density export DHBBMaterial" setup=[PsiExport] begin
+    Δx = 0.2
+    horizon = 3.01Δx
+    E = 210e9
+    nu = 0.25
+    λ = 1.05
+    pos, vol = uniform_box(1,1,1,Δx)
+    mat = DHBBMaterial{NoCorrection}()
+    body = Body(mat, pos, vol)
+    material!(body; horizon, rho=8000, E, nu)
+    params = body.point_params[1]
+    ts = VelocityVerlet(steps=1)
+
+    testcase = "homogeneous isotropic extension"
+    ε = λ - 1
+    F_a = @SMatrix [λ 0 0; 0 λ 0; 0 0 λ]
+    Ψ_a = 9/2 * params.K * ε^2
+    tols = (-0.9, 0.3)
+    test_stendens(body, ts, F_a, Ψ_a, tols; testcase)
+
+    testcase = "pure shear deformation"
+    β = 0.1
+    F_a = @SMatrix [1 β 0; 0 1 0; 0 0 1]
+    Ψ_a = 1/2 * params.G * β^2
+    tols = (-0.9, 0.3)
+    test_stendens(body, ts, F_a, Ψ_a, tols; testcase)
+
+    testcase = "uniform extension in x-direction"
+    λ = 1.1
+    ε = λ - 1
+    F_a = @SMatrix [λ 0 0; 0 1 0; 0 0 1]
+    Ψ_a = 3/5 * params.E * ε^2
+    tols = (-0.9, 0.3)
+    test_stendens(body, ts, F_a, Ψ_a, tols; testcase)
+end
+
 @testitem "Strain energy density export OSBMaterial" setup=[PsiExport] begin
     Δx = 0.2
     horizon = 3.01Δx