contact_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