diff --git a/Project.toml b/Project.toml index 1af551ca..f84e33dd 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Peridynamics" uuid = "4dc47793-80f3-4232-b30e-ca78ca9d621b" authors = ["Kai Friebertshäuser"] -version = "0.1.0" +version = "0.1.1" [deps] AbaqusReader = "bc6b9049-e460-56d6-94b4-a597b2c0390d" diff --git a/docs/src/logo.md b/docs/src/logo.md index f3730815..c83556ff 100644 --- a/docs/src/logo.md +++ b/docs/src/logo.md @@ -75,7 +75,7 @@ matₛ = BondBasedMaterial(; All material points of the spheres have a initial velocity of $-20\, \mathrm{m} \, \mathrm{s}^{-1}$ in $z$-direction. ```julia -icₛ = [VelocityIC(-20.0, 1:(pcₛ₁.n_points), 3)] +icₛ = [VelocityIC(-20.0, 1:pcₛ₁.n_points, 3)] ``` For the contact analysis, every body needs to be specified with a [`BodySetup`](@ref). Time can be saved by using only one sphere for the calculation of the stable time step. diff --git a/examples/Logo.jl b/examples/Logo.jl index 155f6b95..403bb0be 100644 --- a/examples/Logo.jl +++ b/examples/Logo.jl @@ -48,7 +48,7 @@ matₛ = BondBasedMaterial(; ) ##-- 4. DEFINE INITIAL CONDITIONS ------------------------------------------------------ -icₛ = [VelocityIC(-20.0, 1:(pcₛ₁.n_points), 3)] +icₛ = [VelocityIC(-20.0, 1:pcₛ₁.n_points, 3)] ##-- 5. DEFINE THE BODY-SETUP ---------------------------------------------------------- plate = BodySetup(pcₚ, matₚ) diff --git a/src/BondBased.jl b/src/BondBased.jl index 7e0b2329..20a9b8a2 100644 --- a/src/BondBased.jl +++ b/src/BondBased.jl @@ -99,7 +99,7 @@ function BondBasedBody(mat::BondBasedMaterial, pc::PointCloud) owned_bonds = defaultdist(n_bonds, n_threads) bond_failure = ones(Int, n_bonds) n_active_family_members = zeros(Int, (n_points, n_threads)) - @threads :static for tid in 1:n_threads + @threads for tid in 1:n_threads for i in owned_points[tid] n_active_family_members[i, tid] = n_family_members[i] end @@ -138,8 +138,7 @@ create_simmodel(mat::BondBasedMaterial, pc::PointCloud) = BondBasedBody(mat, pc) function compute_forcedensity!(body::BondBasedBody, mat::BondBasedMaterial) body.b_int .= 0.0 body.n_active_family_members .= 0 - @threads :static for _ in 1:(body.n_threads) - tid = threadid() + @threads for tid in 1:body.n_threads for current_bond in body.owned_bonds[tid] (i, j, ξ, failureflag) = body.bond_data[current_bond] ϑx = body.position[1, j] - body.position[1, i] diff --git a/src/ContactBase.jl b/src/ContactBase.jl index 9daa8cf6..da47c1ce 100644 --- a/src/ContactBase.jl +++ b/src/ContactBase.jl @@ -97,7 +97,7 @@ function submit(sim::PDContactAnalysis) log_describe_sim(sim) bodies = [create_simmodel(ps.mat, ps.pc) for ps in sim.body_setup] log_describe_interactions(bodies, sim.es) - for i in 1:(sim.n_bodies) + for i in 1:sim.n_bodies for precrack in sim.body_setup[i].precracks define_precrack!(bodies[i], precrack) end @@ -105,7 +105,7 @@ function submit(sim::PDContactAnalysis) end if sim.td.Δt < 0.0 _Δt = Vector{Float64}(undef, 0) - for i in 1:(sim.n_bodies) + for i in 1:sim.n_bodies if sim.body_setup[i].calc_timestep push!( _Δt, @@ -121,11 +121,11 @@ function submit(sim::PDContactAnalysis) sim.td.Δt = minimum(_Δt) end log_simsetup(sim) - for i in 1:(sim.n_bodies) + for i in 1:sim.n_bodies apply_ics!(bodies[i], sim.body_setup[i].ics) end if sim.es.exportflag - for i in 1:(sim.n_bodies) + for i in 1:sim.n_bodies export_results(bodies[i], string(sim.es.resultfile_prefix, "_b", i), 0, 0.0) end end @@ -137,7 +137,7 @@ end function log_describe_sim(sim::PDContactAnalysis) msg = "Peridynamic multibody contact analysis: " * sim.name * "\nMaterial parameters:\n" - for i in 1:(sim.n_bodies) + for i in 1:sim.n_bodies msg *= @sprintf("Body %d:\n", i) msg *= @sprintf( " - Number of material points [-]: %30g\n", @@ -209,9 +209,9 @@ function velocity_verlet!(body::Vector{<:AbstractPDBody}, sim::PDContactAnalysis enabled=!is_logging(stderr), ) Δt½ = 0.5 * sim.td.Δt - for t in 1:(sim.td.n_timesteps) + for t in 1:sim.td.n_timesteps time = t * sim.td.Δt - for i in 1:(sim.n_bodies) + for i in 1:sim.n_bodies update_velhalf!(body[i], Δt½) apply_bcs!(body[i], sim.body_setup[i].bcs, time) update_disp_and_position!(body[i], sim.td.Δt) @@ -219,11 +219,11 @@ function velocity_verlet!(body::Vector{<:AbstractPDBody}, sim::PDContactAnalysis calc_damage!(body[i]) end compute_contactforcedensity!(body, sim) - for i in 1:(sim.n_bodies) + for i in 1:sim.n_bodies compute_equation_of_motion!(body[i], Δt½, sim.body_setup[i].mat.rho) end if mod(t, sim.es.exportfreq) == 0 - for i in 1:(sim.n_bodies) + for i in 1:sim.n_bodies export_results(body[i], string(sim.es.resultfile_prefix, "_b", i), t, time) end end @@ -242,10 +242,9 @@ function compute_contactforcedensity!( body_b = bodies[jj] r = contact.search_radius C = 9 * contact.spring_constant / (π * sim.body_setup[ii].mat.δ^4) - @threads :static for _ in 1:nthreads() - tid = threadid() + @threads for tid in 1:nthreads() for i in body_a.owned_points[tid] - for j in 1:(body_b.n_points) + for j in 1:body_b.n_points ϑx = body_b.position[1, j] - body_a.position[1, i] ϑy = body_b.position[2, j] - body_a.position[2, i] ϑz = body_b.position[3, j] - body_a.position[3, i] diff --git a/src/PeridynamicsBase.jl b/src/PeridynamicsBase.jl index 0cc5da78..40aac6d6 100644 --- a/src/PeridynamicsBase.jl +++ b/src/PeridynamicsBase.jl @@ -431,8 +431,7 @@ function find_bonds(pc::PointCloud, δ::Float64, owned_points::Vector{UnitRange{ color=:normal, enabled=!is_logging(stderr), ) - @threads :static for _ in 1:n_threads - tid = threadid() + @threads for tid in 1:n_threads local_bond_data = Vector{Tuple{Int,Int,Float64,Bool}}(undef, 0) sizehint!(local_bond_data, pc.n_points * 500) idst = 0.0 @@ -440,7 +439,7 @@ function find_bonds(pc::PointCloud, δ::Float64, owned_points::Vector{UnitRange{ fail = false for a in owned_points[tid] num = 0 - for i in 1:(pc.n_points) + for i in 1:pc.n_points if a !== i idst = sqrt( (pc.position[1, i] - pc.position[1, a])^2 + @@ -475,8 +474,7 @@ function find_unique_bonds(pc::PointCloud, δ::Float64, owned_points::Vector{Uni color=:normal, enabled=!is_logging(stderr), ) - @threads :static for _ in 1:n_threads - tid = threadid() + @threads for tid in 1:n_threads local_bonds_data = Vector{Tuple{Int,Int,Float64,Bool}}(undef, 0) sizehint!(local_bonds_data, pc.n_points * 500) idst = 0.0 @@ -484,7 +482,7 @@ function find_unique_bonds(pc::PointCloud, δ::Float64, owned_points::Vector{Uni fail = false for a in owned_points[tid] num = 0 - for i in 1:(pc.n_points) + for i in 1:pc.n_points if a !== i idst = sqrt( (pc.position[1, i] - pc.position[1, a])^2 + @@ -531,7 +529,7 @@ function defaultdist(sz::Int, nc::Int) return grid else sidx = [1:(sz + 1);] - grid = fill(0:0, nc) + grid = fill(0:-1, nc) for i in 1:sz grid[i] = sidx[i]:(sidx[i + 1] - 1) end @@ -737,8 +735,7 @@ function log_closesim(part::AbstractPDBody, timingsummary::NamedTuple, es::Expor end function define_precrack!(body::AbstractPDBody, precrack::PreCrack) - @threads :static for _ in 1:(body.n_threads) - tid = threadid() + @threads for tid in 1:body.n_threads (@view body.n_active_family_members[:, tid]) .= 0 for current_one_ni in body.owned_bonds[tid] i, j, _, _ = body.bond_data[current_one_ni] @@ -760,8 +757,7 @@ end function calc_stable_timestep(body::AbstractPDBody, rho::Float64, K::Float64, δ::Float64) _Δt = zeros(Float64, body.n_threads) - @threads :static for _ in 1:(body.n_threads) - tid = threadid() + @threads for tid in 1:body.n_threads timesteps = zeros(Float64, body.n_points) dtsum = zeros(Float64, (body.n_points, body.n_threads)) for current_one_ni in body.owned_bonds[tid] @@ -796,8 +792,7 @@ function calc_stable_user_timestep(pc::PointCloud, mat::AbstractPDMaterial, Sf:: bond_data, _ = find_bonds(pc, mat.δ, owned_points) owned_bonds = defaultdist(length(bond_data), nthreads()) _Δt = zeros(Float64, nthreads()) - @threads :static for _ in 1:nthreads() - tid = threadid() + @threads for tid in 1:nthreads() timesteps = zeros(Float64, pc.n_points) dtsum = zeros(Float64, (pc.n_points, nthreads())) for current_one_ni in owned_bonds[tid] @@ -815,8 +810,8 @@ function calc_stable_user_timestep(pc::PointCloud, mat::AbstractPDMaterial, Sf:: end function calc_damage!(body::AbstractPDBody) - @threads :static for _ in 1:(body.n_threads) - for a in body.owned_points[threadid()] + @threads for tid in 1:body.n_threads + for a in body.owned_points[tid] body.n_active_family_members[a, 1] = sum( @view body.n_active_family_members[a, :] ) @@ -897,7 +892,7 @@ end function export_vtk(body::AbstractPDBody, expfile::String, timestep::Int, time::Float64) filename = expfile * "_t" * string(timestep) - cells = [MeshCell(VTKCellTypes.VTK_VERTEX, (j,)) for j in 1:(body.n_points)] + cells = [MeshCell(VTKCellTypes.VTK_VERTEX, (j,)) for j in 1:body.n_points] vtkfile = vtk_grid(filename, body.position, cells) vtkfile["Damage", VTKPointData()] = body.damage vtkfile["Force density", VTKPointData()] = @views body.b_int[:, :, 1] @@ -923,7 +918,7 @@ function velocity_verlet!(body::AbstractPDBody, sim::PDSingleBodyAnalysis) enabled=!is_logging(stderr), ) Δt½ = 0.5 * sim.td.Δt - for t in 1:(sim.td.n_timesteps) + for t in 1:sim.td.n_timesteps time = t * sim.td.Δt update_velhalf!(body, Δt½) apply_bcs!(body, sim.bcs, time) @@ -956,7 +951,7 @@ function dynamic_relaxation_finite_difference!( color=:normal, enabled=!is_logging(stderr), ) - for t in 1:(sim.td.n_timesteps) + for t in 1:sim.td.n_timesteps time = t * sim.td.Δt apply_bcs!(body, sim.bcs, time) compute_forcedensity!(body, sim.mat) @@ -989,7 +984,7 @@ function calc_damping( ) cn1 = 0.0 cn2 = 0.0 - for i in 1:(body.n_points) + for i in 1:body.n_points for d in 1:3 body.b_int[d, i, 1] = sum(@view body.b_int[d, i, :]) if velocity_half_old[d, i] !== 0.0 @@ -1022,7 +1017,7 @@ function finite_difference_first_step!( b_int_old::Matrix{Float64}, Δt::Float64, ) - @threads for i in 1:(body.n_points) + @threads for i in 1:body.n_points for d in 1:3 body.velocity_half[d, i] = 0.5 * Δt / damping_matrix[d, i] * (body.b_int[d, i, 1] + body.b_ext[d, i]) @@ -1044,7 +1039,7 @@ function finite_difference!( Δt::Float64, cn::Float64, ) - @threads for i in 1:(body.n_points) + @threads for i in 1:body.n_points for d in 1:3 body.velocity_half[d, i] = ( @@ -1063,7 +1058,7 @@ function finite_difference!( end function update_velhalf!(body::AbstractPDBody, Δt½::Float64) - @threads for i in 1:(body.n_points) + @threads for i in 1:body.n_points body.velocity_half[1, i] = body.velocity[1, i] + body.acceleration[1, i] * Δt½ body.velocity_half[2, i] = body.velocity[2, i] + body.acceleration[2, i] * Δt½ body.velocity_half[3, i] = body.velocity[3, i] + body.acceleration[3, i] * Δt½ @@ -1072,7 +1067,7 @@ function update_velhalf!(body::AbstractPDBody, Δt½::Float64) end function update_disp_and_position!(body::AbstractPDBody, Δt::Float64) - @threads for i in 1:(body.n_points) + @threads for i in 1:body.n_points body.displacement[1, i] += body.velocity_half[1, i] * Δt body.displacement[2, i] += body.velocity_half[2, i] * Δt body.displacement[3, i] += body.velocity_half[3, i] * Δt @@ -1084,7 +1079,7 @@ function update_disp_and_position!(body::AbstractPDBody, Δt::Float64) end function compute_equation_of_motion!(body::AbstractPDBody, Δt½::Float64, rho::Float64) - @threads for i in 1:(body.n_points) + @threads for i in 1:body.n_points body.b_int[1, i, 1] = sum(@view body.b_int[1, i, :]) body.b_int[2, i, 1] = sum(@view body.b_int[2, i, :]) body.b_int[3, i, 1] = sum(@view body.b_int[3, i, :]) diff --git a/test/bondbased.jl b/test/bondbased.jl index bb2c8406..c083f609 100644 --- a/test/bondbased.jl +++ b/test/bondbased.jl @@ -49,7 +49,11 @@ if Threads.nthreads() <= 2 @test body.n_threads == Threads.nthreads() @test body.unique_bonds == true @test body.b_int == zeros(3, n_points, Threads.nthreads()) - @test body.n_active_family_members == [1; 1;;] + if Threads.nthreads() == 1 + @test body.n_active_family_members == [1; 1;;] + elseif Threads.nthreads() == 2 + @test body.n_active_family_members == [1 0; 0 1] + end @test body.bond_failure == [1] # Boundary Condition: @@ -57,7 +61,11 @@ if Threads.nthreads() <= 2 body.position[1, 2] = 1.0015 Peridynamics.compute_forcedensity!(body, mat) - @test body.n_active_family_members == [0; 0;;] + if Threads.nthreads() == 1 + @test body.n_active_family_members == [0; 0;;] + elseif Threads.nthreads() == 2 + @test body.n_active_family_members == [0 0; 0 0] + end @test body.bond_failure == [0] b¹² = [ @@ -65,8 +73,8 @@ if Threads.nthreads() <= 2 0.0 0.0 ] - Threads.@threads for _ in 1:Threads.nthreads() - for i in body.owned_points[Threads.threadid()] + Threads.@threads for tid in 1:Threads.nthreads() + for i in body.owned_points[tid] body.b_int[1,i,1] = sum(@view body.b_int[1,i,:]) body.b_int[2,i,1] = sum(@view body.b_int[2,i,:]) body.b_int[3,i,1] = sum(@view body.b_int[3,i,:]) diff --git a/test/bonds.jl b/test/bonds.jl index a5645598..06f3c492 100644 --- a/test/bonds.jl +++ b/test/bonds.jl @@ -37,14 +37,22 @@ if Threads.nthreads() <= 2 body = Peridynamics.create_simmodel(mat, pc) precrack = PreCrack([1], [2]) - @test body.n_active_family_members == [1; 1;;] + if Threads.nthreads() == 1 + @test body.n_active_family_members == [1; 1;;] + elseif Threads.nthreads() == 2 + @test body.n_active_family_members == [1 0; 0 1] + end @test body.bond_failure == [1] @test body.damage == [0, 0] Peridynamics.define_precrack!(body, precrack) Peridynamics.calc_damage!(body) - @test body.n_active_family_members == [0; 0;;] + if Threads.nthreads() == 1 + @test body.n_active_family_members == [0; 0;;] + elseif Threads.nthreads() == 2 + @test body.n_active_family_members == [0 0; 0 0] + end @test body.bond_failure == [0] @test body.damage == [1, 1] else diff --git a/test/timeintegration.jl b/test/timeintegration.jl index 0309e369..a7c7fbb1 100644 --- a/test/timeintegration.jl +++ b/test/timeintegration.jl @@ -91,7 +91,7 @@ if Threads.nthreads() <= 2 @test body.velocity_half == [-0.09375 0.09375; 0.0 0.0; 0.0 0.0] @test body.displacement == [-0.09375 0.09375; 0.0 0.0; 0.0 0.0] @test body.position == [-0.09375 1.09375; 0.0 0.0; 0.0 0.0] - @test body.b_int == zeros(3, 2, 1) + @test body.b_int == zeros(3, 2, Threads.nthreads()) @test body.acceleration == zeros(3, 2) @test body.velocity == [-0.046875 0.046875; 0.0 0.0; 0.0 0.0] else diff --git a/test/utilities.jl b/test/utilities.jl index 2a86d0c4..ab17b0c1 100644 --- a/test/utilities.jl +++ b/test/utilities.jl @@ -6,4 +6,4 @@ using Test @test Peridynamics.defaultdist(4, 2) == [1:2, 3:4] @test Peridynamics.defaultdist(3, 2) == [1:2, 3:3] -@test Peridynamics.defaultdist(2, 3) == [1:1, 2:2, 0:0] +@test Peridynamics.defaultdist(2, 3) == [1:1, 2:2, 0:-1]