diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl index 8e47d8db..37bf6e22 100644 --- a/benchmark/benchmarks.jl +++ b/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)) diff --git a/benchmark/mode_i.jl b/benchmark/mode_i.jl index 37442052..b483d3c6 100644 --- a/benchmark/mode_i.jl +++ b/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 diff --git a/docs/src/api_reference.md b/docs/src/api_reference.md index c797e380..ca87e722 100644 --- a/docs/src/api_reference.md +++ b/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 diff --git a/docs/src/literate/tutorial_wave_interface.jl b/docs/src/literate/tutorial_wave_interface.jl index 4fbefb4e..27445c5d 100644 --- a/docs/src/literate/tutorial_wave_interface.jl +++ b/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) diff --git a/src/Peridynamics.jl b/src/Peridynamics.jl index c3c3a4e8..602a057c 100644 --- a/src/Peridynamics.jl +++ b/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,6 +72,11 @@ 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") @@ -72,10 +84,13 @@ 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 diff --git a/src/auxiliary/nans.jl b/src/auxiliary/nans.jl new file mode 100644 index 00000000..45c16887 --- /dev/null +++ b/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 diff --git a/src/auxiliary/precompile_workload.jl b/src/auxiliary/precompile_workload.jl index af9b9b38..45ad8dfb 100644 --- a/src/auxiliary/precompile_workload.jl +++ b/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 diff --git a/src/auxiliary/static_arrays.jl b/src/auxiliary/static_arrays.jl new file mode 100644 index 00000000..c24a3e47 --- /dev/null +++ b/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 diff --git a/src/core/mpi_body_data_handler.jl b/src/core/mpi_body_data_handler.jl index a9c025cb..8406f086 100644 --- a/src/core/mpi_body_data_handler.jl +++ b/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 diff --git a/src/core/parameters.jl b/src/core/parameters.jl index 7446f78c..2a9a5d28 100644 --- a/src/core/parameters.jl +++ b/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 diff --git a/src/core/threads_body_data_handler.jl b/src/core/threads_body_data_handler.jl index 3f63c099..930e8e8e 100644 --- a/src/core/threads_body_data_handler.jl +++ b/src/core/threads_body_data_handler.jl @@ -157,6 +157,17 @@ function exchange_halo_to_loc!(get_field_function::F, dh::ThreadsBodyDataHandler 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 diff --git a/src/core/threads_multibody_data_handler.jl b/src/core/threads_multibody_data_handler.jl index e7e1d425..8fbf8767 100644 --- a/src/core/threads_multibody_data_handler.jl +++ b/src/core/threads_multibody_data_handler.jl @@ -33,6 +33,17 @@ end @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) diff --git a/src/discretization/body.jl b/src/discretization/body.jl index d4c547d2..35549e28 100644 --- a/src/discretization/body.jl +++ b/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 diff --git a/src/discretization/bond_associated_system.jl b/src/discretization/bond_associated_system.jl new file mode 100644 index 00000000..8a3d0a6f --- /dev/null +++ b/src/discretization/bond_associated_system.jl @@ -0,0 +1,170 @@ + +struct BondAssociatedSystem <: AbstractBondSystem + position::Matrix{Float64} + volume::Vector{Float64} + bonds::Vector{Bond} + n_neighbors::Vector{Int} + bond_ids::Vector{UnitRange{Int}} + intersection_bond_ids::Vector{Vector{Int}} + hood_volume::Vector{Float64} + ba_hood_volume::Vector{Float64} + kernels::Vector{Float64} + chunk_handler::ChunkHandler +end + +function BondAssociatedSystem(body::AbstractBody, pd::PointDecomposition, chunk_id::Int) + check_bond_associated_system_compat(body.mat) + loc_points = pd.decomp[chunk_id] + bonds, n_neighbors = find_bonds(body, loc_points) + bond_ids = find_bond_ids(n_neighbors) + intersection_bond_ids = find_intersection_bond_ids(body, loc_points, bonds, bond_ids) + chunk_handler = get_chunk_handler(bonds, pd, chunk_id) + localize!(bonds, chunk_handler.localizer) + position, volume = get_pos_and_vol_chunk(body, chunk_handler.point_ids) + hood_volume = zeros(get_n_points(chunk_handler)) + ba_hood_volume = zeros(length(bonds)) + kernels = find_kernels(body, chunk_handler, bonds, bond_ids) + bas = BondAssociatedSystem(position, volume, bonds, n_neighbors, bond_ids, + intersection_bond_ids, hood_volume, ba_hood_volume, kernels, + chunk_handler) + return bas +end + +function get_system(body::AbstractBody{Material}, pd::PointDecomposition, + chunk_id::Int) where {Material<:AbstractBondAssociatedSystemMaterial} + return BondAssociatedSystem(body, pd, chunk_id) +end + +@inline function system_type(::AbstractBondAssociatedSystemMaterial) + return BondAssociatedSystem +end + +function check_bond_associated_system_compat(::M) where {M<:AbstractMaterial} + msg = "body with material `$(M)` incompatible to `BondAssociatedSystem`!\n" + msg *= "The material has to be a subtype of `AbstractBondAssociatedSystemMaterial`!\n" + return throw(ArgumentError(msg)) +end + +function check_bond_associated_system_compat(::AbstractBondAssociatedSystemMaterial) + return nothing +end + +function find_intersection_bond_ids(body, loc_points, bonds, bond_ids) + intersection_bond_ids = Vector{Vector{Int}}(undef, length(bonds)) + for (li, i) in enumerate(loc_points) + δb = get_point_param(body, :δb, i) + δb² = δb * δb + bond_ids_of_i = bond_ids[li] + for bond_id in bond_ids_of_i + bond = bonds[bond_id] + j = bond.neighbor + Xj = get_coordinates(body, j) + intersecting_bonds = Vector{Int}() + for (ibond_id, bond_id) in enumerate(bond_ids_of_i) + bond = bonds[bond_id] + jj = bond.neighbor + Xjj = get_coordinates(body, jj) + ΔX = Xj - Xjj + L² = dot(ΔX, ΔX) + if L² < δb² + push!(intersecting_bonds, ibond_id) + end + end + intersection_bond_ids[bond_id] = intersecting_bonds + end + end + return intersection_bond_ids +end + +@inline function each_intersecting_bond_idx(system::BondAssociatedSystem, point_id::Int, + bond_id::Int) + return view(each_bond_idx(system, point_id), system.intersection_bond_ids[bond_id]) +end + +function calc_hood_volumes!(chunk::AbstractBodyChunk{<:BondAssociatedSystem}) + (; system) = chunk + (; volume, bonds, hood_volume, ba_hood_volume) = system + + for i in each_point_idx(chunk) + _hood_volume = volume[i] + for bond_idx in each_bond_idx(system, i) + bond = bonds[bond_idx] + j = bond.neighbor + _hood_volume += 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 = i_bond.neighbor + _ba_hood_volume += 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}, + ::AbstractTimeSolver) + @threads :static for chunk in dh.chunks + calc_hood_volumes!(chunk) + end + @threads :static for chunk_id in eachindex(dh.chunks) + exchange_loc_to_halo!(get_hood_volume, dh, chunk_id) + end + return nothing +end + +function initialize!(dh::AbstractMPIBodyDataHandler{<:BondAssociatedSystem}, + ::AbstractTimeSolver) + calc_hood_volumes!(dh.chunk) + exchange_loc_to_halo!(get_hood_volume, dh) + return nothing +end + +@inline function volume_fraction_factor(system::BondAssociatedSystem, point_idx::Int, + bond_idx::Int) + return system.ba_hood_volume[bond_idx] / system.hood_volume[point_idx] +end + +function req_point_data_fields_fracture(::Type{<:AbstractBondAssociatedSystemMaterial}) + return (:damage, :n_active_bonds) +end + +function req_bond_data_fields_fracture(::Type{<:AbstractBondAssociatedSystemMaterial}) + return (:bond_active,) +end + +function req_data_fields_fracture(::Type{<:AbstractBondAssociatedSystemMaterial}) + return () +end + +function required_point_parameters(::Type{<:AbstractBondAssociatedSystemMaterial}) + return (:δ, :δb, :rho, elasticity_parameters()...) +end + +function get_required_point_parameters(::AbstractBondAssociatedSystemMaterial, + p::Dict{Symbol,Any}) + δ_params = get_horizon(p) + δb_params = get_bond_horizon(p, δ_params.δ) + return (; δ_params..., δb_params..., get_density(p)..., get_elastic_params(p)...) +end + +function get_bond_horizon(p::Dict{Symbol,Any}, δ::Float64) + δb::Float64 = float(get(p, :bond_horizon, δ)) + if δb ≤ 0 + throw(ArgumentError("`bond_horizon` should be larger than zero!\n")) + end + if δb < δ + @warn "a small bond horizon < δ will possibly lead to numerical instabilities!" + end + return (; δb) +end + +function allowed_material_kwargs(::AbstractBondAssociatedSystemMaterial) + return (discretization_kwargs()..., elasticity_kwargs()..., fracture_kwargs()..., + :bond_horizon) +end diff --git a/src/discretization/bond_system.jl b/src/discretization/bond_system.jl index 04b3262e..a341536a 100644 --- a/src/discretization/bond_system.jl +++ b/src/discretization/bond_system.jl @@ -4,12 +4,13 @@ struct Bond fail_permit::Bool end -struct BondSystem{Correction<:AbstractCorrection} <: AbstractSystem +struct BondSystem{Correction<:AbstractCorrection} <: AbstractBondSystem position::Matrix{Float64} volume::Vector{Float64} bonds::Vector{Bond} n_neighbors::Vector{Int} bond_ids::Vector{UnitRange{Int}} + kernels::Vector{Float64} correction::Correction chunk_handler::ChunkHandler end @@ -23,7 +24,8 @@ function BondSystem(body::AbstractBody, pd::PointDecomposition, chunk_id::Int) position, volume = get_pos_and_vol_chunk(body, chunk_handler.point_ids) correction = get_correction(body.mat, chunk_handler.n_loc_points, length(chunk_handler.point_ids), length(bonds)) - system = BondSystem(position, volume, bonds, n_neighbors, bond_ids, correction, + kernels = find_kernels(body, chunk_handler, bonds, bond_ids) + system = BondSystem(position, volume, bonds, n_neighbors, bond_ids, kernels, correction, chunk_handler) return system end @@ -146,6 +148,29 @@ function get_chunk_handler(bonds::Vector{Bond}, pd::PointDecomposition, chunk_id localizer) end +function find_kernels(body::AbstractBody, chunk_handler::ChunkHandler, bonds::Vector{Bond}, + bond_ids::Vector{UnitRange{Int}}) + hasproperty(body.mat, :kernel) || return Vector{Float64}() + kernels = zeros(length(bonds)) + for i in each_point_idx(chunk_handler) + params = get_point_param(body, i) + for bond_id in bond_ids[i] + bond = bonds[bond_id] + kernels[bond_id] = get_kernel(body.mat, params, bond.length) + end + end + return kernels +end + +function get_kernel(mat::AbstractMaterial, params::AbstractPointParameters, L) + ω = mat.kernel(params.δ, L) + return ω +end + +@inline function kernel(system::AbstractBondSystem, bond_id::Int) + return system.kernels[bond_id] +end + function find_halo_points(bonds::Vector{Bond}, loc_points::UnitRange{Int}) halo_points = Vector{Int}() for bond in bonds @@ -157,7 +182,7 @@ function find_halo_points(bonds::Vector{Bond}, loc_points::UnitRange{Int}) return halo_points end -@inline each_bond_idx(bd::BondSystem, point_id::Int) = bd.bond_ids[point_id] +@inline each_bond_idx(system::AbstractBondSystem, point_id::Int) = system.bond_ids[point_id] function localize!(bonds::Vector{Bond}, localizer::Dict{Int,Int}) for i in eachindex(bonds) @@ -167,8 +192,8 @@ function localize!(bonds::Vector{Bond}, localizer::Dict{Int,Int}) return nothing end -function break_bonds!(storage::AbstractStorage, system::BondSystem, set_a::Vector{Int}, - set_b::Vector{Int}) +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) @@ -187,26 +212,28 @@ function break_bonds!(storage::AbstractStorage, system::BondSystem, set_a::Vecto return nothing end -function calc_timestep_point(bd::BondSystem, params::AbstractPointParameters, point_id::Int) +function calc_timestep_point(system::AbstractBondSystem, params::AbstractPointParameters, + point_id::Int) dtsum = 0.0 - for bond_id in each_bond_idx(bd, point_id) - bond = bd.bonds[bond_id] - dtsum += bd.volume[bond.neighbor] * params.bc / bond.length + for bond_id in each_bond_idx(system, point_id) + bond = system.bonds[bond_id] + dtsum += system.volume[bond.neighbor] * params.bc / bond.length end return sqrt(2 * params.rho / dtsum) end -function calc_force_density!(chunk::AbstractBodyChunk{S,M}) where {S<:BondSystem,M} +function calc_force_density!(chunk::AbstractBodyChunk{<:AbstractBondSystem}, t, Δt) (; system, mat, paramsetup, storage) = chunk storage.b_int .= 0 storage.n_active_bonds .= 0 for point_id in each_point_idx(chunk) - force_density_point!(storage, system, mat, paramsetup, point_id) + force_density_point!(storage, system, mat, paramsetup, t, Δt, point_id) end + nancheck(chunk, t) return nothing end -@inline function calc_damage!(chunk::AbstractBodyChunk{S,M}) where {S<:BondSystem,M} +@inline function calc_damage!(chunk::AbstractBodyChunk{<:AbstractBondSystem}) (; n_neighbors) = chunk.system (; n_active_bonds, damage) = chunk.storage for point_id in each_point_idx(chunk) @@ -215,7 +242,7 @@ end return nothing end -@inline function stretch_based_failure!(storage::AbstractStorage, ::BondSystem, +@inline function stretch_based_failure!(storage::AbstractStorage, ::AbstractBondSystem, bond::Bond, params::AbstractPointParameters, ε::Float64, i::Int, bond_id::Int) if ε > params.εc && bond.fail_permit @@ -229,8 +256,8 @@ end return storage.bond_active[bond_id] end -function log_system(::Type{B}, options::AbstractJobOptions, - dh::AbstractDataHandler) where {B<:BondSystem} +function log_system(::Type{System}, options::AbstractJobOptions, + dh::AbstractDataHandler) where {System<:AbstractBondSystem} n_bonds = calc_n_bonds(dh) msg = "BOND SYSTEM" body_name = string(get_body_name(dh)) @@ -254,15 +281,15 @@ function calc_n_bonds(dh::AbstractMPIBodyDataHandler) return n_bonds end -function init_field_system(system::BondSystem, ::Val{:bond_active}) +function init_field_system(system::AbstractBondSystem, ::Val{:bond_active}) return ones(Bool, get_n_bonds(system)) end -function init_field_system(system::BondSystem, ::Val{:n_active_bonds}) +function init_field_system(system::AbstractBondSystem, ::Val{:n_active_bonds}) return copy(system.n_neighbors) end -function init_field_system(system::BondSystem, ::Val{:damage}) +function init_field_system(system::AbstractBondSystem, ::Val{:damage}) return zeros(get_n_loc_points(system)) end @@ -290,4 +317,4 @@ function allowed_material_kwargs(::AbstractBondSystemMaterial) return (discretization_kwargs()..., elasticity_kwargs()..., fracture_kwargs()...) end -@inline get_n_bonds(system::BondSystem) = length(system.bonds) +@inline get_n_bonds(system::AbstractBondSystem) = length(system.bonds) diff --git a/src/discretization/interaction_system.jl b/src/discretization/interaction_system.jl index 7a8fc05e..9338d695 100644 --- a/src/discretization/interaction_system.jl +++ b/src/discretization/interaction_system.jl @@ -354,17 +354,17 @@ function calc_timestep_point(system::InteractionSystem, params::AbstractPointPar return sqrt(2 * params.rho / dtsum) end -function calc_force_density!(chunk::AbstractBodyChunk{S,M}) where {S<:InteractionSystem,M} +function calc_force_density!(chunk::AbstractBodyChunk{<:InteractionSystem}, t, Δt) (; system, mat, paramsetup, storage) = chunk storage.b_int .= 0 storage.n_active_one_nis .= 0 for point_id in each_point_idx(chunk) - force_density_point!(storage, system, mat, paramsetup, point_id) + force_density_point!(storage, system, mat, paramsetup, t, Δt, point_id) end return nothing end -@inline function calc_damage!(chunk::AbstractBodyChunk{S,M}) where {S<:InteractionSystem,M} +@inline function calc_damage!(chunk::AbstractBodyChunk{<:InteractionSystem}) (; n_one_nis) = chunk.system (; n_active_one_nis, damage) = chunk.storage for point_id in each_point_idx(chunk) diff --git a/src/discretization/kernels.jl b/src/discretization/kernels.jl new file mode 100644 index 00000000..2c977a4b --- /dev/null +++ b/src/discretization/kernels.jl @@ -0,0 +1,41 @@ +@doc raw""" + linear_kernel(δ, L) + +A linear kernel function ``\omega`` (also called influence function) used for weighting the +bonds in a family. The kernel function is defined as +```math +\omega = \frac{\delta}{L} \; , +``` +with the horizon ``\delta`` and the initial bond length ``L``. +""" +@inline linear_kernel(δ, L) = δ / L + +@doc raw""" + cubic_b_spline_kernel(δ, 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 &= \left\{ + \begin{array}{ll} + \frac{2}{3} - 4 \xi^2 + 4 \xi^3 & \quad \text{if} \; 0 < \xi \leq 0.5 \; , \\[3pt] + \frac{4}{3} - 4 \xi + 4 \xi^2 - \frac{4}{3} \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``. +""" +@inline function cubic_b_spline_kernel(δ, L) + ξ = L / δ + if 0 < ξ ≤ 0.5 + return 2/3 - 4 * ξ^2 + 4 * ξ^3 + elseif 0.5 < ξ ≤ 1 + return 4/3 - 4 * ξ + 4 * ξ^2 - 4/3 * ξ^3 + else + return 0 + end +end diff --git a/src/discretization/zem_stabilization.jl b/src/discretization/zem_stabilization.jl new file mode 100644 index 00000000..7be22cd7 --- /dev/null +++ b/src/discretization/zem_stabilization.jl @@ -0,0 +1,21 @@ +function get_correction(mat::AbstractBondSystemMaterial{<:AbstractZEMStabilization}, + ::Int, ::Int, ::Int) + return mat.zem_stabilization +end + +""" + ZEMSilling(; Cs) + +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 + function ZEMSilling(; Cs::Real=100.0) + return new(Cs) + end +end diff --git a/src/physics/ba_correspondence.jl b/src/physics/ba_correspondence.jl new file mode 100644 index 00000000..098f9491 --- /dev/null +++ b/src/physics/ba_correspondence.jl @@ -0,0 +1,215 @@ +""" + BACMaterial(; kernel, model, maxdmg) + +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.85`) + +# Examples + +```julia-repl +julia> mat = BACMaterial() +CMaterial(maxdmg=0.95, zem_fac=ZEMSilling()) +``` + +--- + +```julia +BACMaterial{CM,K} +``` + +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. + +# Allowed material parameters +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 +- `E::Float64`: Young's modulus +- `nu::Float64`: Poisson's ratio +- `Gc::Float64`: Critical energy release rate +- `epsilon_c::Float64`: Critical strain + +# Allowed export fields +When specifying the `fields` keyword of [`Job`](@ref) for a [`Body`](@ref) with +`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 +- `velocity_half::Matrix{Float64}`: Velocity parameter for Verlet time solver +- `acceleration::Matrix{Float64}`: Acceleration of each point +- `b_int::Matrix{Float64}`: Internal force density of each point +- `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 +""" +struct BACMaterial{CM,K} <: AbstractBondAssociatedSystemMaterial + kernel::K + constitutive_model::CM + maxdmg::Float64 + function BACMaterial(kernel::K, cm::CM, maxdmg::Real) where {K,CM} + return new{CM,K}(kernel, cm, maxdmg) + end +end + +function Base.show(io::IO, @nospecialize(mat::BACMaterial)) + print(io, typeof(mat)) + print(io, msg_fields_in_brackets(mat, (:maxdmg,))) + return nothing +end + +function BACMaterial(; kernel::Function=linear_kernel, + model::AbstractConstitutiveModel=LinearElastic(), + maxdmg::Real=0.85) + return BACMaterial(kernel, model, maxdmg) +end + +struct BACPointParameters <: AbstractPointParameters + δ::Float64 + δb::Float64 + rho::Float64 + E::Float64 + nu::Float64 + G::Float64 + K::Float64 + λ::Float64 + μ::Float64 + Gc::Float64 + εc::Float64 + bc::Float64 +end + +function BACPointParameters(mat::BACMaterial, p::Dict{Symbol,Any}) + (; δ, δb, rho, E, nu, G, K, λ, μ) = get_required_point_parameters(mat, p) + (; Gc, εc) = get_frac_params(p, δ, K) + bc = 18 * K / (π * δ^4) # bond constant + return BACPointParameters(δ, δb, rho, E, nu, G, K, λ, μ, Gc, εc, bc) +end + +@params BACMaterial BACPointParameters + +@storage BACMaterial struct BACStorage + @lthfield position::Matrix{Float64} + @pointfield displacement::Matrix{Float64} + @pointfield velocity::Matrix{Float64} + @pointfield velocity_half::Matrix{Float64} + @pointfield velocity_half_old::Matrix{Float64} + @pointfield acceleration::Matrix{Float64} + @htlfield b_int::Matrix{Float64} + @pointfield b_int_old::Matrix{Float64} + @pointfield b_ext::Matrix{Float64} + @pointfield density_matrix::Matrix{Float64} + @pointfield damage::Vector{Float64} + bond_active::Vector{Bool} + @pointfield n_active_bonds::Vector{Int} + @pointfield stress::Matrix{Float64} + @pointfield von_mises_stress::Vector{Float64} +end + +function init_field(::BACMaterial, ::AbstractTimeSolver, system::BondAssociatedSystem, + ::Val{:b_int}) + return zeros(3, get_n_points(system)) +end + +function init_field(::BACMaterial, ::AbstractTimeSolver, system::BondAssociatedSystem, + ::Val{:stress}) + return zeros(9, get_n_loc_points(system)) +end + +function init_field(::BACMaterial, ::AbstractTimeSolver, system::BondAssociatedSystem, + ::Val{:von_mises_stress}) + 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) + force_density_bond!(storage, system, mat, params, t, Δt, i, bond_idx) + end + return nothing +end + +function force_density_bond!(storage::BACStorage, system::BondAssociatedSystem, + mat::BACMaterial, params::BACPointParameters, t, Δt, i, + bond_idx) + defgrad_res = calc_deformation_gradient(storage, system, mat, params, i, bond_idx) + (; F) = defgrad_res + if containsnan(F) || storage.damage[i] > mat.maxdmg + storage.bond_active[bond_idx] = false + return nothing + end + PKinv = calc_first_piola_kirchhoff!(storage, mat, params, defgrad_res, Δt, i, bond_idx) + + bond = system.bonds[bond_idx] + 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_idx) + + ωij = kernel(system, bond_idx) * storage.bond_active[bond_idx] + ϕi = volume_fraction_factor(system, i, bond_idx) + tij = ϕi * ωij * PKinv * ΔXij + update_add_b_int!(storage, i, tij .* system.volume[j]) + update_add_b_int!(storage, j, -tij .* system.volume[i]) + return nothing +end + +function calc_deformation_gradient(storage::BACStorage, system::BondAssociatedSystem, + mat::BACMaterial, params::BACPointParameters, i, + bond_idx) + (; bonds, volume) = system + (; bond_active) = storage + K = zero(SMatrix{3,3,Float64,9}) + _F = zero(SMatrix{3,3,Float64,9}) + for bond_id in each_intersecting_bond_idx(system, i, bond_idx) + bond = bonds[bond_id] + j = bond.neighbor + ΔXij = get_diff(system.position, i, j) + Δxij = get_diff(storage.position, i, j) + ωijV = kernel(system, bond_id) * volume[j] + ωijωDV = ωijV * bond_active[bond_id] + K += ωijV * (ΔXij * ΔXij') + _F += ωijωDV * (Δxij * ΔXij') + end + Kinv = inv(K) + F = _F * Kinv + return (; F, Kinv) +end + +function calc_first_piola_kirchhoff!(storage::BACStorage, mat::BACMaterial, + params::BACPointParameters, defgrad_res, Δt, i, + bond_idx) + (; F, Kinv) = defgrad_res + P = first_piola_kirchhoff(mat.constitutive_model, storage, params, F) + PKinv = P * Kinv + return PKinv +end diff --git a/src/physics/bond_based.jl b/src/physics/bond_based.jl index e872ac6b..2bb1e003 100644 --- a/src/physics/bond_based.jl +++ b/src/physics/bond_based.jl @@ -108,7 +108,7 @@ end end function force_density_point!(storage::BBStorage, system::BondSystem, ::BBMaterial, - params::BBPointParameters, i::Int) + params::BBPointParameters, t, Δt, i) for bond_id in each_bond_idx(system, i) bond = system.bonds[bond_id] j, L = bond.neighbor, bond.length @@ -125,7 +125,7 @@ function force_density_point!(storage::BBStorage, system::BondSystem, ::BBMateri end function force_density_point!(storage::BBStorage, system::BondSystem, ::BBMaterial, - paramhandler::ParameterHandler, i::Int) + paramhandler::ParameterHandler, t, Δt, i) params_i = get_params(paramhandler, i) for bond_id in each_bond_idx(system, i) bond = system.bonds[bond_id] diff --git a/src/physics/constitutive_models.jl b/src/physics/constitutive_models.jl new file mode 100644 index 00000000..100754f2 --- /dev/null +++ b/src/physics/constitutive_models.jl @@ -0,0 +1,134 @@ +@doc raw""" + LinearElastic + +Linear elastic constitutive model that can be specified when using a [`CMaterial`](@ref) and +[`BACMaterial`](@ref). +The first Piola-Kirchhoff stress ``\boldsymbol{P}`` is given by +```math +\boldsymbol{P} = \mathbb{C} : \boldsymbol{E} \; , +``` +with the elastic stiffness tensor ``\mathbb{C}`` and the Green-Lagrange strain tensor +``\boldsymbol{E}`` with +```math +\boldsymbol{E} = \frac{1}{2} \left( \boldsymbol{F}^{\top} \boldsymbol{F} - \boldsymbol{I} + \right) \; . +``` +""" +struct LinearElastic <: AbstractConstitutiveModel end + +function first_piola_kirchhoff(::LinearElastic, storage::AbstractStorage, + params::AbstractPointParameters, F::SMatrix{3,3,T,9}) where T + 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.μ) + Pvoigt = Cvoigt * Evoigt + P = SMatrix{3,3,Float64,9}(Pvoigt[1], Pvoigt[6], Pvoigt[5], + Pvoigt[6], Pvoigt[2], Pvoigt[4], + Pvoigt[5], Pvoigt[4], Pvoigt[3]) + return P +end + +function get_hooke_matrix(nu, λ, μ) + a = (1 - nu) * λ / nu + CVoigt = SMatrix{6,6,Float64,36}(a, λ, λ, 0, 0, 0, λ, a, λ, 0, 0, 0, λ, λ, a, 0, 0, 0, + 0, 0, 0, μ, 0, 0, 0, 0, 0, 0, μ, 0, 0, 0, 0, 0, 0, μ) + return CVoigt +end + +@doc raw""" + NeoHooke + +Neo-Hookean constitutive model that can be specified when using a [`CMaterial`](@ref) and +[`BACMaterial`](@ref). +The first Piola-Kirchhoff stress ``\boldsymbol{P}`` is given by +```math +\begin{aligned} +\boldsymbol{C} &= \boldsymbol{F}^{\top} \boldsymbol{F} \; , \\ +\boldsymbol{S} &= \mu \left( \boldsymbol{I} - \boldsymbol{C}^{-1} \right) + + \lambda \log(J) \boldsymbol{C}^{-1} \; , \\ +\boldsymbol{P} &= \boldsymbol{F} \, \boldsymbol{S} \; , +\end{aligned} +``` +with the deformation gradient ``\boldsymbol{F}``, the right Cauchy-Green deformation tensor +``\boldsymbol{C}``, the Jacobian ``J = \mathrm{det}(\boldsymbol{F})``, the second +Piola-Kirchhoff stress ``\boldsymbol{S}``, and the first and second Lamé parameters +``\lambda`` and ``\mu``. +""" +struct NeoHooke <: AbstractConstitutiveModel end + +function first_piola_kirchhoff(::NeoHooke, storage::AbstractStorage, + params::AbstractPointParameters, F::SMatrix{3,3,T,9}) where T + J = det(F) + Cinv = inv(F' * F) + S = params.μ * (I - Cinv) + params.λ * log(J) * Cinv + P = F * S + return P +end + +@doc raw""" + MooneyRivlin + +Mooney-Rivlin constitutive model that can be specified when using a [`CMaterial`](@ref) and +[`BACMaterial`](@ref). +The first Piola-Kirchhoff stress ``\boldsymbol{P}`` is given by +```math +\begin{aligned} +\boldsymbol{C} &= \boldsymbol{F}^{\top} \boldsymbol{F} \; , \\ +\boldsymbol{S} &= G \left( \boldsymbol{I} - \frac{1}{3} \mathrm{tr}(\boldsymbol{C}) + \boldsymbol{C}^{-1} \right) \cdot J^{-\frac{2}{3}} + + \frac{K}{4} \left( J^2 - J^{-2} \right) \boldsymbol{C}^{-1} \; , \\ +\boldsymbol{P} &= \boldsymbol{F} \, \boldsymbol{S} \; , +\end{aligned} +``` +with the deformation gradient ``\boldsymbol{F}``, the right Cauchy-Green deformation tensor +``\boldsymbol{C}``, the Jacobian ``J = \mathrm{det}(\boldsymbol{F})``, the second +Piola-Kirchhoff stress ``\boldsymbol{S}``, the shear modulus ``G``, and the +bulk modulus ``K``. + +# Error handling +If the Jacobian ``J`` is smaller than the machine precision `eps()` or a `NaN`, the first +Piola-Kirchhoff stress tensor is defined as ``\boldsymbol{P} = \boldsymbol{0}``. +""" +struct MooneyRivlin <: AbstractConstitutiveModel end + +function first_piola_kirchhoff(::MooneyRivlin, 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}) + C = F' * F + Cinv = inv(C) + S = params.G .* (I - 1 / 3 .* tr(C) .* Cinv) .* J^(-2 / 3) .+ + params.K / 4 .* (J^2 - J^(-2)) .* Cinv + P = F * S + return P +end + +@doc raw""" + SaintVenantKirchhoff + +Saint-Venant-Kirchhoff constitutive model that can be specified when using a +[`CMaterial`](@ref) and [`BACMaterial`](@ref). +The first Piola-Kirchhoff stress ``\boldsymbol{P}`` is given by +```math +\begin{aligned} +\boldsymbol{E} &= \frac{1}{2} \left( \boldsymbol{F}^{\top} \boldsymbol{F} - \boldsymbol{I} + \right) \; , \\ +\boldsymbol{S} &= \lambda \, \mathrm{tr}(\boldsymbol{E}) \, \boldsymbol{I} + + 2 \mu \boldsymbol{E} \; , \\ +\boldsymbol{P} &= \boldsymbol{F} \, \boldsymbol{S} \; , +\end{aligned} +``` +with the deformation gradient ``\boldsymbol{F}``, the Green-Lagrange strain tensor +``\boldsymbol{E}``, the second Piola-Kirchhoff stress ``\boldsymbol{S}``, and the first and +second Lamé parameters ``\lambda`` and ``\mu``. +""" +struct SaintVenantKirchhoff <: AbstractConstitutiveModel end + +function first_piola_kirchhoff(::SaintVenantKirchhoff, storage::AbstractStorage, + params::AbstractPointParameters, F::SMatrix{3,3,T,9}) where T + E = 0.5 .* (F' * F - I) + S = params.λ * tr(E) * I + 2 * params.μ * E + P = F * S + return P +end diff --git a/src/physics/continuum_kinematics_inspired.jl b/src/physics/continuum_kinematics_inspired.jl index c46fa768..83a9c1cc 100644 --- a/src/physics/continuum_kinematics_inspired.jl +++ b/src/physics/continuum_kinematics_inspired.jl @@ -97,15 +97,15 @@ end end function force_density_point!(storage::CKIStorage, system::InteractionSystem, - mat::CKIMaterial, params::AbstractParameterSetup, i::Int) - force_density_point_one_ni!(storage, system, mat, params, i) - has_two_nis(params) && force_density_point_two_ni!(storage, system, mat, params, i) - has_three_nis(params) && force_density_point_three_ni!(storage, system, mat, params, i) + mat::CKIMaterial, params::AbstractParameterSetup, t, Δt, i) + force_density_point_one_ni!(storage, system, mat, params, t, Δt, i) + has_two_nis(params) && force_density_point_two_ni!(storage, system, mat, params, t, Δt, i) + has_three_nis(params) && force_density_point_three_ni!(storage, system, mat, params, t, Δt, i) return nothing end function force_density_point_one_ni!(storage::CKIStorage, system::InteractionSystem, - ::CKIMaterial, params::CKIPointParameters, i::Int) + ::CKIMaterial, params::CKIPointParameters, t, Δt, i) for one_ni_id in each_one_ni_idx(system, i) one_ni = system.one_nis[one_ni_id] j, L = one_ni.neighbor, one_ni.length @@ -121,7 +121,7 @@ function force_density_point_one_ni!(storage::CKIStorage, system::InteractionSys end function force_density_point_one_ni!(storage::CKIStorage, system::InteractionSystem, - ::CKIMaterial, paramhandler::ParameterHandler, i::Int) + ::CKIMaterial, paramhandler::ParameterHandler, t, Δt, i) params_i = get_params(paramhandler, i) for one_ni_id in each_one_ni_idx(system, i) one_ni = system.one_nis[one_ni_id] @@ -139,7 +139,7 @@ function force_density_point_one_ni!(storage::CKIStorage, system::InteractionSys end function force_density_point_two_ni!(storage::CKIStorage, system::InteractionSystem, - ::CKIMaterial, params::CKIPointParameters, i::Int) + ::CKIMaterial, params::CKIPointParameters, t, Δt, i) for two_ni_id in each_two_ni_idx(system, i) two_ni = system.two_nis[two_ni_id] oni_j_id, oni_k_id, surface_ref = two_ni.oni_j, two_ni.oni_k, two_ni.surface @@ -176,7 +176,7 @@ function force_density_point_two_ni!(storage::CKIStorage, system::InteractionSys end function force_density_point_two_ni!(storage::CKIStorage, system::InteractionSystem, - ::CKIMaterial, paramhandler::ParameterHandler, i::Int) + ::CKIMaterial, paramhandler::ParameterHandler, t, Δt, i) params_i = get_params(paramhandler, i) for two_ni_id in each_two_ni_idx(system, i) two_ni = system.two_nis[two_ni_id] @@ -216,7 +216,7 @@ function force_density_point_two_ni!(storage::CKIStorage, system::InteractionSys end function force_density_point_three_ni!(storage::CKIStorage, system::InteractionSystem, - ::CKIMaterial, params::CKIPointParameters, i::Int) + ::CKIMaterial, params::CKIPointParameters, t, Δt, i) for three_ni_id in each_three_ni_idx(system, i) three_ni = system.three_nis[three_ni_id] oni_j_id = three_ni.oni_j @@ -276,7 +276,7 @@ end function force_density_point_three_ni!(storage::CKIStorage, system::InteractionSystem, ::CKIMaterial, paramhandler::ParameterHandler, - i::Int) + t, Δt, i) params_i = get_params(paramhandler, i) for three_ni_id in each_three_ni_idx(system, i) three_ni = system.three_nis[three_ni_id] diff --git a/src/physics/correspondence.jl b/src/physics/correspondence.jl index 8708477b..33a75458 100644 --- a/src/physics/correspondence.jl +++ b/src/physics/correspondence.jl @@ -1,52 +1,69 @@ """ - NOSBMaterial(; maxdmg, maxjacobi, corr) + CMaterial(; kernel, model, zem, 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. # Keywords +- `kernel::Function`: Kernel function used for weighting the interactions between points. \\ + (default: `linear_kernel`) \\ + The following kernels can be used: + - [`linear_kernel`](@ref) + - [`cubic_b_spline_kernel`](@ref) +- `model::AbstractConstitutiveModel`: Constitutive model defining the material behavior. \\ + (default: `LinearElastic()`) \\ + The following models can be used: + - [`LinearElastic`](@ref) + - [`NeoHooke`](@ref) + - [`MooneyRivlin`](@ref) + - [`SaintVenantKirchhoff`](@ref) +- `zem::AbstractZEMStabilization`: Zero-energy mode stabilization. The + stabilization algorithm of Silling (2017) is used as default. \\ + (default: `ZEMSilling()`) - `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`) -- `maxjacobi::Float64`: Maximum value of the Jacobi determinant. If this value is exceeded, - all bonds of that point are broken. - (default: `1.03`) -- `corr::Float64`: Correction factor used for zero-energy mode stabilization. The - stabilization algorithm of Silling (2017) is used. - (default: `100.0`) + possibly contain `NaN` values. \\ + (default: `0.85`) !!! note "Stability of fracture simulations" - This formulation is known to be not suitable for fracture simultations without + This formulation is known to be not suitable for fracture simulations without stabilization of the zero-energy modes. Therefore be careful when doing fracture - simulations and try out different paremeters for `maxdmg`, `maxjacobi`, and `corr`. + simulations and try out different parameters for `maxdmg` and `zem`. # Examples ```julia-repl -julia> mat = NOSBMaterial() -NOSBMaterial(maxdmg=0.95, maxjacobi=1.03, corr=100.0) +julia> mat = CMaterial() +CMaterial(maxdmg=0.85, zem=ZEMSilling()) ``` --- ```julia -NOSBMaterial +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 -- `maxdmg::Float64`: Maximum value of damage a point is allowed to obtain. See the - constructor docs for more informations. -- `maxjacobi::Float64`: Maximum value of the Jacobi determinant. See the constructor docs +- `kernel::Function`: Kernel function used for weighting the interactions between points. + See the constructor docs for more informations. +- `model::AbstractConstitutiveModel`: Constitutive model defining the material behavior. See + the constructor docs for more informations. +- `zem::AbstractZEMStabilization`: Zero-energy mode stabilization. See the constructor docs for more informations. -- `corr::Float64`: Correction factor used for zero-energy mode stabilization. See the +- `maxdmg::Float64`: Maximum value of damage a point is allowed to obtain. See the constructor docs for more informations. # Allowed material parameters -When using [`material!`](@ref) on a [`Body`](@ref) with `NOSBMaterial`, then the following +When using [`material!`](@ref) on a [`Body`](@ref) with `CMaterial`, then the following parameters are allowed: - `horizon::Float64`: Radius of point interactions - `rho::Float64`: Density @@ -57,7 +74,7 @@ parameters are allowed: # Allowed export fields When specifying the `fields` keyword of [`Job`](@ref) for a [`Body`](@ref) with -`NOSBMaterial`, the following fields are allowed: +`CMaterial`, 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 @@ -67,20 +84,32 @@ 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 +- `von_mises_stress::Vector{Float64}`: Von Mises stress of each point """ -Base.@kwdef struct NOSBMaterial <: AbstractBondSystemMaterial{NoCorrection} - maxdmg::Float64 = 0.95 - maxjacobi::Float64 = 1.03 - corr::Float64 = 100.0 +struct CMaterial{CM,ZEM,K} <: AbstractCorrespondenceMaterial{CM,ZEM} + kernel::K + constitutive_model::CM + zem_stabilization::ZEM + maxdmg::Float64 + function CMaterial(kernel::K, cm::CM, zem::ZEM, maxdmg::Real) where {CM,ZEM,K} + return new{CM,ZEM,K}(kernel, cm, zem, maxdmg) + end end -function Base.show(io::IO, @nospecialize(mat::NOSBMaterial)) +function Base.show(io::IO, @nospecialize(mat::CMaterial)) print(io, typeof(mat)) - print(io, msg_fields_in_brackets(mat)) + print(io, msg_fields_in_brackets(mat, (:maxdmg,))) return nothing end -struct NOSBPointParameters <: AbstractPointParameters +function CMaterial(; kernel::Function=linear_kernel, + model::AbstractConstitutiveModel=LinearElastic(), + zem::AbstractZEMStabilization=ZEMSilling(), maxdmg::Real=0.85) + return CMaterial(kernel, model, zem, maxdmg) +end + +struct CPointParameters <: AbstractPointParameters δ::Float64 rho::Float64 E::Float64 @@ -94,16 +123,16 @@ struct NOSBPointParameters <: AbstractPointParameters bc::Float64 end -function NOSBPointParameters(mat::NOSBMaterial, p::Dict{Symbol,Any}) +function CPointParameters(mat::CMaterial, p::Dict{Symbol,Any}) (; δ, rho, E, nu, G, K, λ, μ) = get_required_point_parameters(mat, p) (; Gc, εc) = get_frac_params(p, δ, K) bc = 18 * K / (π * δ^4) # bond constant - return NOSBPointParameters(δ, rho, E, nu, G, K, λ, μ, Gc, εc, bc) + return CPointParameters(δ, rho, E, nu, G, K, λ, μ, Gc, εc, bc) end -@params NOSBMaterial NOSBPointParameters +@params CMaterial CPointParameters -@storage NOSBMaterial struct NOSBStorage <: AbstractStorage +@storage CMaterial struct CStorage @lthfield position::Matrix{Float64} @pointfield displacement::Matrix{Float64} @pointfield velocity::Matrix{Float64} @@ -117,34 +146,86 @@ end @pointfield damage::Vector{Float64} bond_active::Vector{Bool} @pointfield n_active_bonds::Vector{Int} + @pointfield stress::Matrix{Float64} + @pointfield von_mises_stress::Vector{Float64} end -function init_field(::NOSBMaterial, ::AbstractTimeSolver, system::BondSystem, ::Val{:b_int}) +function init_field(::CMaterial, ::AbstractTimeSolver, system::BondSystem, ::Val{:b_int}) return zeros(3, get_n_points(system)) end -function force_density_point!(storage::NOSBStorage, system::BondSystem, mat::NOSBMaterial, - paramhandler::AbstractParameterHandler, i::Int) +function init_field(::CMaterial, ::AbstractTimeSolver, system::BondSystem, ::Val{:stress}) + return zeros(9, get_n_loc_points(system)) +end + +function init_field(::CMaterial, ::AbstractTimeSolver, system::BondSystem, + ::Val{:von_mises_stress}) + return zeros(get_n_loc_points(system)) +end + +function force_density_point!(storage::AbstractStorage, system::AbstractSystem, + mat::AbstractCorrespondenceMaterial, + paramhandler::AbstractParameterHandler, t, Δt, i) params = get_params(paramhandler, i) - force_density_point!(storage, system, mat, params, i) + force_density_point!(storage, system, mat, params, t, Δt, i) return nothing end -function force_density_point!(storage::NOSBStorage, system::BondSystem, mat::NOSBMaterial, - params::NOSBPointParameters, i::Int) - F, Kinv, ω0 = calc_deformation_gradient(storage, system, mat, params, i) - if storage.damage[i] > mat.maxdmg || containsnan(F) - kill_point!(storage, system, i) - return nothing - end - P = calc_first_piola_stress(F, mat, params) - if iszero(P) || containsnan(P) - kill_point!(storage, system, i) - return nothing +function force_density_point!(storage::AbstractStorage, system::AbstractSystem, + mat::AbstractCorrespondenceMaterial, + params::AbstractPointParameters, t, Δt, 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) + zem = mat.zem_stabilization + c_force_density!(storage, system, mat, params, zem, PKinv, defgrad_res, i) + return nothing +end + +function calc_deformation_gradient(storage::CStorage, system::BondSystem, ::CMaterial, + ::CPointParameters, i) + (; bonds, volume) = system + (; bond_active) = storage + K = zero(SMatrix{3,3,Float64,9}) + _F = zero(SMatrix{3,3,Float64,9}) + ω0 = 0.0 + for bond_id in each_bond_idx(system, i) + bond = bonds[bond_id] + j = bond.neighbor + ΔXij = get_diff(system.position, i, j) + Δxij = get_diff(storage.position, i, j) + ωij = kernel(system, bond_id) * bond_active[bond_id] + ω0 += ωij + temp = ωij * volume[j] + ΔXijt = ΔXij' + K += temp * (ΔXij * ΔXijt) + _F += temp * (Δxij * ΔXijt) end + Kinv = inv(K) + F = _F * Kinv + return (; F, Kinv, ω0) +end + +function calc_first_piola_kirchhoff!(storage::CStorage, mat::CMaterial, + params::CPointParameters, defgrad_res, Δt, i) + (; F, Kinv) = defgrad_res + P = first_piola_kirchhoff(mat.constitutive_model, storage, params, F) PKinv = P * Kinv + σ = cauchy_stress(P, F) + update_tensor!(storage.stress, i, σ) + storage.von_mises_stress[i] = von_mises_stress(σ) + return PKinv +end + +function c_force_density!(storage::AbstractStorage, system::AbstractSystem, + ::AbstractCorrespondenceMaterial, params::AbstractPointParameters, + zem_correction::ZEMSilling, PKinv, defgrad_res, i) + (; bonds, volume) = system + (; bond_active) = storage + (; F, ω0) = defgrad_res + (; Cs) = zem_correction for bond_id in each_bond_idx(system, i) - bond = system.bonds[bond_id] + bond = bonds[bond_id] j, L = bond.neighbor, bond.length ΔXij = get_coordinates_diff(system, i, j) Δxij = get_coordinates_diff(storage, i, j) @@ -153,68 +234,26 @@ function force_density_point!(storage::NOSBStorage, system::BondSystem, mat::NOS stretch_based_failure!(storage, system, bond, params, ε, i, bond_id) # stabilization - ωij = influence_function(mat, params, L) * storage.bond_active[bond_id] - Tij = mat.corr .* params.bc * ωij / ω0 .* (Δxij .- F * ΔXij) + ωij = kernel(system, bond_id) * bond_active[bond_id] + tzem = Cs .* params.bc * ωij / ω0 .* (Δxij .- F * ΔXij) # update of force density - tij = ωij * PKinv * ΔXij + Tij - if containsnan(tij) - tij = zero(SMatrix{3,3}) - end - update_add_b_int!(storage, i, tij .* system.volume[j]) - update_add_b_int!(storage, j, -tij .* system.volume[i]) + tij = ωij * PKinv * ΔXij + tzem + update_add_b_int!(storage, i, tij .* volume[j]) + update_add_b_int!(storage, j, -tij .* volume[i]) end return nothing end -@inline function influence_function(::NOSBMaterial, params::NOSBPointParameters, L::Float64) - return params.δ / L -end - -function calc_deformation_gradient(storage::NOSBStorage, system::BondSystem, - mat::NOSBMaterial, params::NOSBPointParameters, i::Int) - K = zeros(SMatrix{3,3}) - _F = zeros(SMatrix{3,3}) - ω0 = 0.0 - for bond_id in each_bond_idx(system, i) - bond = system.bonds[bond_id] - j, L = bond.neighbor, bond.length - ΔXij = get_coordinates_diff(system, i, j) - Δxij = get_coordinates_diff(storage, i, j) - Vj = system.volume[j] - ωij = influence_function(mat, params, L) * storage.bond_active[bond_id] - ω0 += ωij - temp = ωij * Vj - K += temp * ΔXij * ΔXij' - _F += temp * Δxij * ΔXij' - end - Kinv = inv(K) - F = _F * Kinv - return F, Kinv, ω0 -end - -function calc_first_piola_stress(F::SMatrix{3,3}, mat::NOSBMaterial, - params::NOSBPointParameters) - J = det(F) - J < eps() && return zero(SMatrix{3,3}) - J > mat.maxjacobi && return zero(SMatrix{3,3}) - C = F' * F - Cinv = inv(C) - S = params.G .* (I - 1 / 3 .* tr(C) .* Cinv) .* J^(-2 / 3) .+ - params.K / 4 .* (J^2 - J^(-2)) .* Cinv - P = F * S - return P -end - -function containsnan(K::T) where {T<:AbstractArray} - @simd for i in eachindex(K) - isnan(K[i]) && return true +function too_much_damage!(storage::AbstractStorage, system::AbstractSystem, + mat::AbstractCorrespondenceMaterial, defgrad_res, i) + (; F) = defgrad_res + # if storage.n_active_bonds[i] ≤ 3 || storage.damage[i] > mat.maxdmg || containsnan(F) + if 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 return false end - -function kill_point!(s::AbstractStorage, bd::BondSystem, i::Int) - s.bond_active[each_bond_idx(bd, i)] .= false - s.n_active_bonds[i] = 0 - return nothing -end diff --git a/src/physics/ordinary_state_based.jl b/src/physics/ordinary_state_based.jl index d92d86e8..b38ff773 100644 --- a/src/physics/ordinary_state_based.jl +++ b/src/physics/ordinary_state_based.jl @@ -54,9 +54,15 @@ When specifying the `fields` keyword of [`Job`](@ref) for a [`Body`](@ref) with - `damage::Vector{Float64}`: Damage of each point - `n_active_bonds::Vector{Int}`: Number of intact bonds of each point """ -struct OSBMaterial{Correction} <: AbstractBondSystemMaterial{Correction} end +struct OSBMaterial{Correction,K} <: AbstractBondSystemMaterial{Correction} + kernel::K + function OSBMaterial{Correction}(kernel::K) where {Correction,K} + return new{Correction,K}(kernel) + end +end -OSBMaterial() = OSBMaterial{NoCorrection}() +OSBMaterial{C}(; kernel::F=linear_kernel) where{C,F} = OSBMaterial{C}(kernel) +OSBMaterial(; kwargs...) = OSBMaterial{NoCorrection}(; kwargs...) struct OSBPointParameters <: AbstractPointParameters δ::Float64 @@ -102,7 +108,7 @@ function init_field(::OSBMaterial, ::AbstractTimeSolver, system::BondSystem, ::V end function force_density_point!(storage::OSBStorage, system::BondSystem, mat::OSBMaterial, - params::OSBPointParameters, i::Int) + params::OSBPointParameters, t, Δt, i) wvol = calc_weighted_volume(storage, system, mat, params, i) iszero(wvol) && return nothing dil = calc_dilatation(storage, system, mat, params, wvol, i) @@ -115,7 +121,7 @@ function force_density_point!(storage::OSBStorage, system::BondSystem, mat::OSBM l = norm(Δxij) ε = (l - L) / L stretch_based_failure!(storage, system, bond, params, ε, i, bond_id) - p_int = influence_function(mat, params, L) * bond_failure(storage, bond_id) * + p_int = kernel(system, bond_id) * bond_failure(storage, bond_id) * surface_correction_factor(system.correction, bond_id) * (c2 * L + c1 * (l - L)) / l .* Δxij update_add_b_int!(storage, i, p_int .* system.volume[j]) @@ -125,7 +131,7 @@ function force_density_point!(storage::OSBStorage, system::BondSystem, mat::OSBM end function force_density_point!(storage::OSBStorage, system::BondSystem, mat::OSBMaterial, - paramhandler::ParameterHandler, i::Int) + paramhandler::ParameterHandler, t, Δt, i) params_i = get_params(paramhandler, i) wvol = calc_weighted_volume(storage, system, mat, params_i, i) iszero(wvol) && return nothing @@ -140,7 +146,7 @@ function force_density_point!(storage::OSBStorage, system::BondSystem, mat::OSBM params_j = get_params(paramhandler, j) c1 = 15.0 * (params_i.G + params_j.G) / (2 * wvol) c2 = dil * (3.0 * (params_i.K + params_j.K) / (2 * wvol) - c1 / 3.0) - p_int = influence_function(mat, params_i, L) * bond_failure(storage, bond_id) * + p_int = kernel(system, bond_id) * bond_failure(storage, bond_id) * surface_correction_factor(system.correction, bond_id) * (c2 * L + c1 * (l - L)) / l .* Δxij update_add_b_int!(storage, i, p_int .* system.volume[j]) @@ -149,12 +155,8 @@ function force_density_point!(storage::OSBStorage, system::BondSystem, mat::OSBM return nothing end -@inline function influence_function(::OSBMaterial, params::OSBPointParameters, L::Float64) - return params.δ / L -end - function calc_weighted_volume(storage::OSBStorage, system::BondSystem, mat::OSBMaterial, - params::OSBPointParameters, i::Int) + params::OSBPointParameters, i) wvol = 0.0 for bond_id in each_bond_idx(system, i) bond = system.bonds[bond_id] @@ -162,14 +164,14 @@ function calc_weighted_volume(storage::OSBStorage, system::BondSystem, mat::OSBM ΔXij = get_coordinates_diff(system, i, j) ΔXij_sq = dot(ΔXij, ΔXij) scfactor = surface_correction_factor(system.correction, bond_id) - ωij = influence_function(mat, params, L) * storage.bond_active[bond_id] * scfactor + ωij = kernel(system, bond_id) * storage.bond_active[bond_id] * scfactor wvol += ωij * ΔXij_sq * system.volume[j] end return wvol end function calc_dilatation(storage::OSBStorage, system::BondSystem, mat::OSBMaterial, - params::OSBPointParameters, wvol::Float64, i::Int) + params::OSBPointParameters, wvol, i) dil = 0.0 c1 = 3.0 / wvol for bond_id in each_bond_idx(system, i) @@ -178,7 +180,7 @@ function calc_dilatation(storage::OSBStorage, system::BondSystem, mat::OSBMateri Δxij = get_coordinates_diff(storage, i, j) l = norm(Δxij) scfactor = surface_correction_factor(system.correction, bond_id) - ωij = influence_function(mat, params, L) * storage.bond_active[bond_id] * scfactor + ωij = kernel(system, bond_id) * storage.bond_active[bond_id] * scfactor dil += ωij * c1 * L * (l - L) * system.volume[j] end return dil diff --git a/src/physics/stress.jl b/src/physics/stress.jl new file mode 100644 index 00000000..504e7dd7 --- /dev/null +++ b/src/physics/stress.jl @@ -0,0 +1,102 @@ +function von_mises_stress(σ) + σx, σy, σz = σ[1], σ[5], σ[9] + τxy, τxz, τyz = σ[4], σ[7], σ[8] + a = σx * σx + σy * σy + σz * σz + b = - σx * σy - σx * σz - σy * σz + c = 3 * (τxy * τxy + τxz * τxz + τyz * τyz) + σvm = √(a + b + c) + return σvm +end + +function cauchy_stress(P, F) + J = det(F) + σ = 1/J .* P * F' + return σ +end + +function first_piola_kirchhoff(σ, F) + J = det(F) + P = J * σ * inv(F)' + return P +end + +function init_stress_rotation!(storage::AbstractStorage, F, Ḟ, Δt, i) + # inverse of the deformation gradient + F⁻¹ = inv(F) + if containsnan(F⁻¹) + OTensor = zero(SMatrix{3,3,Float64,9}) + update_tensor!(storage.rotation, i, OTensor) + update_tensor!(storage.left_stretch, i, OTensor) + return OTensor + end + + # Eulerian velocity gradient [FT87, eq. (3)] + L = Ḟ * F⁻¹ + + # rate-of-deformation tensor D + D = 0.5 .* (L + L') + + # spin tensor W + W = 0.5 .* (L - L') + + # left stretch V + V = get_tensor(storage.left_stretch, i) + + # vector z [FT87, eq. (13)] + z_x = - V[1,3] * D[2,1] - V[2,3] * D[2,2] - + V[3,3] * D[2,3] + V[1,2] * D[3,1] + + V[2,2] * D[3,2] + V[3,2] * D[3,3] + z_y = V[1,3] * D[1,1] + V[2,3] * D[1,2] + + V[3,3] * D[1,3] - V[1,1] * D[3,1] - + V[2,1] * D[3,2] - V[3,1] * D[3,3] + z_z = - V[1,2] * D[1,1] - V[2,2] * D[1,2] - + V[3,2] * D[1,3] + V[1,1] * D[2,1] + + V[2,1] * D[2,2] + V[3,1] * D[2,3] + z = SVector{3}(z_x, z_y, z_z) + + # w = -1/2 * \epsilon_{ijk} * W_{jk} [FT87, eq. (11)] + w = 0.5 .* SVector{3}(W[3,2] - W[2,3], W[1,3] - W[3,1], W[2,1] - W[1,2]) + + # ω = w + (I * tr(V) - V)^(-1) * z [FT87, eq. (12)] + ω = w + inv(I * tr(V) - V) * z + + # Ω [FT87, eq. (10)] + Ωtens = SMatrix{3,3}(0.0, -ω[3], ω[2], ω[3], 0.0, -ω[1], -ω[2], ω[1], 0.0) + Ω² = dot(ω, ω) + Ω = sqrt(Ω²) + + # compute Q with [FT87, eq. (44)] + if Ω² > 1e-30 # avoid a potential divide-by-zero + fac1 = sin(Δt * Ω) / Ω + fac2 = -(1.0 - cos(Δt * Ω)) / Ω² + Ωtens² = Ωtens * Ωtens + Q = I + fac1 .* Ωtens + fac2 .* Ωtens² + else + Q = SMatrix{3,3}(1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0) + end + + # compute Rotation of new step [FT87, eq. (36)] + R = get_tensor(storage.rotation, i) + Rₙ₊₁ = Q * R + + # compute step 4 of [FT87] + V̇ = L * V - V * Ωtens + + # compute step 5 of [FT87] + Vₙ₊₁ = V + Δt * V̇ + + # compute the unrotated rate of deformation + Dᵣ = Rₙ₊₁' * D * Rₙ₊₁ + + # update rotation and left stretch + update_tensor!(storage.rotation, i, Rₙ₊₁) + update_tensor!(storage.left_stretch, i, Vₙ₊₁) + # update_tensor!(storage.unrot_rate_of_deformation, i, Dᵣ) + return Dᵣ +end + +function rotate_stress(storage::AbstractStorage, σ, i) + R = get_tensor(storage.rotation, i) + T = R * σ * R' + return T +end diff --git a/src/time_solvers/dynamic_relaxation.jl b/src/time_solvers/dynamic_relaxation.jl index 9e0e648f..cdb661d4 100644 --- a/src/time_solvers/dynamic_relaxation.jl +++ b/src/time_solvers/dynamic_relaxation.jl @@ -165,20 +165,15 @@ end end function relaxation_timestep!(dh::AbstractThreadsBodyDataHandler, - options::AbstractJobOptions, - Δt::Float64, n::Int) + options::AbstractJobOptions, Δt, n) t = n * Δt @threads :static for chunk_id in eachindex(dh.chunks) chunk = dh.chunks[chunk_id] apply_boundary_conditions!(chunk, t) update_disp_and_pos!(chunk, Δt) end + calc_force_density!(dh, 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]) - end - @threads :static for chunk_id in eachindex(dh.chunks) - exchange_halo_to_loc!(dh, chunk_id) chunk = dh.chunks[chunk_id] calc_damage!(chunk) cn = calc_damping(chunk, Δt) @@ -192,6 +187,24 @@ function relaxation_timestep!(dh::AbstractThreadsBodyDataHandler, return nothing end +function relaxation_timestep!(dh::AbstractMPIBodyDataHandler, + options::AbstractJobOptions, Δt, n) + t = n * Δt + (; chunk) = dh + @timeit_debug TO "apply_boundary_conditions!" apply_boundary_conditions!(chunk, t) + @timeit_debug TO "update_disp_and_pos!" update_disp_and_pos!(chunk, Δt) + calc_force_density!(dh, t, Δt) + @timeit_debug TO "calc_damage!" calc_damage!(chunk) + @timeit_debug TO "relaxation_step!" if n == 1 + relaxation_first_step!(chunk, Δt) + else + cn = calc_damping(chunk, Δt) + relaxation_step!(chunk, Δt, cn) + end + @timeit_debug TO "export_results" export_results(dh, options, n, t) + return nothing +end + function calc_damping(chunk::AbstractBodyChunk, Δt::Float64) s = chunk.storage cn1 = 0.0 diff --git a/src/time_solvers/velocity_verlet.jl b/src/time_solvers/velocity_verlet.jl index f62b0b99..edaf4253 100644 --- a/src/time_solvers/velocity_verlet.jl +++ b/src/time_solvers/velocity_verlet.jl @@ -187,12 +187,8 @@ function verlet_timestep!(dh::AbstractThreadsBodyDataHandler, options::AbstractJ apply_boundary_conditions!(chunk, t) update_disp_and_pos!(chunk, Δt) end + calc_force_density!(dh, 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]) - end - @threads :static for chunk_id in eachindex(dh.chunks) - exchange_halo_to_loc!(dh, chunk_id) chunk = dh.chunks[chunk_id] calc_damage!(chunk) update_acc_and_vel!(chunk, Δt½) @@ -202,7 +198,7 @@ function verlet_timestep!(dh::AbstractThreadsBodyDataHandler, options::AbstractJ end function verlet_timestep!(dh::AbstractThreadsMultibodyDataHandler, - options::AbstractJobOptions, Δt::Float64, Δt½::Float64, n::Int) + options::AbstractJobOptions, Δt, Δt½, n) t = n * Δt for body_idx in each_body_idx(dh) body_dh = get_body_dh(dh, body_idx) @@ -212,16 +208,12 @@ function verlet_timestep!(dh::AbstractThreadsMultibodyDataHandler, apply_boundary_conditions!(chunk, t) update_disp_and_pos!(chunk, Δt) end - @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]) - end end + calc_force_density!(dh, t, Δt) update_caches!(dh) calc_contact_force_densities!(dh) for body_idx in each_body_idx(dh) body_dh = get_body_dh(dh, body_idx) - body_name = get_body_name(dh, body_idx) @threads :static for chunk_id in eachindex(body_dh.chunks) exchange_halo_to_loc!(body_dh, chunk_id) chunk = body_dh.chunks[chunk_id] @@ -240,9 +232,7 @@ function verlet_timestep!(dh::AbstractMPIBodyDataHandler, options::AbstractJobOp @timeit_debug TO "update_vel_half!" update_vel_half!(chunk, Δt½) @timeit_debug TO "apply_boundary_conditions!" apply_boundary_conditions!(chunk, t) @timeit_debug TO "update_disp_and_pos!" update_disp_and_pos!(chunk, Δt) - @timeit_debug TO "exchange_loc_to_halo!" exchange_loc_to_halo!(dh) - @timeit_debug TO "calc_force_density!" calc_force_density!(chunk) - @timeit_debug TO "exchange_halo_to_loc!" exchange_halo_to_loc!(dh) + calc_force_density!(dh, t, Δt) @timeit_debug TO "calc_damage!" calc_damage!(chunk) @timeit_debug TO "update_acc_and_vel!" update_acc_and_vel!(chunk, Δt½) @timeit_debug TO "export_results" export_results(dh, options, n, t) diff --git a/test/VtkReader/test_read_vtk.jl b/test/VtkReader/test_read_vtk.jl index 9ba7bee6..21a3ec2e 100644 --- a/test/VtkReader/test_read_vtk.jl +++ b/test/VtkReader/test_read_vtk.jl @@ -3,9 +3,9 @@ @testitem "read complete results" begin using Peridynamics.WriteVTK - bname = joinpath(@__DIR__, "vtk_test_1") + root = mktempdir() + bname = joinpath(root, "vtk_test_1") name = bname * ".vtu" - isfile(name) && rm(name) n_points = 10 position = rand(3, n_points) cells = [MeshCell(VTKCellTypes.VTK_VERTEX, (j,)) for j in 1:n_points] @@ -26,17 +26,15 @@ @test result[:Time] == [time] @test result[:Damage] == damage @test result[:Displacement] == displacement - - rm(name) end #-- read incomplete results @testitem "read incomplete results" begin using Peridynamics.WriteVTK - bname = joinpath(@__DIR__, "vtk_test_2") + root = mktempdir() + bname = joinpath(root, "vtk_test_2") name = bname * ".vtu" - isfile(name) && rm(name) n_points = 10 position = rand(3, n_points) cells = [MeshCell(VTKCellTypes.VTK_VERTEX, (j,)) for j in 1:n_points] @@ -58,8 +56,6 @@ end @test result[:damage] == damage @test result[:random_point_data] == random_point_data @test result[:random_field_data] == random_field_data - - rm(name) end @testitem "wrong file type" begin @@ -69,9 +65,9 @@ end @testitem "corrupt file raw encoding" begin using Peridynamics.WriteVTK - bname = joinpath(@__DIR__, "vtk_test_3") + root = mktempdir() + bname = joinpath(root, "vtk_test_3") name = bname * ".vtu" - isfile(name) && rm(name) n_points = 10 position = rand(3, n_points) cells = [MeshCell(VTKCellTypes.VTK_VERTEX, (j,)) for j in 1:n_points] @@ -92,16 +88,14 @@ end end @test_throws ErrorException read_vtk(name) - - rm(name) end @testitem "corrupt file offset marker" begin using Peridynamics.WriteVTK - bname = joinpath(@__DIR__, "vtk_test_4") + root = mktempdir() + bname = joinpath(root, "vtk_test_4") name = bname * ".vtu" - isfile(name) && rm(name) n_points = 10 position = rand(3, n_points) cells = [MeshCell(VTKCellTypes.VTK_VERTEX, (j,)) for j in 1:n_points] @@ -122,17 +116,15 @@ end end @test_throws ErrorException read_vtk(name) - - rm(name) end #-- corrupt file @testitem "corrupt file appended data" begin using Peridynamics.WriteVTK - bname = joinpath(@__DIR__, "vtk_test_5") + root = mktempdir() + bname = joinpath(root, "vtk_test_5") name = bname * ".vtu" - isfile(name) && rm(name) n_points = 10 position = rand(3, n_points) cells = [MeshCell(VTKCellTypes.VTK_VERTEX, (j,)) for j in 1:n_points] @@ -153,18 +145,14 @@ end end @test_throws ErrorException read_vtk(name) - - rm(name) end @testitem "read pvtu file" begin using Peridynamics.WriteVTK - root = joinpath(@__DIR__, "tempdir_1") - isdir(root) && rm(root; recursive=true, force=true) + root = mktempdir() bname = joinpath(root, "vtk_test") name = bname * ".pvtu" - isfile(name) && rm(name) n_points = 10 n_parts = 3 position = [rand(3, n_points) for _ in 1:n_parts] @@ -189,18 +177,14 @@ end @test result[:displacement] ≈ reduce(hcat, displacement) @test result[:integer_test] ≈ reduce(vcat, integer_test) @test result[:time] ≈ [time] - - rm(root; recursive=true, force=true) end @testitem "read pvtu files with nothing exported" begin using Peridynamics.WriteVTK - root = joinpath(@__DIR__, "tempdir_2") - isdir(root) && rm(root; recursive=true, force=true) + root = mktempdir() bname = joinpath(root, "vtk_test") name = bname * ".pvtu" - isfile(name) && rm(name) n_points = 10 n_parts = 3 position = [rand(3, n_points) for _ in 1:n_parts] @@ -213,16 +197,14 @@ end result = read_vtk(name) @test result[:position] ≈ reduce(hcat, position) - - rm(root; recursive=true, force=true) end @testitem "read vtu files with nothing exported" begin using Peridynamics.WriteVTK - bname = joinpath(@__DIR__, "vtk_test_6") + root = mktempdir() + bname = joinpath(root, "vtk_test_6") name = bname * ".vtu" - isfile(name) && rm(name) n_points = 10 position = rand(3, n_points) cells = [MeshCell(VTKCellTypes.VTK_VERTEX, (j,)) for j in 1:n_points] @@ -232,6 +214,4 @@ end @test typeof(result) == Dict{Symbol,VecOrMat{Float64}} @test result[:position] == position - - rm(name) end diff --git a/test/auxiliary/test_io.jl b/test/auxiliary/test_io.jl index 1a8fe196..dc5d3303 100644 --- a/test/auxiliary/test_io.jl +++ b/test/auxiliary/test_io.jl @@ -209,17 +209,17 @@ # setup bbb = Body(BBMaterial(), rand(3, 10), rand(10)) bosb = Body(OSBMaterial(), rand(3, 10), rand(10)) - bnosb = Body(NOSBMaterial(), rand(3, 10), rand(10)) - ms = MultibodySetup(:a => bbb, :b => bosb, :c => bnosb) + bcc = Body(CMaterial(), rand(3, 10), rand(10)) + ms = MultibodySetup(:a => bbb, :b => bosb, :c => bcc) vv = VelocityVerlet(steps=1) dr = DynamicRelaxation(steps=1) tests_body_specific(bbb, vv) tests_body_specific(bosb, vv) - tests_body_specific(bnosb, vv) + tests_body_specific(bcc, vv) tests_body_specific(bbb, dr) tests_body_specific(bosb, dr) - tests_body_specific(bnosb, dr) + tests_body_specific(bcc, dr) tests_multibody_specific(ms, vv) @@ -244,8 +244,7 @@ end @testitem "export_results Body" begin - temp_root = joinpath(@__DIR__, "temp_root_export_results_test_body") - rm(temp_root; recursive=true, force=true) + temp_root = mktempdir() pos, vol = uniform_box(1, 1, 1, 0.5) body = Body(BBMaterial(), pos, vol) @@ -283,13 +282,10 @@ end @test r[:displacement] ≈ u_rand @test r[:damage] == zeros(8) @test r[:time] == [t] - - rm(temp_root; recursive=true, force=true) end @testitem "export_results MultibodySetup" begin - temp_root = joinpath(@__DIR__, "temp_root_export_results_test_multibody") - rm(temp_root; recursive=true, force=true) + temp_root = mktempdir() pos_a, vol_a = uniform_box(1, 1, 1, 0.5) body_a = Body(BBMaterial(), pos_a, vol_a) @@ -359,6 +355,4 @@ end @test r_b[:displacement] ≈ u_rand_b @test r_b[:damage] == zeros(8) @test r_b[:time] == [t] - - rm(temp_root; recursive=true, force=true) end diff --git a/test/auxiliary/test_process_each_export.jl b/test/auxiliary/test_process_each_export.jl index 9aceacaf..e97545d7 100644 --- a/test/auxiliary/test_process_each_export.jl +++ b/test/auxiliary/test_process_each_export.jl @@ -1,6 +1,5 @@ @testitem "process_each_export" begin - root = joinpath(@__DIR__, "temp_root_proc_each_export") - isdir(root) && rm(root; recursive=true, force=true) + root = mktempdir() root_post_threads = joinpath(root, "post_threads") mkpath(root_post_threads) root_post_serial = joinpath(root, "post_serial") @@ -83,17 +82,12 @@ file_3_mpi = joinpath(root_post_mpi, "max_displacement_3.txt") @test isfile(file_3_mpi) @test contains(read(file_3_mpi, String), "maximum displacement x: 2.4") - - rm(root; recursive=true, force=true) end @testitem "find_vtk_files" begin - root = joinpath(@__DIR__, "temp_root_find_vtk_files") - isdir(root) && rm(root; recursive=true, force=true) - + root = mktempdir() vtk_files_unsorted = ["timestep_123.pvtu", "abcd_timestep_000005.pvtu", "timestep_02.pvtu", "_1.pvtu"] - mkpath(root) for file in vtk_files_unsorted open(joinpath(root, file), "w+") do io write(io, "") @@ -103,6 +97,4 @@ end vtk_files = Peridynamics.find_vtk_files(root) @test basename.(vtk_files) == ["_1.pvtu", "timestep_02.pvtu", "abcd_timestep_000005.pvtu", "timestep_123.pvtu"] - - rm(root; recursive=true, force=true) end diff --git a/test/discretization/test_bond_associated_system.jl b/test/discretization/test_bond_associated_system.jl new file mode 100644 index 00000000..d734146b --- /dev/null +++ b/test/discretization/test_bond_associated_system.jl @@ -0,0 +1,72 @@ +@testitem "bond-associated family" begin + # setup + position = [0.0 1.0 0.0 0.0 + 0.0 0.0 1.0 0.0 + 0.0 0.0 0.0 1.0] + volume = [1.1, 1.2, 1.3, 1.4] + mat = BACMaterial() + body = Body(mat, position, volume) + material!(body, horizon=2, rho=1, E=1, nu=0.25, Gc=1) + pd = Peridynamics.PointDecomposition(body, 1) + + # 1 + system = Peridynamics.BondAssociatedSystem(body, pd, 1) + + @test system.position == position + @test system.volume == volume + @test system.bonds == [ + Peridynamics.Bond(2, 1.0, true), + Peridynamics.Bond(3, 1.0, true), + Peridynamics.Bond(4, 1.0, true), + Peridynamics.Bond(1, 1.0, true), + Peridynamics.Bond(3, √2, true), + Peridynamics.Bond(4, √2, true), + Peridynamics.Bond(1, 1.0, true), + Peridynamics.Bond(2, √2, true), + Peridynamics.Bond(4, √2, true), + Peridynamics.Bond(1, 1.0, true), + Peridynamics.Bond(2, √2, true), + Peridynamics.Bond(3, √2, true), + ] + @test system.n_neighbors == [3, 3, 3, 3] + @test system.bond_ids == [1:3, 4:6, 7:9, 10:12] + +end + +@testitem "intersection bond ids" begin + pos, vol = uniform_box(1, 0.25, 0.25, 0.25) + body = Body(BACMaterial(), pos, vol) + material!(body, horizon=0.5, bond_horizon=0.5, rho=1, E=1, nu=0.25, Gc=1) + + pd = Peridynamics.PointDecomposition(body, 1) + + bas = Peridynamics.get_system(body, pd, 1) + (; bonds, bond_ids, intersection_bond_ids) = bas + + @test bonds == [ + Peridynamics.Bond(2, 0.25, true) + Peridynamics.Bond(3, 0.5, true) + Peridynamics.Bond(1, 0.25, true) + Peridynamics.Bond(3, 0.25, true) + Peridynamics.Bond(4, 0.5, true) + Peridynamics.Bond(1, 0.5, true) + Peridynamics.Bond(2, 0.25, true) + Peridynamics.Bond(4, 0.25, true) + Peridynamics.Bond(2, 0.5, true) + Peridynamics.Bond(3, 0.25, true) + ] + @test bond_ids == [1:2, 3:5, 6:8, 9:10] + i_bond_ids = [[1, 2], [1, 2], [1], [2, 3], [2, 3], [1, 2], [1, 2], [3], [1, 2], [1, 2]] + @test intersection_bond_ids == i_bond_ids + + # for i in Peridynamics.each_point_idx(system) + # for bond_idx in Peridynamics.each_bond_idx(system, i) + # # bond = bonds[bond_idx] + # # j = bond.neighbor + # for babond_idx in Peridynamics.each_intersecting_bond_idx(system, i, bond_idx) + + # end + # end + # end + +end diff --git a/test/integration/b_int_bond_based.jl b/test/integration/b_int_bond_based.jl index c0745eea..08956302 100644 --- a/test/integration/b_int_bond_based.jl +++ b/test/integration/b_int_bond_based.jl @@ -19,7 +19,7 @@ # Point 2 with v_z = 1 m/s with Δt = 0.0015 s position[1, 2] = 1.0015 - Peridynamics.calc_force_density!(chunk) + Peridynamics.calc_force_density!(chunk, 0, 0) b12 = 18 * E / (3 * (1 - 2 * 0.25)) / (π * δ^4) * 1.0015 * 0.0015/1.0015 * 1.0 @test b_int ≈ [b12 -b12; 0.0 0.0; 0.0 0.0] @@ -51,7 +51,7 @@ end # Point 2 with v_z = 1 m/s with Δt = 0.0015 s position[1, 2] = 1.0015 - Peridynamics.calc_force_density!(chunk) + Peridynamics.calc_force_density!(chunk, 0, 0) b12 = 18 * E / (3 * (1 - 2 * 0.25)) / (π * δ^4) * 1.0015 * 0.0015/1.0015 * 1.0 @test b_int ≈ [b12 -b12; 0.0 0.0; 0.0 0.0] diff --git a/test/integration/b_int_continuum_inspired.jl b/test/integration/b_int_continuum_inspired.jl index 7760ff6a..2513bc5e 100644 --- a/test/integration/b_int_continuum_inspired.jl +++ b/test/integration/b_int_continuum_inspired.jl @@ -21,7 +21,7 @@ # Point 2 with v_z = 1 m/s with Δt = 0.0015 s position[1, 2] = 1.0015 - Peridynamics.calc_force_density!(chunk) + Peridynamics.calc_force_density!(chunk, 0, 0) @test b_int[:,1] ≈ [1.0000000000000625e9, 4.0060000000002503e8, 4.0060000000002503e8] @test b_int[:,2] ≈ [-1.449731150462047e9, 2.245287820579052e8, 2.245287820579052e8] @@ -54,7 +54,7 @@ end # Point 2 with v_z = 1 m/s with Δt = 0.0015 s position[1, 2] = 1.0015 - Peridynamics.calc_force_density!(chunk) + Peridynamics.calc_force_density!(chunk, 0, 0) @test b_int[:,1] ≈ [1.0000000000000625e9, 4.0060000000002503e8, 4.0060000000002503e8] @test b_int[:,2] ≈ [-1.449731150462047e9, 2.245287820579052e8, 2.245287820579052e8] diff --git a/test/integration/b_int_correspondence.jl b/test/integration/b_int_correspondence.jl index d64d26c1..cf811cf4 100644 --- a/test/integration/b_int_correspondence.jl +++ b/test/integration/b_int_correspondence.jl @@ -4,7 +4,7 @@ 0.0 0.0 0.0 1.0 2.0] volume = fill(1.0, 5) δ = 1.5 - body = Body(NOSBMaterial(), ref_position, volume) + body = Body(CMaterial(model=MooneyRivlin()), ref_position, volume) material!(body, horizon=δ, rho=1, E=1, nu=0.25, Gc=1.0) failure_permit!(body, false) @@ -21,7 +21,7 @@ # Point 2 with v_z = 1 m/s with Δt = 0.0015 s position[1, 2] = 1.0015 - Peridynamics.calc_force_density!(chunk) + Peridynamics.calc_force_density!(chunk, 0, 0) @test b_int[:,1] ≈ [0.007185432432164514, 0.0023974081843325646, 0.0023974081843347313] @test b_int[:,2] ≈ [-0.007185432432123352, -9.047513544240966e-15, -6.0453612662025915e-15] @@ -36,7 +36,7 @@ end 0.0 0.0 0.0 1.0 2.0] volume = fill(1.0, 5) δ = 1.5 - body = Body(NOSBMaterial(), ref_position, volume) + body = Body(CMaterial(model=MooneyRivlin()), ref_position, volume) point_set!(body, :a, [1]) point_set!(body, :b, [2,3,4,5]) material!(body, :a, horizon=δ, rho=1, E=1, nu=0.25, Gc=1.0) @@ -54,7 +54,7 @@ end # Point 2 with v_z = 1 m/s with Δt = 0.0015 s position[1, 2] = 1.0015 - Peridynamics.calc_force_density!(chunk) + Peridynamics.calc_force_density!(chunk, 0, 0) @test b_int[:,1] ≈ [0.007185432432164514, 0.0023974081843325646, 0.0023974081843347313] @test b_int[:,2] ≈ [-0.007185432432123352, -9.047513544240966e-15, -6.0453612662025915e-15] diff --git a/test/integration/b_int_ordinary_state_based.jl b/test/integration/b_int_ordinary_state_based.jl index fbc530ce..d8debcb8 100644 --- a/test/integration/b_int_ordinary_state_based.jl +++ b/test/integration/b_int_ordinary_state_based.jl @@ -20,7 +20,7 @@ # Point 2 with v_z = 1 m/s with Δt = 0.0015 s position[1, 2] = 1.0015 - Peridynamics.calc_force_density!(chunk) + Peridynamics.calc_force_density!(chunk, 0, 0) b12 = 3.780000000000143e9 @test b_int ≈ [b12 -b12; 0.0 0.0; 0.0 0.0] @@ -51,7 +51,7 @@ end # Point 2 with v_z = 1 m/s with Δt = 0.0015 s position[1, 2] = 1.0015 - Peridynamics.calc_force_density!(chunk) + Peridynamics.calc_force_density!(chunk, 0, 0) b1 = 4.252500000000161e9 b2 = 2.1262500000000806e9 diff --git a/test/integration/material_interface.jl b/test/integration/material_interface.jl index 2b795e1d..2e9fae73 100644 --- a/test/integration/material_interface.jl +++ b/test/integration/material_interface.jl @@ -539,7 +539,7 @@ end @testitem "TestMaterial full simulation" begin include("test_material.jl") - root = joinpath(@__DIR__, "temp_testmaterial_full_simulation") + root = mktempdir() path_vtk = joinpath(root, "vtk") N = 10 @@ -562,6 +562,4 @@ end @test isdir(path_vtk) vtk_files = Peridynamics.find_vtk_files(path_vtk) @test length(vtk_files) == 3 - - rm(root; recursive=true, force=true) end diff --git a/test/integration/mpi_threads_comparison.jl b/test/integration/mpi_threads_comparison.jl index ea5505fa..8fdf5ab7 100644 --- a/test/integration/mpi_threads_comparison.jl +++ b/test/integration/mpi_threads_comparison.jl @@ -1,5 +1,5 @@ @testitem "MPI-Threads comparison BBMaterial{NoCorrection}" begin - root = joinpath(@__DIR__, "temp_mpi_threads_comparison_bbnc") + root = mktempdir() path_threads = joinpath(root, "results_threads") path_threads_vtk = joinpath(path_threads, "vtk") path_mpi = joinpath(root, "results_mpi") @@ -61,12 +61,82 @@ @test res_threads[key] ≈ res_mpi[key] end end +end + +@testitem "MPI-Threads comparison BBMaterial{NoCorrection} DynamicRelaxation" begin + root = mktempdir() + path_threads = joinpath(root, "results_threads") + path_threads_vtk = joinpath(path_threads, "vtk") + path_mpi = joinpath(root, "results_mpi") + path_mpi_vtk = joinpath(path_mpi, "vtk") + + function sim_bb(N::Int, path::String) + l, Δx, δ, a = 100.0, 100.0/N, 3.015/N, 50.0 + pos, vol = uniform_box(l, l, 0.1l, Δx) + ids = sortperm(pos[2,:]) + body = Body(BBMaterial(), pos[:, ids], vol[ids]) + failure_permit!(body, false) + material!(body; horizon=3.015Δx, E=2.1e5, rho=8e-6, Gc=2.7) + point_set!(p -> p[1] ≤ -l/2+a && 0 ≤ p[2] ≤ 2δ, body, :set_a) + point_set!(p -> p[1] ≤ -l/2+a && -2δ ≤ p[2] < 0, body, :set_b) + precrack!(body, :set_a, :set_b) + point_set!(p -> p[2] > l/2-Δx, body, :set_top) + point_set!(p -> p[2] < -l/2+Δx, body, :set_bottom) + forcedensity_bc!(t -> -1e3, body, :set_bottom, :y) + forcedensity_bc!(t -> 1e3, body, :set_top, :y) + solver = DynamicRelaxation(steps=200, damping_factor=0.05) + job = Job(body, solver; path=path, freq=200) + submit(job) + return nothing + end + sim_bb(30, path_threads) - rm(root; recursive=true, force=true) + mpi_cmd = """ + using Peridynamics + function sim_bb(N::Int, path::String) + l, Δx, δ, a = 100.0, 100.0/N, 3.015/N, 50.0 + pos, vol = uniform_box(l, l, 0.1l, Δx) + ids = sortperm(pos[2,:]) + body = Body(BBMaterial(), pos[:, ids], vol[ids]) + failure_permit!(body, false) + material!(body; horizon=3.015Δx, E=2.1e5, rho=8e-6, Gc=2.7) + point_set!(p -> p[1] ≤ -l/2+a && 0 ≤ p[2] ≤ 2δ, body, :set_a) + point_set!(p -> p[1] ≤ -l/2+a && -2δ ≤ p[2] < 0, body, :set_b) + precrack!(body, :set_a, :set_b) + point_set!(p -> p[2] > l/2-Δx, body, :set_top) + point_set!(p -> p[2] < -l/2+Δx, body, :set_bottom) + forcedensity_bc!(t -> -1e3, body, :set_bottom, :y) + forcedensity_bc!(t -> 1e3, body, :set_top, :y) + solver = DynamicRelaxation(steps=200, damping_factor=0.05) + job = Job(body, solver; path=path, freq=200) + submit(job) + return nothing + end + sim_bb(30, "$path_mpi") + """ + run(`$(Peridynamics.MPI.mpiexec()) -n 2 $(Base.julia_cmd()) --project -e $(mpi_cmd)`) + + @test isdir(path_threads_vtk) + @test isdir(path_mpi_vtk) + vtk_files_threads = Peridynamics.find_vtk_files(path_threads_vtk) + vtk_files_mpi = Peridynamics.find_vtk_files(path_mpi_vtk) + @test length(vtk_files_mpi) == length(vtk_files_threads) == 2 + for i in eachindex(vtk_files_threads, vtk_files_mpi) + res_threads = read_vtk(vtk_files_threads[i]) + res_mpi = read_vtk(vtk_files_mpi[i]) + for key in keys(res_threads) + res_threads_qty = res_threads[key] + res_mpi_qty = res_mpi[key] + Δe = maximum(abs.(res_threads_qty .- res_mpi_qty)) + #TODO: why is this error so high? + # See also: https://github.com/kaipartmann/Peridynamics.jl/issues/187 + @test Δe < 0.03 + end + end end @testitem "MPI-Threads comparison BBMaterial{EnergySurfaceCorrection}" begin - root = joinpath(@__DIR__, "temp_mpi_threads_comparison_bbesc") + root = mktempdir() path_threads = joinpath(root, "results_threads") path_threads_vtk = joinpath(path_threads, "vtk") path_mpi = joinpath(root, "results_mpi") @@ -128,12 +198,10 @@ end @test res_threads[key] ≈ res_mpi[key] end end - - rm(root; recursive=true, force=true) end @testitem "MPI-Threads comparison OSBMaterial" begin - root = joinpath(@__DIR__, "temp_mpi_threads_comparison_osb") + root = mktempdir() path_threads = joinpath(root, "results_threads") path_threads_vtk = joinpath(path_threads, "vtk") path_mpi = joinpath(root, "results_mpi") @@ -195,22 +263,20 @@ end @test res_threads[key] ≈ res_mpi[key] end end - - rm(root; recursive=true, force=true) end -@testitem "MPI-Threads comparison NOSBMaterial" begin - root = joinpath(@__DIR__, "temp_mpi_threads_comparison_cc") +@testitem "MPI-Threads comparison CMaterial" begin + root = mktempdir() path_threads = joinpath(root, "results_threads") path_threads_vtk = joinpath(path_threads, "vtk") path_mpi = joinpath(root, "results_mpi") path_mpi_vtk = joinpath(path_mpi, "vtk") - function sim_nosb(N::Int, path::String) + function sim_cc(N::Int, path::String) l, Δx, δ, a = 1.0, 1/N, 3.015/N, 0.5 pos, vol = uniform_box(l, l, 0.1l, Δx) ids = sortperm(pos[2,:]) - b = Body(NOSBMaterial(), pos[:, ids], vol[ids]) + b = Body(CMaterial(), pos[:, ids], vol[ids]) material!(b; horizon=3.015Δx, E=2.1e5, nu=0.25, rho=8e-6, Gc=2.7) point_set!(p -> p[1] ≤ -l/2+a && 0 ≤ p[2] ≤ 2δ, b, :set_a) point_set!(p -> p[1] ≤ -l/2+a && -2δ ≤ p[2] < 0, b, :set_b) @@ -224,15 +290,15 @@ end submit(job) return nothing end - sim_nosb(30, path_threads) + sim_cc(30, path_threads) mpi_cmd = """ using Peridynamics - function sim_nosb(N::Int, path::String) + function sim_cc(N::Int, path::String) l, Δx, δ, a = 1.0, 1/N, 3.015/N, 0.5 pos, vol = uniform_box(l, l, 0.1l, Δx) ids = sortperm(pos[2,:]) - b = Body(NOSBMaterial(), pos[:, ids], vol[ids]) + b = Body(CMaterial(), pos[:, ids], vol[ids]) material!(b; horizon=3.015Δx, E=2.1e5, nu=0.25, rho=8e-6, Gc=2.7) point_set!(p -> p[1] ≤ -l/2+a && 0 ≤ p[2] ≤ 2δ, b, :set_a) point_set!(p -> p[1] ≤ -l/2+a && -2δ ≤ p[2] < 0, b, :set_b) @@ -246,7 +312,7 @@ end submit(job) return nothing end - sim_nosb(30, "$path_mpi") + sim_cc(30, "$path_mpi") """ run(`$(Peridynamics.MPI.mpiexec()) -n 2 $(Base.julia_cmd()) --project -e $(mpi_cmd)`) @@ -262,6 +328,4 @@ end @test res_threads[key] ≈ res_mpi[key] end end - - rm(root; recursive=true, force=true) end diff --git a/test/integration/symmetry_bond_based.jl b/test/integration/symmetry_bond_based.jl index d3325f23..86dfe151 100644 --- a/test/integration/symmetry_bond_based.jl +++ b/test/integration/symmetry_bond_based.jl @@ -16,8 +16,7 @@ velocity_bc!(t -> 100, body, :set_a, :z) velocity_bc!(t -> -100, body, :set_b, :z) vv = VelocityVerlet(steps=100) - temp_root = joinpath(@__DIR__, "temp_root_symmetry_test_bb_vv") - rm(temp_root; recursive=true, force=true) + temp_root = mktempdir() job = Job(body, vv; path=temp_root, freq=10) dh = Peridynamics.submit_threads(job, 1) @@ -62,9 +61,6 @@ @test abs(end_displacement[d, i]) > 0 end end - - # delete all simulation files - rm(temp_root; recursive=true, force=true) end @testitem "symmetry BBMaterial DynamicRelaxation" begin @@ -85,8 +81,7 @@ end forcedensity_bc!(t -> 1e10, body, :set_a, :z) forcedensity_bc!(t -> -1e10, body, :set_b, :z) dr = DynamicRelaxation(steps=100) - temp_root = joinpath(@__DIR__, "temp_root_symmetry_test_bb_dr") - rm(temp_root; recursive=true, force=true) + temp_root = mktempdir() job = Job(body, dr; path=temp_root, freq=10) dh = Peridynamics.submit_threads(job, 1) @@ -121,7 +116,4 @@ end @test abs(end_displacement[d, i]) > 0 end end - - # delete all simulation files - rm(temp_root; recursive=true, force=true) end diff --git a/test/integration/symmetry_continuum_inspired.jl b/test/integration/symmetry_continuum_inspired.jl index 5b12131e..8cf8ccde 100644 --- a/test/integration/symmetry_continuum_inspired.jl +++ b/test/integration/symmetry_continuum_inspired.jl @@ -17,8 +17,7 @@ velocity_bc!(t -> 100, body, :set_a, :z) velocity_bc!(t -> -100, body, :set_b, :z) vv = VelocityVerlet(steps=20) - temp_root = joinpath(@__DIR__, "temp_root_symmetry_test_cki_vv") - rm(temp_root; recursive=true, force=true) + temp_root = mktempdir() job = Job(body, vv; path=temp_root, freq=10) dh = Peridynamics.submit_threads(job, 1) @@ -59,9 +58,6 @@ @test abs(end_displacement[d, i]) > 0 end end - - # delete all simulation files - rm(temp_root; recursive=true, force=true) end @testitem "symmetry CKIMaterial DynamicRelaxation" begin @@ -83,8 +79,7 @@ end forcedensity_bc!(t -> 1e10, body, :set_a, :z) forcedensity_bc!(t -> -1e10, body, :set_b, :z) dr = DynamicRelaxation(steps=20) - temp_root = joinpath(@__DIR__, "temp_root_symmetry_test_cki_dr") - rm(temp_root; recursive=true, force=true) + temp_root = mktempdir() job = Job(body, dr; path=temp_root, freq=10) dh = Peridynamics.submit_threads(job, 1) @@ -125,7 +120,4 @@ end @test abs(end_displacement[d, i]) > 0 end end - - # delete all simulation files - rm(temp_root; recursive=true, force=true) end diff --git a/test/integration/symmetry_correspondence.jl b/test/integration/symmetry_correspondence.jl index d6c7237f..4c0841f5 100644 --- a/test/integration/symmetry_correspondence.jl +++ b/test/integration/symmetry_correspondence.jl @@ -1,4 +1,4 @@ -@testitem "symmetry NOSBMaterial VelocityVerlet" begin +@testitem "symmetry CMaterial VelocityVerlet" begin # simulation Δx = 0.2 width = 1 @@ -8,16 +8,15 @@ pos = hcat(([x;y;z] for x in grid for y in grid for z in grid)...) n_points = size(pos, 2) vol = fill(Δx^3, n_points) - body = Body(NOSBMaterial(), pos, vol) + body = Body(CMaterial(model=MooneyRivlin()), pos, vol) failure_permit!(body, false) material!(body, horizon=3.015Δx, rho=7850, E=210e9, nu=0.25, Gc=1) point_set!(z -> z > width/2 - 0.6Δx, body, :set_a) point_set!(z -> z < -width/2 + 0.6Δx, body, :set_b) - velocity_bc!(t -> 100, body, :set_a, :z) - velocity_bc!(t -> -100, body, :set_b, :z) + velocity_bc!(t -> 10, body, :set_a, :z) + velocity_bc!(t -> -10, body, :set_b, :z) vv = VelocityVerlet(steps=100) - temp_root = joinpath(@__DIR__, "temp_root_symmetry_test_cc_vv") - rm(temp_root; recursive=true, force=true) + temp_root = mktempdir() job = Job(body, vv; path=temp_root, freq=10) dh = Peridynamics.submit_threads(job, 1) @@ -52,12 +51,9 @@ @test abs(end_displacement[d, i]) > 0 end end - - # delete all simulation files - rm(temp_root; recursive=true, force=true) end -@testitem "symmetry NOSBMaterial DynamicRelaxation" begin +@testitem "symmetry CMaterial DynamicRelaxation" begin # simulation Δx = 0.2 width = 1 @@ -67,7 +63,7 @@ end pos = hcat(([x;y;z] for x in grid for y in grid for z in grid)...) n_points = size(pos, 2) vol = fill(Δx^3, n_points) - body = Body(NOSBMaterial(), pos, vol) + body = Body(CMaterial(), pos, vol) failure_permit!(body, false) material!(body, horizon=3.015Δx, rho=7850, E=210e9, nu=0.25, Gc=1) point_set!(z -> z > width/2 - 0.6Δx, body, :set_a) @@ -75,8 +71,7 @@ end forcedensity_bc!(t -> 1e10, body, :set_a, :z) forcedensity_bc!(t -> -1e10, body, :set_b, :z) dr = DynamicRelaxation(steps=100) - temp_root = joinpath(@__DIR__, "temp_root_symmetry_test_cc_dr") - rm(temp_root; recursive=true, force=true) + temp_root = mktempdir() job = Job(body, dr; path=temp_root, freq=10) dh = Peridynamics.submit_threads(job, 1) @@ -111,7 +106,4 @@ end @test abs(end_displacement[d, i]) > 0 end end - - # delete all simulation files - rm(temp_root; recursive=true, force=true) end diff --git a/test/integration/symmetry_ordinary_state_based.jl b/test/integration/symmetry_ordinary_state_based.jl index e4edfbb4..92d6d5fc 100644 --- a/test/integration/symmetry_ordinary_state_based.jl +++ b/test/integration/symmetry_ordinary_state_based.jl @@ -16,8 +16,7 @@ velocity_bc!(t -> 100, body, :set_a, :z) velocity_bc!(t -> -100, body, :set_b, :z) vv = VelocityVerlet(steps=100) - temp_root = joinpath(@__DIR__, "temp_root_symmetry_test_osb_vv") - rm(temp_root; recursive=true, force=true) + temp_root = mktempdir() job = Job(body, vv; path=temp_root, freq=10) dh = Peridynamics.submit_threads(job, 1) @@ -52,9 +51,6 @@ @test abs(end_displacement[d, i]) > 0 end end - - # delete all simulation files - rm(temp_root; recursive=true, force=true) end @testitem "symmetry OSBMaterial DynamicRelaxation" begin @@ -75,8 +71,7 @@ end forcedensity_bc!(t -> 1e10, body, :set_a, :z) forcedensity_bc!(t -> -1e10, body, :set_b, :z) dr = DynamicRelaxation(steps=100) - temp_root = joinpath(@__DIR__, "temp_root_symmetry_test_osb_dr") - rm(temp_root; recursive=true, force=true) + temp_root = mktempdir() job = Job(body, dr; path=temp_root, freq=10) dh = Peridynamics.submit_threads(job, 1) @@ -111,7 +106,4 @@ end @test abs(end_displacement[d, i]) > 0 end end - - # delete all simulation files - rm(temp_root; recursive=true, force=true) end diff --git a/test/integration/test_material.jl b/test/integration/test_material.jl index 4d9c68de..d3954b25 100644 --- a/test/integration/test_material.jl +++ b/test/integration/test_material.jl @@ -40,7 +40,7 @@ function Peridynamics.init_field(::TestMaterial, ::Peridynamics.VelocityVerlet, end function Peridynamics.force_density_point!(storage::TestStorage, system::Peridynamics.BondSystem, ::TestMaterial, - params::TestPointParameters, i::Int) + params::TestPointParameters, t, Δt, i) for bond_id in Peridynamics.each_bond_idx(system, i) bond = system.bonds[bond_id] j, L = bond.neighbor, bond.length