packing.jl (18774B)
1 ## Functions for creating grain packings 2 import Random 3 using LinearAlgebra 4 using Random 5 6 export regularPacking! 7 """ 8 9 regularPacking!(simulation, n, r_min, r_max[, tiling, padding_factor, 10 size_distribution, size_distribution_parameter, seed]) 11 12 Create a grid-based regular packing with grain numbers along each axis specified 13 by the `n` vector. 14 15 # Arguments 16 * `simulation::Simulation`: simulation object where the grains are inserted, 17 preferably not containing prior grains. 18 * `n::Vector{Integer}`: 2-element vector determining number of grains along the 19 `x` and `y` axes. 20 * `r_min::Real`: minimum desired grain radius. 21 * `r_max::Real`: maximum desired grain radius. 22 * `tiling::String`: the packing method to use, valid values are `"square"` 23 (default) and `"triangular"` (see 24 [Wikipedia](https://en.wikipedia.org/wiki/Circle_packing#Uniform_packings)). 25 * `padding_factor::Real`: percentage-wise padding around each grain to allow for 26 random perturbations to grain position (default = 0.0). 27 * `origo::Vector{Real}`: spatial offset for the packing (default `[0.0, 0.0]`). 28 * `size_distribution::String`: grain-size distribution to sample. Valid values 29 are "powerlaw" and "uniform". 30 * `size_distribution_parameter::Real`: parameter to pass to the grain-size 31 distribution generating function. 32 * `seed::Integer`: seed value to the pseudo-random number generator. 33 """ 34 function regularPacking!(simulation::Simulation, 35 n::Vector{Int}, 36 r_min::Real, 37 r_max::Real; 38 tiling::String="square", 39 padding_factor::Real=0.0, 40 origo::Vector{Float64}=[0.0, 0.0], 41 size_distribution::String="powerlaw", 42 size_distribution_parameter::Real=-1.8, 43 seed::Integer=1) 44 45 r_rand = 0. 46 pos = zeros(2) 47 h = .5 # disc tickness 48 Random.seed!(seed) 49 50 if tiling == "square" 51 dx = r_max * 2. * (1. + padding_factor) # cell size 52 dx_padding = r_max * 2. * padding_factor 53 for iy in 1:n[2] 54 for ix in 1:n[1] 55 56 if size_distribution == "powerlaw" 57 r_rand = Granular.randpower(1, 58 size_distribution_parameter, 59 r_min, r_max) 60 elseif size_distribution == "uniform" 61 r_rand = rand()*(r_max - r_min) + r_min 62 end 63 64 # Determine position from grid index and sample randomly from 65 # within padding 66 pos = [ix*dx - .5*dx, iy*dx - .5*dx] .+ 67 rand(2) .* dx_padding .- .5*dx_padding .+ origo 68 69 addGrainCylindrical!(simulation, pos, r_rand, h, verbose=false) 70 end 71 end 72 73 elseif tiling == "triangular" 74 75 dx = r_max * 2. * (1. + padding_factor) # cell size 76 dy = r_max * 2. * (1. + padding_factor) * sin(π/3) 77 dx_padding = r_max * 2. * padding_factor 78 for iy in 1:n[2] 79 for ix in 1:n[1] 80 81 if size_distribution == "powerlaw" 82 r_rand = Granular.randpower(1, 83 size_distribution_parameter, 84 r_min, r_max) 85 elseif size_distribution == "uniform" 86 r_rand = rand()*(r_max - r_min) + r_min 87 end 88 89 # Determine position from grid index and sample randomly from 90 # within padding 91 if iy%2 == 0 92 pos = [ix*dx - .5*dx, (iy - 1)*dy + .5*dx] .+ 93 rand(2) .* dx_padding .- .5*dx_padding .+ origo 94 else 95 pos = [ix*dx, (iy - 1)*dy + .5*dx] .+ 96 rand(2) .* dx_padding .- .5*dx_padding .+ origo 97 end 98 99 addGrainCylindrical!(simulation, pos, r_rand, h, verbose=false) 100 end 101 end 102 103 else 104 error("tiling method '$tiling' not understood") 105 end 106 107 end 108 109 """ 110 Return random point in spherical annulus between `(r_i + r_j)` and `2.*(r_i + 111 r_j)` around `x_i`. Note: there is slightly higher point density towards (r_i + 112 r_j). 113 """ 114 function generateNeighboringPoint(x_i::Vector, r_i::Real, 115 r_max::Real, r_min::Real; 116 padding::Real=0.) 117 118 if padding > 0. 119 r_j = r_min + (rand()*0.5 + 0.5)*(r_max - r_min) 120 else 121 r_j = r_min + rand()*(r_max - r_min) # between r_min and r_max 122 end 123 #r_j = r_min + rand()*(r_i - r_min) # between r_min and r_i 124 #R = rand() * (r_i + r_j) * max_padding_factor + 2. * (r_i + r_j) 125 R = r_i + r_j + padding 126 T = rand() * 2. * pi 127 return [x_i[1] + R * sin(T), x_i[2] + R * cos(T), x_i[3]], r_j 128 end 129 130 function generateRandomDirection() 131 return rand() * 2. * pi 132 end 133 134 function getPositionDistancedFromPoint(T::Real, x_i::Vector, dist::Real) 135 return [x_i[1] + dist * sin(T), x_i[2] + dist * cos(T), x_i[3]] 136 end 137 138 export irregularPacking! 139 """ 140 irregularPacking!(simulation[, radius_max, radius_min, sample_limit, 141 padding_factor, binary_radius_search, 142 binary_sampling_quality, thickness, seed, 143 plot_during_packing, verbose) 144 145 Generate a dense disc packing in 2D using Poisson disc sampling with O(N) 146 complexity, as described by [Robert Bridson (2007) "Fast Poisson disk sampling 147 in arbitrary dimensions"](https://doi.org/10.1145/1278780.1278807). The 148 `simulation` can be empty or already contain grains. However, an 149 `simulation.ocean` or `simulation.atmosphere` grid is required. 150 151 # Arguments 152 * `simulation::Simulation`: simulation object where grains are inserted. 153 * `radius_max::Real`: largest grain radius to use. 154 * `radius_min::Real`: smallest grain radius to use. 155 * `sample_limit::Integer=30`: number of points to sample around each grain 156 before giving up. 157 * `padding_factor::Real=0.`: if positive and `binary_radius_search = false`, try to 158 add an occasional grain from the current active grain 159 (`radius_max*padding_factor`). 160 * `binary_radius_search::Bool=false`: use a binary radius-sampling procedure to 161 fit the largest possible grains into the packing. This option will create 162 the highest packing density. 163 * `binary_sampling_quality::Real=100.`: the quality to enforce during the binary 164 radius search when `binary_radius_search = true`. Larger values create 165 denser packings but take longer to complete. 166 * `seed::Integer`: seed value to the pseudo-random number generator. 167 * `plot_during_packing::Bool=false`: produce successive plots as the packing is 168 generated. Requires gnuplot (default). 169 * `verbose::Bool=true`: show diagnostic information to stdout. 170 """ 171 function irregularPacking!(simulation::Simulation; 172 radius_max::Real=.1, 173 radius_min::Real=.1, 174 sample_limit::Integer=100, 175 padding_factor::Real=0., 176 binary_radius_search::Bool=false, 177 binary_sampling_quality::Real=100., 178 thickness::Real=1., 179 seed::Integer=1, 180 plot_during_packing::Bool=false, 181 verbose::Bool=true) 182 Random.seed!(seed) 183 184 active_list = Int[] # list of points to originate search from 185 i = 0 186 187 # Step 0: Use existing `grid` (ocean or atmosphere) for contact search 188 if typeof(simulation.ocean.input_file) != Bool 189 grid = simulation.ocean 190 elseif typeof(simulation.atmosphere.input_file) != Bool 191 grid = simulation.atmosphere 192 else 193 error("irregularPacking requires an ocean or atmosphere grid") 194 end 195 sortGrainsInGrid!(simulation, grid) 196 # save grid boundaries 197 sw, se, ne, nw = getGridCornerCoordinates(grid.xq, grid.yq) 198 width_x = se[1] - sw[1] # assume regular grid 199 width_y = nw[2] - sw[2] # assume regular grid 200 201 # Step 1: If grid is empty: select random initial sample and save its index 202 # to the background grid. Otherwise mark all existing grains as active 203 np_init = length(simulation.grains) 204 if isempty(simulation.grains) 205 r = rand()*(radius_max - radius_min) + radius_min 206 x0 = rand(2).*[width_x, width_y] + sw 207 addGrainCylindrical!(simulation, x0, r, thickness, color=1, 208 verbose=false) 209 sortGrainsInGrid!(simulation, grid) 210 push!(active_list, 1) 211 else 212 for idx=1:length(simulation.grains) 213 simulation.grains[idx].color = 1 214 push!(active_list, idx) 215 end 216 end 217 218 # Step 2: While the active list is not empty, choose a random index `i` from 219 # it. Generate up to `sample_limit` points chosen uniformly from the 220 # distance `(r_i + r_j)` around `x_i`. 221 # For each point in turn, check if it is within distance r of existing 222 # samples (using the background grid to only test nearby samples). If a 223 # point is adequately far from existing samples, emit it as the next sample 224 # and add it to the active list. If after `sample_limit` attempts no such 225 # point is found, instead remove `i` from the active list. 226 j = 0 227 x_active = zeros(3) 228 x_candidate = zeros(3) 229 r_active = 0.0 230 r_candidate = 0.0 231 T = 0.0 232 n = 0 233 neighbor_found = false 234 i_last_active = 0 235 if verbose 236 println("") 237 end 238 239 while !isempty(active_list) 240 241 # Draw a random grain from the list of active grains 242 i = active_list[rand(1:length(active_list))] 243 i_last_active = i 244 245 x_active = simulation.grains[i].lin_pos 246 r_active = simulation.grains[i].contact_radius 247 248 # Did the algoritm find a neighbor to the current active grain `i`? 249 neighbor_found = false 250 251 for j = 1:sample_limit 252 253 # Generate a candidate point 254 if binary_radius_search 255 # Generate a point positioned at r_active + radius_max from the 256 # position x_active. 257 T = generateRandomDirection() 258 r_candidate = radius_max 259 x_candidate = getPositionDistancedFromPoint(T, x_active, 260 r_active + r_candidate) 261 else 262 if j <= 2 # generate large grains during the first two samples 263 x_candidate, r_candidate = generateNeighboringPoint( 264 x_active, 265 r_active, 266 radius_max, 267 radius_min, 268 padding=padding_factor*radius_max) 269 else 270 x_candidate, r_candidate = generateNeighboringPoint( 271 x_active, 272 r_active, 273 radius_max, 274 radius_min) 275 end 276 end 277 278 # Make sure that the point is within the grid limits 279 if !(isPointInGrid(grid, x_candidate)) 280 continue # skip this candidate 281 end 282 283 # If the binary_radius_search is selected, try to adjust the radius 284 # to a value as large as possible 285 sortGrainsInGrid!(simulation, grid) 286 if binary_radius_search 287 288 # first test the maximum radius. If unsuccessful, iteratively 289 # find the optimal radius using binary searches 290 if !checkForContacts(simulation, grid, x_candidate, r_candidate, 291 return_when_overlap_found=true) 292 293 # 1. Set L to min and R to max of sampling range 294 r_L = radius_min 295 r_R = radius_max 296 297 # size of radius sampling step 298 dr = (r_R - r_L)/binary_sampling_quality 299 300 # 2. If L > R, the search terminates as unsuccessful 301 while r_L < r_R 302 303 # 3. Set r to the middle of the current range 304 r_candidate = (r_L + r_R)/2.0 305 x_candidate = getPositionDistancedFromPoint(T, x_active, 306 r_active + r_candidate) 307 #println("[$r_L, \t $r_candidate, \t $r_R]") 308 309 # 4. If r < target, set L to r+dr and go to step 2 310 if checkForContacts(simulation, grid, x_candidate, 311 r_candidate) <= 1 312 r_L = r_candidate + dr 313 314 else # 5. If r > target, set R to r-dr and go to step 2 315 r_R = r_candidate - dr 316 end 317 end 318 end 319 320 end 321 322 # if the grain candidate doesn't overlap with any other grains, 323 # add it and mark it as active 324 sortGrainsInGrid!(simulation, grid) 325 if checkForContacts(simulation, grid, x_candidate, r_candidate, 326 return_when_overlap_found=true) 327 #println("Added grain from parent $i") 328 addGrainCylindrical!(simulation, x_candidate, r_candidate, 329 thickness, color=1, verbose=false) 330 sortGrainsInGrid!(simulation, grid) 331 push!(active_list, length(simulation.grains)) 332 simulation.grains[i].color = 1 333 break 334 end 335 336 if j == sample_limit 337 # If no neighbors were found, delete the grain `i` from the list 338 # of active grains 339 simulation.grains[i].color = 0 340 filter!(f->f≠i, active_list) 341 end 342 end 343 if verbose 344 print("\rActive points: $(length(active_list)) ") 345 #println(active_list) 346 end 347 348 if plot_during_packing 349 n += 1 350 color = simulation.grains[i_last_active].color 351 simulation.grains[i_last_active].color = 2 352 filepostfix = @sprintf("packing.%05d.png", n) 353 plotGrains(simulation, filetype=filepostfix, show_figure=false, 354 palette_scalar="color", cbrange=[0.,2.]) 355 simulation.grains[i_last_active].color = color 356 end 357 358 end # end while !isempty(active_list) 359 360 if verbose 361 println("") 362 @info "Generated $(length(simulation.grains) - np_init) points" 363 end 364 end 365 366 export rasterPacking! 367 function rasterPacking!(sim::Simulation, 368 r_min::Real, 369 r_max::Real; 370 padding_factor::Real=0.1, 371 size_distribution::String="powerlaw", 372 size_distribution_parameter::Real=-1.8, 373 seed::Integer=1, 374 verbose::Bool=true) 375 376 r_rand = 0. 377 h = .5 # disc tickness 378 dx = r_max * 2. * (1. + padding_factor) # cell size 379 dx_padding = r_max * 2. * padding_factor 380 Random.seed!(seed) 381 382 np_init = length(sim.grains) 383 384 # Generate a grid spanning the entire domain, with cell width corresponding 385 # to the largest grain to be inserted 386 occupied = rasterMap(sim, dx) 387 388 # Add grains in unoccupied places 389 pos = zeros(2) 390 for ix=1:size(occupied, 1) 391 for iy=1:size(occupied, 2) 392 393 if occupied[ix,iy] 394 continue 395 end 396 397 if size_distribution == "powerlaw" 398 r_rand = Granular.randpower(1, size_distribution_parameter, 399 r_min, r_max) 400 elseif size_distribution == "uniform" 401 r_rand = rand()*(r_max - r_min) + r_min 402 end 403 404 # Determine position from grid index and sample randomly from within 405 # padding 406 pos = [ix*dx - .5*dx, iy*dx - .5*dx] .+ 407 rand(2) .* dx_padding .- .5*dx_padding 408 409 addGrainCylindrical!(sim, pos, r_rand, h, verbose=false) 410 411 end 412 end 413 if verbose 414 @info "Generated $(length(sim.grains) - np_init) points" 415 end 416 end 417 418 """ 419 rasterMap(sim, dx) 420 421 Returns a rasterized map of grain extent in the domain with length `L` and a 422 pixel size of `dx`. The function will return a map of `Bool` type with size 423 `floor.(L./dx)`. 424 425 * Arguments 426 - `sim::Simulation`: simulation object containing the grains. 427 - `dx::Real`: pixel size to use for the raster map. 428 429 """ 430 function rasterMap(sim::Simulation, dx::Real) 431 432 # Use existing `grid` (ocean or atmosphere) for contact search 433 if typeof(sim.ocean.input_file) != Bool 434 grid = sim.ocean 435 elseif typeof(sim.atmosphere.input_file) != Bool 436 grid = sim.atmosphere 437 else 438 error("rasterMap(...) requires an ocean or atmosphere grid") 439 end 440 # save grid boundaries 441 if grid.regular_grid 442 L = grid.L[1:2] 443 origo = grid.origo 444 else 445 sw, se, ne, nw = getGridCornerCoordinates(grid.xq, grid.yq) 446 L = [se[1] - sw[1], nw[2] - sw[2]] 447 origo = [sw[1], sw[2]] 448 end 449 dims = floor.(L./dx) 450 occupied = zeros(Bool, convert(Dims, (dims[1], dims[2]))) 451 452 # Loop over existing grains and mark their extent in the `occupied` array 453 i = 0; j = 0 454 min_i = 0; min_j = 0 455 max_i = 0; max_j = 0 456 cell_mid_point = zeros(2) 457 dist = sqrt(2.0*(dx/2.0)^2.) 458 for grain in sim.grains 459 460 # Find center position in `occupied` grid 461 #i, j = Int.(floor.((grain.lin_pos - origo) ./ dx)) + [1,1] 462 463 # Find corner indexes for box spanning the grain 464 min_i, min_j = Int.(floor.((grain.lin_pos[1:2] - origo .- 465 grain.contact_radius) ./ dx)) .+ [1,1] 466 max_i, max_j = Int.(floor.((grain.lin_pos[1:2] - origo .+ 467 grain.contact_radius) ./ dx)) .+ [1,1] 468 469 # For each cell in box, check if the grain is contained 470 for i = min_i:max_i 471 for j = min_j:max_j 472 cell_mid_point = dx .* Vector{Float64}([i,j]) .- 0.5 * dx 473 474 if (norm(cell_mid_point - grain.lin_pos[1:2]) - 475 grain.contact_radius < dist) 476 occupied[i,j] = true 477 end 478 end 479 end 480 end 481 return occupied 482 end