Add custom fields for some materials

Author: kaipartmann

src/auxiliary/io.jl

@@ -160,8 +160,9 @@ end
 end
 
 function export_fields!(vtk, chunk, fields::Vector{Symbol}, t)
+    (; mat, system, storage, paramsetup) = chunk
     for field in fields
-        point_data = export_field(Val(field), chunk.mat, chunk.system, chunk.storage, t)
+        point_data = export_field(Val(field), mat, system, storage, paramsetup, t)
         vtk[string(field), VTKPointData()] = point_data
     end
     return nothing
@@ -174,7 +175,7 @@ function export_fields!(vtk, chunk, fields_spec::Dict{Symbol,Vector{Symbol}}, t)
 end
 
 # this function can be specialized for each field, even custom export fields can be written!
-function export_field(::Val{field}, mat, system, storage, t) where {field}
+function export_field(::Val{field}, mat, system, storage, paramsetup, t) where {field}
     return get_loc_point_data(storage, system, field)
 end
 

src/physics/bond_based.jl

@@ -125,30 +125,36 @@ end
     @pointfield b_ext::Matrix{Float64}
     @pointfield density_matrix::Matrix{Float64}
     @pointfield damage::Vector{Float64}
-    bond_stretch::Vector{Float64}
+    bond_length::Vector{Float64}
     bond_active::Vector{Bool}
     @pointfield n_active_bonds::Vector{Int}
+    @pointfield strain_energy_density::Vector{Float64}
 end
 
-function init_field(::BBMaterial, ::AbstractTimeSolver, system::BondSystem,
-                    ::Val{:bond_stretch})
+function init_field(::AbstractBondBasedMaterial, ::AbstractTimeSolver, system::BondSystem,
+                    ::Val{:bond_length})
     return zeros(get_n_bonds(system))
 end
 
+function init_field(::AbstractBondBasedMaterial, ::AbstractTimeSolver, system::BondSystem,
+                    ::Val{:strain_energy_density})
+    return zeros(get_n_loc_points(system))
+end
+
 # Customized calc_failure to save the bond stretch ε for force density calculation
 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
+    (; position, n_active_bonds, bond_active, bond_length) = storage
     (; bonds) = system
     for bond_id in each_bond_idx(system, i)
         bond = bonds[bond_id]
         j, L = bond.neighbor, bond.length
         Δxij = get_vector_diff(position, i, j)
         l = norm(Δxij)
         ε = (l - L) / L
-        bond_stretch[bond_id] = ε / l # note that this is  ε / l!
+        bond_length[bond_id] = l # store current bond length
         if ε > εc && bond.fail_permit
             bond_active[bond_id] = false
         end
@@ -159,34 +165,66 @@ end
 
 function force_density_point!(storage::BBStorage, system::BondSystem, ::BBMaterial,
                               params::StandardPointParameters, t, Δt, i)
-    (; position, bond_stretch, bond_active, b_int) = storage
+    (; position, bond_length, bond_active, b_int) = storage
     (; bonds, correction, volume) = system
     for bond_id in each_bond_idx(system, i)
         bond = bonds[bond_id]
-        j = bond.neighbor
+        j, L = bond.neighbor, bond.length
         Δxij = get_vector_diff(position, i, j)
-        ε = bond_stretch[bond_id]
+        l = bond_length[bond_id]
+        ε = (l - L) / L
         ω = bond_active[bond_id] * surface_correction_factor(correction, bond_id)
-        b = ω * params.bc * ε * volume[j] .* Δxij
+        b = ω * params.bc * ε * volume[j] .* Δxij / l
         update_add_vector!(b_int, i, b)
     end
     return nothing
 end
 
 function force_density_point!(storage::BBStorage, system::BondSystem, ::BBMaterial,
                               paramhandler::ParameterHandler, t, Δt, i)
-    (; position, bond_stretch, bond_active, b_int) = storage
+    (; position, bond_length, bond_active, b_int) = storage
     (; bonds, correction, volume) = system
     params_i = get_params(paramhandler, i)
     for bond_id in each_bond_idx(system, i)
         bond = bonds[bond_id]
-        j = bond.neighbor
+        j, L = bond.neighbor, bond.length
         Δxij = get_vector_diff(position, i, j)
-        ε = bond_stretch[bond_id]
+        l = bond_length[bond_id]
+        ε = (l - L) / L
         params_j = get_params(paramhandler, j)
         ω = bond_active[bond_id] * surface_correction_factor(correction, bond_id)
-        b = ω * (params_i.bc + params_j.bc) / 2 * ε * volume[j] .* Δxij
+        b = ω * (params_i.bc + params_j.bc) / 2 * ε * volume[j] .* Δxij / l
         update_add_vector!(b_int, i, b)
     end
     return nothing
 end
+
+function strain_energy_density_point!(storage::AbstractStorage, system::BondSystem,
+                                      ::AbstractBondBasedMaterial,
+                                      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.25 * ωij * bc * ε * ε * L * volume[j]
+    end
+    strain_energy_density[i] = Ψ
+    return nothing
+end
+
+function export_field(::Val{:strain_energy_density}, mat::AbstractBondBasedMaterial,
+                      system::BondSystem, storage::AbstractStorage,
+                      paramsetup::AbstractParameterSetup, t)
+    for i in each_point_idx(system)
+        strain_energy_density_point!(storage, system, mat, paramsetup, i)
+    end
+    return storage.strain_energy_density
+end

src/physics/correspondence.jl

@@ -318,7 +318,7 @@ end
 
 # calculate the von Mises stress from the Cauchy stress tensor just when exporting
 function export_field(::Val{:von_mises_stress}, ::CMaterial, system::BondSystem,
-                      storage::AbstractStorage, t)
+                      storage::AbstractStorage, ::AbstractParameterSetup, t)
     for i in each_point_idx(system)
         σ = get_tensor(storage.cauchy_stress, i)
         storage.von_mises_stress[i] = von_mises_stress(σ)

src/physics/correspondence_rotated.jl

@@ -252,7 +252,7 @@ end
 
 # calculate the von Mises stress from the Cauchy stress tensor just when exporting
 function export_field(::Val{:von_mises_stress}, ::CRMaterial, system::BondSystem,
-                      storage::AbstractStorage, t)
+                      storage::AbstractStorage, ::AbstractParameterSetup, t)
     for i in each_point_idx(system)
         σ = get_tensor(storage.unrotated_stress, i)
         T = rotate_stress(storage::AbstractStorage, σ, i)

src/physics/dh_bond_based.jl

@@ -103,31 +103,24 @@ end
     @pointfield b_ext::Matrix{Float64}
     @pointfield density_matrix::Matrix{Float64}
     @pointfield damage::Vector{Float64}
-    bond_stretch::Vector{Float64}
+    bond_length::Vector{Float64}
     bond_active::Vector{Bool}
     @pointfield n_active_bonds::Vector{Int}
-end
-
-function init_field(::DHBBMaterial, ::AbstractTimeSolver, system::BondSystem, ::Val{:b_int})
-    return zeros(3, get_n_points(system))
-end
-
-function init_field(::DHBBMaterial, ::AbstractTimeSolver, system::BondSystem,
-                    ::Val{:bond_stretch})
-    return zeros(get_n_bonds(system))
+    @pointfield strain_energy_density::Vector{Float64}
 end
 
 function force_density_point!(storage::DHBBStorage, system::BondSystem, ::DHBBMaterial,
                               params::StandardPointParameters, t, Δt, i)
-    (; position, bond_stretch, bond_active, b_int) = storage
+    (; position, bond_length, bond_active, b_int) = storage
     (; bonds, correction, volume) = system
     for bond_id in each_bond_idx(system, i)
         bond = bonds[bond_id]
-        j = bond.neighbor
+        j, L = bond.neighbor, bond.length
         Δxij = get_vector_diff(position, i, j)
-        ε = bond_stretch[bond_id]
+        l = bond_length[bond_id]
+        ε = (l - L) / L
         ω = bond_active[bond_id] * surface_correction_factor(correction, bond_id)
-        b = ω * params.bc * ε .* Δxij
+        b = ω * params.bc * ε .* Δxij / l
         update_add_vector!(b_int, i, b * volume[j])
         update_add_vector!(b_int, j, -b * volume[i])
     end
@@ -136,17 +129,18 @@ end
 
 function force_density_point!(storage::DHBBStorage, system::BondSystem, ::DHBBMaterial,
                               paramhandler::ParameterHandler, t, Δt, i)
-    (; position, bond_stretch, bond_active, b_int) = storage
+    (; position, bond_length, bond_active, b_int) = storage
     (; bonds, correction, volume) = system
     params_i = get_params(paramhandler, i)
     for bond_id in each_bond_idx(system, i)
         bond = bonds[bond_id]
-        j = bond.neighbor
+        j, L = bond.neighbor, bond.length
         Δxij = get_vector_diff(position, i, j)
-        ε = bond_stretch[bond_id]
+        l = bond_length[bond_id]
+        ε = (l - L) / L
         params_j = get_params(paramhandler, j)
         ω = bond_active[bond_id] * surface_correction_factor(correction, bond_id)
-        b = ω * (params_i.bc + params_j.bc) / 2 * ε .* Δxij
+        b = ω * (params_i.bc + params_j.bc) / 2 * ε .* Δxij / l
         update_add_vector!(b_int, i, b * volume[j])
         update_add_vector!(b_int, j, -b * volume[i])
     end

src/physics/ordinary_state_based.jl

@@ -109,6 +109,7 @@ OSBMaterial(; kwargs...) = OSBMaterial{NoCorrection}(; kwargs...)
     bond_length::Vector{Float64}
     bond_active::Vector{Bool}
     @pointfield n_active_bonds::Vector{Int}
+    @pointfield strain_energy_density::Vector{Float64}
 end
 
 function init_field(::OSBMaterial, ::AbstractTimeSolver, system::BondSystem, ::Val{:b_int})
@@ -120,6 +121,11 @@ function init_field(::OSBMaterial, ::AbstractTimeSolver, system::BondSystem,
     return zeros(get_n_bonds(system))
 end
 
+function init_field(::OSBMaterial, ::AbstractTimeSolver, system::BondSystem,
+                    ::Val{:strain_energy_density})
+    return zeros(get_n_loc_points(system))
+end
+
 # Customized calc_failure to save the bond length for force density calculation
 function calc_failure!(storage::OSBStorage, system::BondSystem,
                        ::OSBMaterial, ::CriticalStretch,
@@ -219,3 +225,34 @@ function calc_dilatation(storage::OSBStorage, system::BondSystem, mat::OSBMateri
     end
     return dil
 end
+
+function export_field(::Val{:strain_energy_density}, mat::OSBMaterial, system::BondSystem,
+                      storage::AbstractStorage, paramsetup::AbstractParameterSetup, t)
+    (; bond_active, bond_length, strain_energy_density) = storage
+    (; bonds, correction, volume) = system
+    strain_energy_density .= 0.0
+    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] = Ψ
+    end
+    return strain_energy_density
+end

test/physics/test_strain_energy_density.jl

@@ -0,0 +1,117 @@
+@testsnippet PsiExport begin
+    using Peridynamics.StaticArrays, Peridynamics.Printf
+    mean(x) = sum(x) / length(x)
+    function test_stendens(body, ts, F_a, Ψ_a, tols; testcase="")
+        dh = Peridynamics.threads_data_handler(body, ts, 1)
+        Peridynamics.initialize!(dh, ts)
+        (; mat, system, storage, paramsetup) = dh.chunks[1]
+        field = Val{:strain_energy_density}()
+
+        # no deformation -> no strain energy density
+        storage.position .= system.position
+        Ψ_0 = Peridynamics.export_field(field, mat, system, storage, paramsetup, 0.0)
+        @test all(isapprox.(Ψ_0, 0.0; atol=eps()))
+
+        # apply deformation gradient F_a
+        for i in Peridynamics.each_point_idx(system)
+            Xi = Peridynamics.get_vector(system.position, i)
+            xi = F_a * Xi
+            Peridynamics.update_vector!(storage.position, i, xi)
+        end
+        Peridynamics.calc_force_density!(dh, 0.0, 0.0)
+
+        Ψ_pd = Peridynamics.export_field(field, mat, system, storage, paramsetup, 0.0)
+        Ψ̂_pd = mean(Ψ_pd)
+        ΔΨ = (Ψ_pd .- Ψ_a) ./ Ψ_a
+        ΔΨ_max, ΔΨ_min = maximum(ΔΨ), minimum(ΔΨ)
+
+        # write detailed printed output that shows the calculation results
+        println("-"^20 * "  Ψ TEST  " * "-"^20)
+        if !isempty(testcase)
+            @printf("Test case: %s\n", testcase)
+        end
+        @printf("material: %s\n", nameof(typeof(mat)))
+        @printf("Expected Ψ:        %.6g\n", Ψ_a)
+        @printf("Calc. mean(Ψ):     %.6g\n", Ψ̂_pd)
+        @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]
+        return nothing
+    end
+end
+
+@testitem "Strain energy density export BBMaterial" 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 = BBMaterial{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
+    E = 210e9
+    nu = 0.25
+    λ = 1.05
+    pos, vol = uniform_box(1,1,1,Δx)
+    mat = OSBMaterial{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 = (-1e-5, 1e-5) # should be very accurate for a homogeneous iso. extension
+    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.4, 0.4) # higher errors
+    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 = 1/2 * params.λ * ε^2 + params.μ * ε^2
+    tols = (-0.9, 0.3)
+    test_stendens(body, ts, F_a, Ψ_a, tols; testcase)
+end