Granular.jl

Julia package for granular dynamics simulation
git clone git://src.adamsgaard.dk/Granular.jl
Log | Files | Refs | README | LICENSE

commit eeb31284d08583ecaa038ef8fb55d77692c48e71
parent a79f1155bfb16f8da74ac94efe92e7d9eff79fb3
Author: Anders Damsgaard <andersd@riseup.net>
Date:   Sun, 30 Apr 2017 14:25:22 -0400

add cell-based contact search

Diffstat:
Msrc/contact_search.jl | 47+++++++++++++++++++++++++++++++++++++++++++++++
Msrc/ocean.jl | 6+++---
2 files changed, 50 insertions(+), 3 deletions(-)

diff --git a/src/contact_search.jl b/src/contact_search.jl @@ -47,6 +47,8 @@ end export findContactsAllToAll! """ + findContactsAllToAll!(simulation) + Perform an O(n^2) all-to-all contact search between all ice floes in the `simulation` object. Contacts between fixed ice floes are ignored. """ @@ -77,3 +79,48 @@ function findContactsAllToAll!(simulation::Simulation) end end end + +export findContactsOceanGrid! +""" + findContactsOceanGrid!(simulation) + +Perform an O(n*log(n)) cell-based contact search between all ice floes in the +`simulation` object. Contacts between fixed ice floes are ignored. +""" +function findContactsOceanGrid!(simulation::Simulation) + + for idx_i = 1:length(simulation.ice_floes) + + grid_pos = simulation.ice_floes[idx_i].ocean_grid_pos + nx, ny = size(simulation.ocean.xh) + + for i=(grid_pos[1] - 1):(grid_pos[1] + 1) + for j=(grid_pos[2] - 1):(grid_pos[2] + 1) + + # only check for contacts within grid boundaries + if i < 1 || i > nx || j < 1 || j > ny + continue + end + + for idx_j in simulation.ocean.ice_floe_list[i, j] + + if idx_i < idx_j + + # Inter-grain position vector and grain overlap + position_ij = interIceFloePositionVector(simulation, + idx_i, idx_j) + overlap_ij = findOverlap(simulation, idx_i, idx_j, + position_ij) + + # Check if grains overlap (overlap when negative) + if overlap_ij < 0.0 + push!(simulation.contact_pairs, [idx_i, idx_j]) + push!(simulation.overlaps, + overlap_ij*position_ij/norm(position_ij)) + end + end + end + end + end + end +end diff --git a/src/ocean.jl b/src/ocean.jl @@ -279,12 +279,12 @@ function applyOceanDragToIceFloe!(ice_floe::IceFloeCylindrical, freeboard = .1*ice_floe.thickness # height above water rho_o = 1000. # ocean density draft = ice_floe.thickness - freeboard # height of submerged thickness - c_o_v = 123 # ocean drag coefficient, vertical - c_o_h = 123 # ocean drag coefficient, horizontal + c_o_v = .85 # ocean drag coefficient, vertical, Hunke and Comeau 2011 + c_o_h = 5e-4 # ocean drag coefficient, horizontal length = ice_floe.areal_radius*2. width = ice_floe.areal_radius*2. ice_floe.force += rho_o * (.5*c_o_v*width*draft*freeboard + c_o_h*length*width) * - ([u, v] - ice_floe.vel)*norm([u, v] - ice_floe.vel) + ([u, v] - ice_floe.lin_vel)*norm([u, v] - ice_floe.lin_vel) end