Add halo exchange for write fields

src/core/threads_data_handler.jl

@@ -41,7 +41,7 @@ function exchange_read_fields!(dh::ThreadsDataHandler, chunk_id::Int)
 end
 
 function exchange_write_fields!(dh::ThreadsDataHandler, chunk_id::Int)
-    halo_exchange!(dh, dh.write_halo_exs[chunk_id])
+    halo_exchange_add!(dh, dh.write_halo_exs[chunk_id])
     return nothing
 end
 
@@ -54,6 +54,15 @@ function halo_exchange!(dh::ThreadsDataHandler, halo_exs::Vector{HaloExchange})
     return nothing
 end
 
+function halo_exchange_add!(dh::ThreadsDataHandler, halo_exs::Vector{HaloExchange})
+    for he in halo_exs
+        src_field = get_exchange_field(dh.chunks[he.src_chunk_id], he.field)
+        dest_field = get_exchange_field(dh.chunks[he.dest_chunk_id], he.field)
+        exchange_add!(dest_field, src_field, he.dest_idxs, he.src_idxs)
+    end
+    return nothing
+end
+
 @inline function get_exchange_field(b::AbstractBodyChunk, fieldname::Symbol)
     return getfield(b.store, fieldname)::Matrix{Float64}
 end
@@ -67,3 +76,13 @@ function exchange!(dest::Matrix{T}, src::Matrix{T}, dest_idxs::Vector{Int},
     end
     return nothing
 end
+
+function exchange_add!(dest::Matrix{T}, src::Matrix{T}, dest_idxs::Vector{Int},
+                   src_idxs::Vector{Int}) where {T<:Number}
+    for i in eachindex(dest_idxs)
+        for d in axes(dest, 1)
+            @inbounds dest[d, dest_idxs[i]] += src[d, src_idxs[i]]
+        end
+    end
+    return nothing
+end

test/integration/material_interface.jl

@@ -0,0 +1,129 @@
+
+@testitem "TestMaterial halo exchange" begin
+    include("test_material.jl")
+
+    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 = TestMaterial()
+    body = Body(mat, position, volume)
+    material!(body, horizon=2, rho=1, E=1, nu=0.25, Gc=1)
+    point_set!(body, :a, 1:2)
+    point_set!(body, :b, 3:4)
+    velocity_ic!(body, :a, :x, 1.0)
+    velocity_bc!(t -> t, body, :a, :x)
+    forcedensity_bc!(t -> t, body, :a, :x)
+    precrack!(body, :a, :b)
+    ts = VelocityVerlet(steps=10)
+    point_decomp = Peridynamics.PointDecomposition(body, 2)
+    tdh = Peridynamics.ThreadsDataHandler(body, ts, point_decomp)
+
+    b1 = tdh.chunks[1]
+    @test b1 isa Peridynamics.BodyChunk
+    @test b1.store.position == position
+    @test b1.store.displacement == zeros(3, 2)
+    @test b1.store.velocity == [1.0 1.0; 0.0 0.0; 0.0 0.0]
+    @test b1.store.velocity_half == zeros(3, 2)
+    @test b1.store.acceleration == zeros(3, 2)
+    @test b1.store.b_int == zeros(3, 4)
+    @test b1.store.b_ext == zeros(3, 2)
+    @test b1.store.damage ≈ [2/3, 2/3]
+    @test b1.store.bond_active == [1, 0, 0, 1, 0, 0]
+    @test b1.store.n_active_bonds == [1, 1]
+    b2 = tdh.chunks[2]
+    @test b2 isa Peridynamics.BodyChunk
+    @test b2.store.position == position[:, [3, 4, 1, 2]]
+    @test b2.store.displacement == zeros(3, 2)
+    @test b2.store.velocity == [0.0 0.0; 0.0 0.0; 0.0 0.0]
+    @test b2.store.velocity_half == zeros(3, 2)
+    @test b2.store.acceleration == zeros(3, 2)
+    @test b2.store.b_int == zeros(3, 4)
+    @test b2.store.b_ext == zeros(3, 2)
+    @test b2.store.damage ≈ [2/3, 2/3]
+    @test b2.store.bond_active == [0, 0, 1, 0, 0, 1]
+    @test b2.store.n_active_bonds == [1, 1]
+
+    randpos = rand(3, 4)
+    tdh.chunks[2].store.position .= randpos
+    Peridynamics.exchange_read_fields!(tdh, 1)
+
+    @test b1.store.position[:,1:2] ≈ [0.0 1.0; 0.0 0.0; 0.0 0.0]
+    @test b1.store.position[:,3:4] ≈ randpos[:,1:2]
+    @test b1.store.displacement == zeros(3, 2)
+    @test b1.store.velocity == [1.0 1.0; 0.0 0.0; 0.0 0.0]
+    @test b1.store.velocity_half == zeros(3, 2)
+    @test b1.store.acceleration == zeros(3, 2)
+    @test b1.store.b_int == zeros(3, 4)
+    @test b1.store.b_ext == zeros(3, 2)
+    @test b1.store.damage ≈ [2/3, 2/3]
+    @test b1.store.bond_active == [1, 0, 0, 1, 0, 0]
+    @test b1.store.n_active_bonds == [1, 1]
+
+    @test b2.store.position ≈ randpos
+    @test b2.store.displacement == zeros(3, 2)
+    @test b2.store.velocity == [0.0 0.0; 0.0 0.0; 0.0 0.0]
+    @test b2.store.velocity_half == zeros(3, 2)
+    @test b2.store.acceleration == zeros(3, 2)
+    @test b2.store.b_int == zeros(3, 4)
+    @test b2.store.b_ext == zeros(3, 2)
+    @test b2.store.damage ≈ [2/3, 2/3]
+    @test b2.store.bond_active == [0, 0, 1, 0, 0, 1]
+    @test b2.store.n_active_bonds == [1, 1]
+
+    randbint = rand(3, 4)
+    b2.store.b_int .= randbint
+    b1.store.b_int .+= 1
+    Peridynamics.exchange_write_fields!(tdh, 1)
+
+    @test b1.store.position[:,1:2] ≈ [0.0 1.0; 0.0 0.0; 0.0 0.0]
+    @test b1.store.position[:,3:4] ≈ randpos[:,1:2]
+    @test b1.store.displacement == zeros(3, 2)
+    @test b1.store.velocity == [1.0 1.0; 0.0 0.0; 0.0 0.0]
+    @test b1.store.velocity_half == zeros(3, 2)
+    @test b1.store.acceleration == zeros(3, 2)
+    @test b1.store.b_int[:,1:2] ≈ 1 .+ randbint[:,3:4]
+    @test b1.store.b_int[:,3:4] ≈ ones(3, 2)
+    @test b1.store.b_ext == zeros(3, 2)
+    @test b1.store.damage ≈ [2/3, 2/3]
+    @test b1.store.bond_active == [1, 0, 0, 1, 0, 0]
+    @test b1.store.n_active_bonds == [1, 1]
+
+    @test b2.store.position ≈ randpos
+    @test b2.store.displacement == zeros(3, 2)
+    @test b2.store.velocity == [0.0 0.0; 0.0 0.0; 0.0 0.0]
+    @test b2.store.velocity_half == zeros(3, 2)
+    @test b2.store.acceleration == zeros(3, 2)
+    @test b2.store.b_int ≈ randbint
+    @test b2.store.b_ext == zeros(3, 2)
+    @test b2.store.damage ≈ [2/3, 2/3]
+    @test b2.store.bond_active == [0, 0, 1, 0, 0, 1]
+    @test b2.store.n_active_bonds == [1, 1]
+
+    Peridynamics.exchange_write_fields!(tdh, 2)
+
+    @test b1.store.position[:,1:2] ≈ [0.0 1.0; 0.0 0.0; 0.0 0.0]
+    @test b1.store.position[:,3:4] ≈ randpos[:,1:2]
+    @test b1.store.displacement == zeros(3, 2)
+    @test b1.store.velocity == [1.0 1.0; 0.0 0.0; 0.0 0.0]
+    @test b1.store.velocity_half == zeros(3, 2)
+    @test b1.store.acceleration == zeros(3, 2)
+    @test b1.store.b_int[:,1:2] ≈ 1 .+ randbint[:,3:4]
+    @test b1.store.b_int[:,3:4] ≈ ones(3, 2)
+    @test b1.store.b_ext == zeros(3, 2)
+    @test b1.store.damage ≈ [2/3, 2/3]
+    @test b1.store.bond_active == [1, 0, 0, 1, 0, 0]
+    @test b1.store.n_active_bonds == [1, 1]
+
+    @test b2.store.position ≈ randpos
+    @test b2.store.displacement == zeros(3, 2)
+    @test b2.store.velocity == [0.0 0.0; 0.0 0.0; 0.0 0.0]
+    @test b2.store.velocity_half == zeros(3, 2)
+    @test b2.store.acceleration == zeros(3, 2)
+    @test b2.store.b_int[:,1:2] ≈ 1 .+ randbint[:,1:2]
+    @test b2.store.b_int[:,3:4] ≈ randbint[:,3:4]
+    @test b2.store.b_ext == zeros(3, 2)
+    @test b2.store.damage ≈ [2/3, 2/3]
+    @test b2.store.bond_active == [0, 0, 1, 0, 0, 1]
+    @test b2.store.n_active_bonds == [1, 1]
+end

test/integration/test_material.jl

@@ -0,0 +1,61 @@
+struct TestMaterial <: Peridynamics.AbstractMaterial end
+struct TestPointParameters <: Peridynamics.AbstractPointParameters
+    δ::Float64
+    rho::Float64
+    E::Float64
+    nu::Float64
+    G::Float64
+    K::Float64
+    λ::Float64
+    μ::Float64
+    Gc::Float64
+    εc::Float64
+end
+Peridynamics.point_param_type(::TestMaterial) = TestPointParameters
+Peridynamics.allowed_material_kwargs(::TestMaterial) = Peridynamics.DEFAULT_POINT_KWARGS
+function Peridynamics.get_point_params(::TestMaterial, p::Dict{Symbol,Any})
+    δ = Peridynamics.get_horizon(p)
+    rho = Peridynamics.get_density(p)
+    E, nu, G, K, λ, μ = Peridynamics.get_elastic_params(p)
+    Gc, εc = Peridynamics.get_frac_params(p, δ, K)
+    return TestPointParameters(δ, rho, E, nu, G, K, λ, μ, Gc, εc)
+end
+Peridynamics.discretization_type(::TestMaterial) = Peridynamics.BondDiscretization
+function Peridynamics.init_discretization(body::Body{TestMaterial}, args...)
+    return Peridynamics.init_bond_discretization(body, args...)
+end
+struct TestVerletStorage <: Peridynamics.AbstractStorage
+    position::Matrix{Float64}
+    displacement::Matrix{Float64}
+    velocity::Matrix{Float64}
+    velocity_half::Matrix{Float64}
+    acceleration::Matrix{Float64}
+    b_int::Matrix{Float64}
+    b_ext::Matrix{Float64}
+    damage::Vector{Float64}
+    bond_active::Vector{Bool}
+    n_active_bonds::Vector{Int}
+end
+function Peridynamics.storage_type(::TestMaterial, ::Peridynamics.VelocityVerlet)
+    return TestVerletStorage
+end
+function Peridynamics.init_storage(::Body{TestMaterial}, ::Peridynamics.VelocityVerlet,
+                                    bd::Peridynamics.BondDiscretization,
+                                    ch::Peridynamics.ChunkHandler)
+    n_loc_points = length(ch.loc_points)
+    position = copy(bd.position)
+    displacement = zeros(3, n_loc_points)
+    velocity = zeros(3, n_loc_points)
+    velocity_half = zeros(3, n_loc_points)
+    acceleration = zeros(3, n_loc_points)
+    b_int = zeros(3, length(ch.point_ids))
+    b_ext = zeros(3, n_loc_points)
+    damage = zeros(n_loc_points)
+    bond_active = ones(Bool, length(bd.bonds))
+    n_active_bonds = copy(bd.n_neighbors)
+    return TestVerletStorage(position, displacement, velocity, velocity_half,
+                                acceleration, b_int, b_ext, damage, bond_active,
+                                n_active_bonds)
+end
+Peridynamics.reads_from_halo(::Type{TestMaterial}) = (:position,)
+Peridynamics.writes_to_halo(::Type{TestMaterial}) = (:b_int,)