commit bcc9daa32ee17a31b2576133f546cb19890c539e
parent f7da3ebd8db4dbfa61ca5c0ae22391cff43e8644
Author: Anders Damsgaard <andersd@riseup.net>
Date: Wed, 24 Jan 2018 11:03:42 -0500
Add function to produce rasterized map of grain coverage
Diffstat:
| M | src/packing.jl | | | 65 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| M | test/packing.jl | | | 43 | +++++++++++++++++++++++++++++++++++++++++++ |
2 files changed, 108 insertions(+), 0 deletions(-)
diff --git a/src/packing.jl b/src/packing.jl
@@ -307,3 +307,68 @@ function irregularPacking!(simulation::Simulation;
info("Generated $(length(simulation.grains) - np_init) points")
end
end
+
+"""
+ rasterMap(sim, dx)
+
+Returns a rasterized map of grain extent in the domain with length `L` and a
+pixel size of `dx`. The function will return a map of `Bool` type with size
+`floor.(L./dx)`.
+
+* Arguments
+- `sim::Simulation`: simulation object containing the grains.
+- `dx::Real`: pixel size to use for the raster map.
+
+"""
+function rasterMap(sim::Simulation, dx::Real)
+
+ # Use existing `grid` (ocean or atmosphere) for contact search
+ if typeof(sim.ocean.input_file) != Bool
+ grid = sim.ocean
+ elseif typeof(sim.atmosphere.input_file) != Bool
+ grid = sim.atmosphere
+ else
+ error("rasterMap(...) requires an ocean or atmosphere grid")
+ end
+ # save grid boundaries
+ if grid.regular_grid
+ L = grid.L[1:2]
+ origo = grid.origo
+ else
+ sw, se, ne, nw = getGridCornerCoordinates(grid.xq, grid.yq)
+ L = [se[1] - sw[1], nw[2] - sw[2]]
+ origo = [sw[1], sw[2]]
+ end
+ const dims = floor.(L./dx)
+ occupied = zeros(Bool, dims[1], dims[2])
+
+ # Loop over existing grains and mark their extent in the `occupied` array
+ i = 0; j = 0
+ min_i = 0; min_j = 0
+ max_i = 0; max_j = 0
+ cell_mid_point = zeros(2)
+ for grain in sim.grains
+
+ # Find center position in `occupied` grid
+ #i, j = Int.(floor.((grain.lin_pos - origo) ./ dx)) + [1,1]
+
+ # Find corner indexes for box spanning the grain
+ min_i, min_j = Int.(floor.((grain.lin_pos - origo -
+ grain.contact_radius) ./ dx)) + [1,1]
+ max_i, max_j = Int.(floor.((grain.lin_pos - origo +
+ grain.contact_radius) ./ dx)) + [1,1]
+
+ # For each cell in box, check if the grain is contained
+ for i = min_i:max_i
+ for j = min_j:max_j
+ 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)
+ occupied[i,j] = true
+ end
+ end
+ end
+ end
+ return occupied
+end
diff --git a/test/packing.jl b/test/packing.jl
@@ -82,3 +82,46 @@ Granular.irregularPacking!(sim,
plot_during_packing=plot_packings,
verbose=verbose)
@test length(sim.grains) > 280
+
+
+info("Testing raster-based packing 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)
+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 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 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])
+@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 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 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])
+@test occupied == occupied_ans
+
+#Granular.rasterPacking!(sim,
+ #radius_max=.01,
+ #radius_min=.01,
+ #verbose=verbose)
+
