Granular.jl

commit 6ccf07e651d07718aa2129c1cbb56dc4975b5e45
parent bcc9daa32ee17a31b2576133f546cb19890c539e
Author: Anders Damsgaard <andersd@riseup.net>
Date:   Wed, 24 Jan 2018 11:45:50 -0500

Tightening criterion for cell mapping to ensure no overlaps

Diffstat:
M
src/packing.jl
 | 
55
++++++++++++++++++++++++++++++++++++++++++++++++++++++-
M
test/packing.jl
 | 
35
++++++++++++++++++++++-------------

2 files changed, 76 insertions(+), 14 deletions(-)
diff --git a/src/packing.jl b/src/packing.jl
@@ -308,6 +308,58 @@ function irregularPacking!(simulation::Simulation;
     end
 end
 
+export rasterPacking!
+function rasterPacking!(sim::Simulation,
+                        r_min::Real,
+                        r_max::Real;
+                        padding_factor::Real=0.1,
+                        size_distribution::String="powerlaw",
+                        size_distribution_parameter::Real=-1.8,
+                        seed::Integer=1,
+                        verbose::Bool=true)
+
+    r_rand = 0.
+    const h = .5   # disc tickness
+    const dx = r_max * 2. * (1. + padding_factor)  # cell size
+    const dx_padding = r_max * 2. * padding_factor
+    srand(seed)
+
+    const np_init = length(sim.grains)
+
+    # Generate a grid spanning the entire domain, with cell width corresponding
+    # to the largest grain to be inserted
+    const occupied = rasterMap(sim, dx)
+
+    # Add grains in unoccupied places
+    pos = zeros(2)
+    for ix=1:size(occupied, 1)
+        for iy=1:size(occupied, 2)
+
+            if occupied[ix,iy]
+                continue
+            end
+
+            if size_distribution == "powerlaw"
+                r_rand = Granular.randpower(1, size_distribution_parameter,
+                                            r_min, r_max)
+            elseif size_distribution == "uniform"
+                r_rand = rand()*(r_max - r_min) + r_min
+            end
+
+            # Determine position from grid index and sample randomly from within
+            # padding
+            pos = [ix*dx - .5*dx, iy*dx - .5*dx] .+
+                rand(2) .* dx_padding .- .5*dx_padding
+
+            addGrainCylindrical!(sim, pos, r_rand, h, verbose=false)
+
+        end
+    end
+    if verbose
+        info("Generated $(length(sim.grains) - np_init) points")
+    end
+end
+
 """
     rasterMap(sim, dx)
 
@@ -347,6 +399,7 @@ function rasterMap(sim::Simulation, dx::Real)
     min_i = 0; min_j = 0
     max_i = 0; max_j = 0
     cell_mid_point = zeros(2)
+    const dist = sqrt(2.0*(dx/2.0)^2.)
     for grain in sim.grains
         
         # Find center position in `occupied` grid
@@ -364,7 +417,7 @@ function rasterMap(sim::Simulation, dx::Real)
                 cell_mid_point = dx .* Vector{Float64}([i,j]) - 0.5 * dx
 
                 if (norm(cell_mid_point - grain.lin_pos) -
-                    grain.contact_radius < 0.5*dx)
+                    grain.contact_radius < dist)
                     occupied[i,j] = true
                 end
             end
diff --git a/test/packing.jl b/test/packing.jl
@@ -84,7 +84,7 @@ Granular.irregularPacking!(sim,
 @test length(sim.grains) > 280
 
 
-info("Testing raster-based packing algorithm")
+info("Testing raster-based mapping algorithm")
 sim = Granular.createSimulation("raster-packing1")
 sim.ocean = Granular.createRegularOceanGrid([1, 1, 1], [1., 1., 1.])
 Granular.addGrainCylindrical!(sim, [0.5, 0.5], 0.4, 1.0)
@@ -92,36 +92,45 @@ occupied = Granular.rasterMap(sim, 0.08)
 occupied_ans = Array{Bool}([
 0 0 0 0 0 0 0 0 0 0 0 0;
 0 0 0 1 1 1 1 1 1 0 0 0;
-0 0 1 1 1 1 1 1 1 1 0 0;
-0 1 1 1 1 1 1 1 1 1 1 0;
+0 0 1 1 1 1 1 1 1 1 1 0;
 0 1 1 1 1 1 1 1 1 1 1 0;
 0 1 1 1 1 1 1 1 1 1 1 1;
 0 1 1 1 1 1 1 1 1 1 1 1;
 0 1 1 1 1 1 1 1 1 1 1 1;
+0 1 1 1 1 1 1 1 1 1 1 1;
 0 1 1 1 1 1 1 1 1 1 1 0;
 0 0 1 1 1 1 1 1 1 1 1 0;
-0 0 0 1 1 1 1 1 1 1 0 0;
-0 0 0 0 0 1 1 1 0 0 0 0])
+0 0 1 1 1 1 1 1 1 1 0 0;
+0 0 0 0 1 1 1 1 0 0 0 0])
 @test occupied == occupied_ans
 Granular.addGrainCylindrical!(sim, [0.03, 0.03], 0.02, 1.0)
 occupied = Granular.rasterMap(sim, 0.08)
 occupied_ans = Array{Bool}([
 1 0 0 0 0 0 0 0 0 0 0 0;
 0 0 0 1 1 1 1 1 1 0 0 0;
-0 0 1 1 1 1 1 1 1 1 0 0;
-0 1 1 1 1 1 1 1 1 1 1 0;
+0 0 1 1 1 1 1 1 1 1 1 0;
 0 1 1 1 1 1 1 1 1 1 1 0;
 0 1 1 1 1 1 1 1 1 1 1 1;
 0 1 1 1 1 1 1 1 1 1 1 1;
 0 1 1 1 1 1 1 1 1 1 1 1;
+0 1 1 1 1 1 1 1 1 1 1 1;
 0 1 1 1 1 1 1 1 1 1 1 0;
 0 0 1 1 1 1 1 1 1 1 1 0;
-0 0 0 1 1 1 1 1 1 1 0 0;
-0 0 0 0 0 1 1 1 0 0 0 0])
+0 0 1 1 1 1 1 1 1 1 0 0;
+0 0 0 0 1 1 1 1 0 0 0 0])
 @test occupied == occupied_ans
+sim_init = deepcopy(sim)
+
+info("Testing raster-based mapping algorithm (power law GSD)")
+sim = deepcopy(sim_init)
+np_init = length(sim.grains)
+Granular.rasterPacking!(sim, 0.02, 0.04, verbose=verbose)
+@test np_init < length(sim.grains)
 
-#Granular.rasterPacking!(sim,
-                           #radius_max=.01,
-                           #radius_min=.01,
-                           #verbose=verbose)
+info("Testing raster-based mapping algorithm (uniform GSD)")
+sim = deepcopy(sim_init)
+np_init = length(sim.grains)
+Granular.rasterPacking!(sim, 0.02, 0.04, size_distribution="uniform",
+                        verbose=verbose)
+@test np_init < length(sim.grains)