Merge pull request #222 from kaipartmann/contact-dynamic-relaxation

Added support for dynamic relaxation solver and multibody simulations

Author: kaipartmann

src/time_solvers/dynamic_relaxation.jl

@@ -113,6 +113,17 @@ function init_density_matrix!(dh::AbstractThreadsBodyDataHandler, dr::DynamicRel
     return nothing
 end
 
+function init_density_matrix!(dh::AbstractThreadsMultibodyDataHandler,
+                              dr::DynamicRelaxation)
+    for body_idx in each_body_idx(dh)
+        body_dh = get_body_dh(dh, body_idx)
+        @threads :static for chunk in body_dh.chunks
+            _init_density_matrix!(chunk, dr, chunk.paramsetup)
+        end
+    end
+    return nothing
+end
+
 function init_density_matrix!(dh::AbstractMPIBodyDataHandler, dr::DynamicRelaxation)
     _init_density_matrix!(dh.chunk, dr, dh.chunk.paramsetup)
     return nothing
@@ -216,6 +227,36 @@ function relaxation_timestep!(dh::AbstractMPIBodyDataHandler,
     return nothing
 end
 
+function relaxation_timestep!(dh::AbstractThreadsMultibodyDataHandler,
+                              options::AbstractJobOptions, Δt, n)
+    t = n * Δ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)
+            chunk = body_dh.chunks[chunk_id]
+            apply_boundary_conditions!(chunk, t)
+            update_disp_and_pos!(chunk, Δt)
+        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)
+        @threads :static for chunk_id in eachindex(body_dh.chunks)
+            chunk = body_dh.chunks[chunk_id]
+            cn = calc_damping(chunk, Δt)
+            if n == 1
+                relaxation_first_step!(chunk, Δt)
+            else
+                relaxation_step!(chunk, Δt, cn)
+            end
+            export_results(body_dh, options, chunk_id, n, t)
+        end
+    end
+    return nothing
+end
+
 function calc_damping(chunk::AbstractBodyChunk, Δt::Float64)
     s = chunk.storage
     cn1 = 0.0

test/time_solvers/test_dynamic_relaxation.jl

@@ -55,3 +55,57 @@ end
     msg = String(take!(io))
     @test contains(msg, "DynamicRelaxation:\n  n_steps  1\n  Δt       1\n  Λ        1")
 end
+
+@testitem "init_density_matrix! and relaxation_timestep! for all data handlers" begin
+    import Peridynamics: init_density_matrix!, relaxation_timestep!
+    # Setup single body
+    position = [0.0 1.0 0.0 0.0;
+                0.0 0.0 0.0 1.0;
+                0.0 0.0 1.0 0.0]
+    volume = fill(1.0, 4)
+    body = Body(BBMaterial(), position, volume)
+    material!(body; horizon=1.5, rho=8e-6, E=210e3)
+    point_set!(body, :a, 1:2)
+    material!(body, :a; horizon=3, rho=7.8e-6, E=200e3)
+    forcedensity_bc!((p,t) -> 1e5 * p[1], body, :a, :x)
+    dr = DynamicRelaxation(steps=2)
+    options = Peridynamics.JobOptions(body)
+
+    # ThreadsBodyDataHandler
+    dh_threads = Peridynamics.threads_data_handler(body, dr, 1)
+    init_density_matrix!(dh_threads, dr)
+    @test all(dh_threads.chunks[1].storage.density_matrix .>= 0)
+    @test dh_threads.chunks[1].storage.position ≈ position
+    relaxation_timestep!(dh_threads, options, dr.Δt, 1)
+    relaxation_timestep!(dh_threads, options, dr.Δt, 2)
+    # Check that the density matrix is updated correctly
+    @test dh_threads.chunks[1].storage.position[1, 2] ≈ 1.05
+
+    # ThreadsMultibodyDataHandler
+    b2 = Body(BBMaterial(), position, volume)
+    material!(b2; horizon=1.5, rho=8e-6, E=210e3)
+    ms = MultibodySetup(:b1 => body, :b2 => b2)
+    dh_multibody = Peridynamics.threads_data_handler(ms, dr, 1)
+    options_multi = Peridynamics.JobOptions(ms)
+    init_density_matrix!(dh_multibody, dr)
+    for body_idx in Peridynamics.each_body_idx(dh_multibody)
+        for chunk in Peridynamics.get_body_dh(dh_multibody, body_idx).chunks
+            @test all(chunk.storage.density_matrix .>= 0)
+            @test chunk.storage.position ≈ position
+        end
+    end
+    relaxation_timestep!(dh_multibody, options_multi, dr.Δt, 1)
+    relaxation_timestep!(dh_multibody, options_multi, dr.Δt, 2)
+    # Check that the density matrix is updated correctly
+    @test dh_multibody.body_dhs[1].chunks[1].storage.position[1, 2] ≈ 1.05
+
+    # AbstractMPIBodyDataHandler
+    dh_mpi = Peridynamics.mpi_data_handler(body, dr)
+    init_density_matrix!(dh_mpi, dr)
+    @test all(dh_mpi.chunk.storage.density_matrix .>= 0)
+    @test dh_mpi.chunk.storage.position ≈ position
+    relaxation_timestep!(dh_mpi, options, dr.Δt, 1)
+    relaxation_timestep!(dh_mpi, options, dr.Δt, 2)
+    # Check that the density matrix is updated correctly
+    @test dh_mpi.chunk.storage.position[1, 2] ≈ 1.05
+end