tadd tests for disabled ice floes - Granular.jl - Julia package for granular dy… | |
git clone git://src.adamsgaard.dk/Granular.jl | |
Log | |
Files | |
Refs | |
README | |
LICENSE | |
--- | |
commit 2208ff9aec927ef82e2e88cfaed96bbd066aaf93 | |
parent fdcbb84010a27836443ff76524129e0dfa956f48 | |
Author: Anders Damsgaard <[email protected]> | |
Date: Mon, 1 May 2017 15:51:01 -0400 | |
add tests for disabled ice floes | |
Diffstat: | |
M examples/nares_strait.jl | 2 +- | |
M src/contact_search.jl | 2 +- | |
M src/datatypes.jl | 1 + | |
M src/interaction.jl | 65 +++++++++++++++++++++++++----… | |
M src/simulation.jl | 7 ++++++- | |
M test/contact-search-and-geometry.jl | 12 ++++++++++++ | |
6 files changed, 73 insertions(+), 16 deletions(-) | |
--- | |
diff --git a/examples/nares_strait.jl b/examples/nares_strait.jl | |
t@@ -10,7 +10,7 @@ L = [Lx, Lx*1.5, 1e3] | |
Ly_constriction = L[2]*.33 | |
#n = [100, 100, 2] | |
#n = [50, 50, 2] | |
-n = [25, 25, 2] | |
+n = [6, 6, 2] | |
sim.ocean = SeaIce.createRegularOceanGrid(n, L, name="poiseuille_flow") | |
sim.ocean.v[:, :, 1, 1] = 1e-8*((sim.ocean.xq - Lx/2.).^2 - Lx^2./4.) | |
diff --git a/src/contact_search.jl b/src/contact_search.jl | |
t@@ -122,7 +122,7 @@ function checkAndAddContact!(simulation::Simulation, i::In… | |
if overlap_ij < 0.0 | |
push!(simulation.contact_pairs, [i, j]) | |
push!(simulation.overlaps, overlap_ij*position_ij/norm(position_ij… | |
- #push!(simulation.contact_parallel_displacement, zeros(2)) | |
+ push!(simulation.contact_parallel_displacement, zeros(2)) | |
end | |
end | |
end | |
diff --git a/src/datatypes.jl b/src/datatypes.jl | |
t@@ -173,6 +173,7 @@ type Simulation | |
ice_floes::Array{IceFloeCylindrical, 1} | |
contact_pairs::Array{Array{Int, 1}, 1} | |
overlaps::Array{Array{Float64, 1}, 1} | |
+ contact_parallel_displacement::Array{Array{Float64, 1}, 1} | |
ocean::Ocean | |
end | |
diff --git a/src/interaction.jl b/src/interaction.jl | |
t@@ -10,8 +10,10 @@ function interact!(simulation::Simulation) | |
while !isempty(simulation.contact_pairs) | |
contact_pair = pop!(simulation.contact_pairs) | |
overlap_vector = pop!(simulation.overlaps) | |
+ contact_parallel_displacement = | |
+ pop!(simulation.contact_parallel_displacement) | |
interactIceFloes!(simulation, contact_pair[1], contact_pair[2], | |
- overlap_vector) | |
+ overlap_vector, contact_parallel_displacement) | |
end | |
end | |
t@@ -23,28 +25,36 @@ function adds the compressive force of the interaction to … | |
""" | |
function interactIceFloes!(simulation::Simulation, | |
i::Integer, j::Integer, | |
- overlap_vector::Array{Float64, 1}; | |
- contact_normal::String = "LinearElastic") | |
+ overlap_vector::Array{Float64, 1}, | |
+ contact_parallel_displacement::Array{Float64, 1}; | |
+ contact_normal::String = "LinearElasticNoTangential… | |
force = zeros(2) | |
if contact_normal == "None" | |
- # do nothing | |
+ return nothing | |
- elseif contact_normal == "LinearElastic" | |
- force = interactNormalLinearElastic(simulation, i, j, overlap_vector) | |
+ elseif contact_normal == "LinearElasticNoTangential" || | |
+ contact_normal == "LinearElastic" | |
+ force_n = interactNormalLinearElastic(simulation, i, j, overlap_vector) | |
+ | |
+ if contact_normal == "LinearElastic" | |
+ force_t, torque = interactTangentialLinearElastic(simulation, i, j, | |
+ overlap_vector, | |
+ contact_parallel… | |
+ end | |
else | |
error("Unknown contact_normal interaction model '$contact_normal'") | |
end | |
- simulation.ice_floes[i].force += force; | |
- simulation.ice_floes[j].force -= force; | |
+ simulation.ice_floes[i].force += force_n; | |
+ simulation.ice_floes[j].force -= force_n; | |
simulation.ice_floes[i].pressure += | |
- norm(force)/simulation.ice_floes[i].side_surface_area; | |
+ norm(force_n)/simulation.ice_floes[i].side_surface_area; | |
simulation.ice_floes[j].pressure += | |
- norm(force)/simulation.ice_floes[j].side_surface_area; | |
+ norm(force_n)/simulation.ice_floes[j].side_surface_area; | |
end | |
export interactNormalLinearElastic | |
t@@ -56,9 +66,38 @@ function interactNormalLinearElastic(simulation::Simulation, | |
i::Integer, j::Integer, | |
overlap_vector::vector) | |
- k_n_i = simulation.ice_floes[i].contact_stiffness_normal | |
- k_n_j = simulation.ice_floes[j].contact_stiffness_normal | |
- k_n_harmonic_mean = 2.*k_n_i*k_n_j/(k_n_i + k_n_j) | |
+ k_n_harmonic_mean = | |
+ harmonicMean(simulation.ice_floes[i].contact_stiffness_normal, | |
+ simulation.ice_floes[j].contact_stiffness_normal) | |
return k_n_harmonic_mean * overlap_vector | |
end | |
+ | |
+export interactTangentialLinearElastic | |
+""" | |
+Resolves linear-elastic interaction between two ice floes in the | |
+contact-parallel (tangential) direction. | |
+""" | |
+function interactTangentialLinearElastic(simulation::Simulation, | |
+ i::Integer, j::Integer, | |
+ overlap_vector::vector, | |
+ contact_parallel_displacement::vector) | |
+ | |
+ k_t_harmonic_mean = | |
+ harmonicMean(simulation.ice_floes[i].contact_stiffness_tangential, | |
+ simulation.ice_floes[j].contact_stiffness_tangential) | |
+ | |
+ # correct displacement for contact rotation | |
+ n = overlap_vector/norm(overlap_vector) | |
+ contact_parallel_displacement -= (n * dot(n, contact_parallel_displacement… | |
+ | |
+ | |
+ force_t = k_t_harmonic_mean * contact_parallel_displacement | |
+ | |
+ return force, torque, contact_parallel_displacement | |
+ | |
+end | |
+ | |
+function harmonicMean(a::Any, b::Any) | |
+ return 2.*a*b/(a + b) | |
+end | |
diff --git a/src/simulation.jl b/src/simulation.jl | |
t@@ -12,6 +12,8 @@ export createSimulation | |
ice_floes=Array{IceFloeCylindrical, 1}[], | |
contact_pairs=Array{Int64, 1}[], | |
overlaps=Array{Array{Float64, 1}, 1}[], | |
+ contact_parallel_displacement= | |
+ Array{Array{Float64, 1}, 1}[]) | |
Create a simulation object containing all relevant variables such as temporal | |
parameters, and lists of ice floes and contacts. | |
t@@ -30,6 +32,8 @@ function createSimulation(;id::String="unnamed", | |
ice_floes=Array{IceFloeCylindrical, 1}[], | |
contact_pairs=Array{Int64, 1}[], | |
overlaps=Array{Array{Float64, 1}, 1}[], | |
+ contact_parallel_displacement= | |
+ Array{Array{Float64, 1}, 1}[], | |
ocean::Ocean=createEmptyOcean()) | |
return Simulation(id, | |
t@@ -43,6 +47,7 @@ function createSimulation(;id::String="unnamed", | |
ice_floes, | |
contact_pairs, | |
overlaps, | |
+ contact_parallel_displacement, | |
ocean) | |
end | |
t@@ -165,7 +170,7 @@ function disableIceFloe!(simulation::Simulation, i::Intege… | |
error("Index must be greater than 0 (i = $i)") | |
end | |
- simulation.ice_floes[i].enable = false | |
+ simulation.ice_floes[i].enabled = false | |
end | |
export zeroForcesAndTorques! | |
diff --git a/test/contact-search-and-geometry.jl b/test/contact-search-and-geom… | |
t@@ -54,6 +54,18 @@ SeaIce.findContacts!(sim) | |
@test 0 == length(sim.overlaps) | |
@test 0 == length(sim.contact_pairs) | |
+sim = deepcopy(sim_copy) | |
+SeaIce.disableIceFloe!(sim, 1) | |
+SeaIce.findContacts!(sim) | |
+@test 0 == length(sim.overlaps) | |
+@test 0 == length(sim.contact_pairs) | |
+ | |
+sim = deepcopy(sim_copy) | |
+SeaIce.disableIceFloe!(sim, 1) | |
+SeaIce.disableIceFloe!(sim, 2) | |
+SeaIce.findContacts!(sim) | |
+@test 0 == length(sim.overlaps) | |
+@test 0 == length(sim.contact_pairs) | |
info("Testing if interact(...) removes contacts correctly") | |
sim = deepcopy(sim_copy) |