tcontact_search.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
---
tcontact_search.jl (12084B)
---
1 using LinearAlgebra
2
3 ## Contact mapping
4 export findContacts!
5 """
6 findContacts!(simulation[, method])
7
8 Top-level function to perform an inter-grain contact search, based on grain
9 linear positions and contact radii.
10
11 The simplest contact search algorithm (`method="all to all"`) is the most
12 computationally expensive (O(n^2)). The method "ocean grid" bins the grains
13 into their corresponding cells on the ocean grid and searches for contacts only
14 within the vicinity. When this method is applied, it is assumed that the
15 `contact_radius` values of the grains are *smaller than half the cell size*.
16
17 # Arguments
18 * `simulation::Simulation`: the simulation object containing the grains.
19 * `method::String`: the contact-search method to apply. Valid options are "all
20 to all" and "ocean grid".
21 """
22 function findContacts!(simulation::Simulation;
23 method::String = "all to all")
24
25 if method == "all to all"
26 findContactsAllToAll!(simulation)
27
28 elseif method == "ocean grid"
29 findContactsInGrid!(simulation, simulation.ocean)
30
31 elseif method == "atmosphere grid"
32 findContactsInGrid!(simulation, simulation.atmosphere)
33
34 else
35 error("Unknown contact search method '$method'")
36 end
37 nothing
38 end
39
40 export interGrainPositionVector
41 """
42 interGrainPositionVector(simulation, i, j)
43
44 Returns a `vector` pointing from grain `i` to grain `j` in the
45 `simulation`.
46
47 # Arguments
48 * `simulation::Simulation`: the simulation object containing the grains.
49 * `i::Int`: index of the first grain.
50 * `j::Int`: index of the second grain.
51 """
52 function interGrainPositionVector(simulation::Simulation,
53 i::Int, j::Int)
54 @inbounds return simulation.grains[i].lin_pos -
55 simulation.grains[j].lin_pos
56 end
57
58 """
59 position_ij is the inter-grain position vector, and can be found with
60 interGrainPositionVector().
61 """
62 function findOverlap(simulation::Simulation, i::Int, j::Int,
63 position_ij::Vector{Float64})
64 @inbounds return norm(position_ij) - (simulation.grains[i].contact_radius
65 + simulation.grains[j].contact_radius)
66 end
67
68 export findContactsAllToAll!
69 """
70 findContactsAllToAll!(simulation)
71
72 Perform an O(n^2) all-to-all contact search between all grains in the
73 `simulation` object.
74 """
75 function findContactsAllToAll!(simulation::Simulation)
76
77 if simulation.ocean.bc_west > 1 ||
78 simulation.ocean.bc_east > 1 ||
79 simulation.ocean.bc_north > 1 ||
80 simulation.ocean.bc_south > 1
81 error("Ocean boundary conditions do not work with all-to-all " *
82 "contact search")
83 end
84 if simulation.atmosphere.bc_west > 1 ||
85 simulation.atmosphere.bc_east > 1 ||
86 simulation.atmosphere.bc_north > 1 ||
87 simulation.atmosphere.bc_south > 1
88 error("Atmopshere boundary conditions do not work with " *
89 "all-to-all contact search")
90 end
91
92 @inbounds for i = 1:length(simulation.grains)
93
94 # Check contacts with other grains
95 for j = 1:length(simulation.grains)
96 checkAndAddContact!(simulation, i, j)
97 end
98 end
99 nothing
100 end
101
102 export findContactsInGrid!
103 """
104 findContactsInGrid!(simulation)
105
106 Perform an O(n*log(n)) cell-based contact search between all grains in the
107 `simulation` object.
108 """
109 function findContactsInGrid!(simulation::Simulation, grid::Any)
110
111 distance_modifier = [0., 0., 0.]
112 i_corrected = 0
113 j_corrected = 0
114
115 for idx_i = 1:length(simulation.grains)
116
117 if typeof(grid) == Ocean
118 grid_pos = simulation.grains[idx_i].ocean_grid_pos
119 elseif typeof(grid) == Atmosphere
120 grid_pos = simulation.grains[idx_i].atmosphere_grid_pos
121 else
122 error("Grid type not understood")
123 end
124 nx, ny = size(grid.xh)
125
126 for i = (grid_pos[1] - 1):(grid_pos[1] + 1)
127 for j = (grid_pos[2] - 1):(grid_pos[2] + 1)
128
129 # correct indexes if necessary
130 i_corrected, j_corrected, distance_modifier =
131 periodicBoundaryCorrection!(grid, i, j)
132
133 # skip iteration if target still falls outside grid after
134 # periodicity correction
135 if i_corrected < 1 || i_corrected > nx ||
136 j_corrected < 1 || j_corrected > ny
137 continue
138 end
139
140 @inbounds for idx_j in grid.grain_list[i_corrected, j_corrected]
141 checkAndAddContact!(simulation, idx_i, idx_j,
142 distance_modifier)
143 end
144 end
145 end
146 end
147 nothing
148 end
149
150
151 export checkForContacts
152 """
153 checkForContacts(grid, position, radius)
154
155 Perform an O(n*log(n)) cell-based contact search between a candidate grain with
156 position `position` and `radius`, against all grains registered in the `grid`.
157 Returns the number of contacts that were found as an `Integer` value, unless
158 `return_when_overlap_found` is `true`.
159
160 # Arguments
161 * `simulation::Simulation`: Simulation object containing grain positions.
162 * `grid::Any`: `Ocean` or `Atmosphere` grid containing sorted particles.
163 * `x_candidate::Vector{Float64}`: Candidate center position to probe for
164 contacts with existing grains [m].
165 * `r_candidate::Float64`: Candidate radius [m].
166 * `return_when_overlap_found::Bool` (default: `false`): Return `true` if no
167 contacts are found, or return `false` as soon as a contact is found.
168 """
169 function checkForContacts(simulation::Simulation,
170 grid::Any,
171 x_candidate::Vector{Float64},
172 r_candidate::Float64;
173 return_when_overlap_found::Bool=false)
174
175 distance_modifier = zeros(3)
176 overlaps_found::Integer = 0
177
178 # Inter-grain position vector and grain overlap
179 ix, iy = findCellContainingPoint(grid, x_candidate)
180
181 # Check for overlap with existing grains
182 for ix_=(ix - 1):(ix + 1)
183 for iy_=(iy - 1):(iy + 1)
184
185 # correct indexes if necessary
186 ix_corrected, iy_corrected, distance_modifier =
187 periodicBoundaryCorrection!(grid, ix_, iy_)
188
189 # skip iteration if target still falls outside grid after
190 # periodicity correction
191 if ix_corrected < 1 || ix_corrected > size(grid.xh)[1] ||
192 iy_corrected < 1 || iy_corrected > size(grid.xh)[2]
193 continue
194 end
195
196 @inbounds for idx in grid.grain_list[ix_corrected, iy_corrected]
197 if norm(x_candidate - simulation.grains[idx].lin_pos +
198 distance_modifier) -
199 (simulation.grains[idx].contact_radius + r_candidate) < 0.0
200
201 if return_when_overlap_found
202 return false
203 end
204 overlaps_found += 1
205 end
206 end
207 end
208 end
209 if return_when_overlap_found
210 return true
211 else
212 return overlaps_found
213 end
214 end
215
216 """
217 periodicBoundaryCorrection!(grid::Any, i::Integer, j::Integer,
218 i_corrected::Integer, j_corrected::Integer)
219
220 Determine the geometric correction and grid-index adjustment required across
221 periodic boundaries.
222 """
223 function periodicBoundaryCorrection!(grid::Any, i::Integer, j::Integer)
224
225 # vector for correcting inter-particle distance in case of
226 # boundary periodicity
227 distance_modifier = zeros(3)
228
229 # i and j are not corrected for periodic boundaries
230 i_corrected = i
231 j_corrected = j
232
233 # only check for contacts within grid boundaries, and wrap
234 # around if they are periodic
235 if i < 1 || i > size(grid.xh)[1] || j < 1 || j > size(grid.xh)[2]
236
237 if i < 1 && grid.bc_west == 2 # periodic -x
238 distance_modifier[1] = grid.xq[end] - grid.xq[1]
239 i_corrected = size(grid.xh)[1]
240 elseif i > size(grid.xh)[1] && grid.bc_east == 2 # periodic +x
241 distance_modifier[1] = -(grid.xq[end] - grid.xq[1])
242 i_corrected = 1
243 end
244
245 if j < 1 && grid.bc_south == 2 # periodic -y
246 distance_modifier[2] = grid.yq[end] - grid.yq[1]
247 j_corrected = size(grid.xh)[2]
248 elseif j > size(grid.xh)[2] && grid.bc_north == 2 # periodic +y
249 distance_modifier[2] = -(grid.yq[end] - grid.yq[1])
250 j_corrected = 1
251 end
252 end
253
254 return i_corrected, j_corrected, distance_modifier
255 end
256
257 export checkAndAddContact!
258 """
259 checkAndAddContact!(simulation, i, j)
260
261 Check for contact between two grains and register the interaction in the
262 `simulation` object. The indexes of the two grains is stored in
263 `simulation.contact_pairs` as `[i, j]`. The overlap vector is parallel to a
264 straight line connecting the grain centers, points away from grain `i` and
265 towards `j`, and is stored in `simulation.overlaps`. A zero-length vector is
266 written to `simulation.contact_parallel_displacement`.
267
268 # Arguments
269 * `simulation::Simulation`: the simulation object containing the grains.
270 * `i::Int`: index of the first grain.
271 * `j::Int`: index of the second grain.
272 * `distance_Modifier::Vector{Float64}`: vector modifying percieved
273 inter-particle distance, which is used for contact search across periodic
274 boundaries.
275 """
276 function checkAndAddContact!(sim::Simulation, i::Int, j::Int,
277 distance_modifier::Vector{Float64} = [0., 0., 0.])
278 if i < j
279
280 @inbounds if !sim.grains[i].enabled || !sim.grains[j].enabled
281 return
282 end
283
284 # Inter-grain position vector and grain overlap
285 position_ij = interGrainPositionVector(sim, i, j) + distance_modifier
286 overlap_ij = findOverlap(sim, i, j, position_ij)
287
288 contact_found = false
289
290 # Check if contact is already registered, and update position if so
291 for ic=1:sim.Nc_max
292 @inbounds if sim.grains[i].contacts[ic] == j
293 contact_found = true
294 @inbounds sim.grains[i].position_vector[ic] = position_ij
295 nothing # contact already registered
296 end
297 end
298
299 # Check if grains overlap (overlap when negative)
300 if overlap_ij < 0.
301
302 # Register as new contact in first empty position
303 if !contact_found
304
305 for ic=1:(sim.Nc_max + 1)
306
307 # Test if this contact exceeds the number of contacts
308 if ic == (sim.Nc_max + 1)
309 for ic=1:sim.Nc_max
310 @warn "grains[$i].contacts[$ic] = " *
311 "$(sim.grains[i].contacts[ic])"
312 @warn "grains[$i].contact_age[$ic] = " *
313 "$(sim.grains[i].contact_age[ic])"
314 end
315 error("contact $i-$j exceeds max. number of " *
316 "contacts (sim.Nc_max = $(sim.Nc_max)) " *
317 "for grain $i")
318 end
319
320 # Register as new contact
321 @inbounds if sim.grains[i].contacts[ic] == 0 # empty
322 @inbounds sim.grains[i].n_contacts += 1
323 @inbounds sim.grains[j].n_contacts += 1
324 @inbounds sim.grains[i].contacts[ic] = j
325 @inbounds sim.grains[i].position_vector[ic] =
326 position_ij
327 @inbounds fill!(sim.grains[i].
328 contact_parallel_displacement[ic] , 0.)
329 @inbounds fill!(sim.grains[i].contact_rotation[ic] , 0.)
330 @inbounds sim.grains[i].contact_age[ic] = 0.
331 @inbounds sim.grains[i].contact_area[ic] = 0.
332 @inbounds sim.grains[i].compressive_failure[ic] = false
333 break
334 end
335 end
336 end
337 end
338 end
339 nothing
340 end