tpacking.jl - Granular.jl - Julia package for granular dynamics simulation
HTML git clone git://src.adamsgaard.dk/Granular.jl
DIR Log
DIR Files
DIR Refs
DIR README
DIR LICENSE
---
tpacking.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