Granular.jl

commit e647a43e390dad3f690891f5f9cfffa924f092f5
parent b0918c98604e04ea2c4f983c73e589cf7b7e9646
Author: Anders Damsgaard <andersd@riseup.net>
Date:   Wed,  1 Nov 2017 22:16:48 -0400

store inter-particle position vector on contact

Diffstat:
M
src/contact_search.jl
 | 
18
+++++++++---------
M
src/datatypes.jl
 | 
1
+
M
src/grain.jl
 | 
10
+++++++---
M
src/interaction.jl
 | 
17
++++-------------

4 files changed, 21 insertions(+), 25 deletions(-)
diff --git a/src/contact_search.jl b/src/contact_search.jl
@@ -189,18 +189,18 @@ function checkAndAddContact!(sim::Simulation, i::Int, j::Int,
 
         contact_found = false
 
+        # Check if contact is already registered, and update position if so
+        for ic=1:sim.Nc_max
+            @inbounds if sim.grains[i].contacts[ic] == j
+                contact_found = true
+                @inbounds sim.grains[i].position_vector[ic] .= position_ij
+                nothing  # contact already registered
+            end
+        end
+
         # Check if grains overlap (overlap when negative)
         if overlap_ij < 0.
 
-            # Check if contact is already registered, and update position if so
-            for ic=1:sim.Nc_max
-                @inbounds if sim.grains[i].contacts[ic] == j
-                    contact_found = true
-                    @inbounds sim.grains[i].position_vector[ic] .= position_ij
-                    break  # contact already registered
-                end
-            end
-
             # Register as new contact in first empty position
             if !contact_found
 
diff --git a/src/datatypes.jl b/src/datatypes.jl
@@ -58,6 +58,7 @@ mutable struct GrainCylindrical
     ocean_grid_pos::Vector{Int}
     atmosphere_grid_pos::Vector{Int}
     contacts::Vector{Int}
+    position_vector::Vector{Vector{Float64}}
     contact_parallel_displacement::Vector{Vector{Float64}}
     contact_age::Vector{Float64}
 
diff --git a/src/grain.jl b/src/grain.jl
@@ -195,10 +195,11 @@ function addGrainCylindrical!(simulation::Simulation,
     end
 
     contacts::Array{Int, 1} = zeros(Int, simulation.Nc_max)
-    contact_parallel_displacement =
-        Vector{Vector{Float64}}(simulation.Nc_max)
-        contact_age::Vector{Float64} = zeros(Float64, simulation.Nc_max)
+    position_vector = Vector{Vector{Float64}}(simulation.Nc_max)
+    contact_parallel_displacement = Vector{Vector{Float64}}(simulation.Nc_max)
+    contact_age::Vector{Float64} = zeros(Float64, simulation.Nc_max)
     for i=1:simulation.Nc_max
+        position_vector[i] = zeros(2)
         contact_parallel_displacement[i] = zeros(2)
     end
 
@@ -252,6 +253,7 @@ function addGrainCylindrical!(simulation::Simulation,
                                  ocean_grid_pos,
                                  atmosphere_grid_pos,
                                  contacts,
+                                 position_vector,
                                  contact_parallel_displacement,
                                  contact_age,
 
@@ -548,6 +550,7 @@ function printGrainInfo(f::GrainCylindrical)
     println("  pressure:   $(f.pressure) Pa")
     println("  n_contacts: $(f.n_contacts)")
     println("  contacts:   $(f.contacts)")
+    println("  pos_ij:     $(f.position_vector)\n")
     println("  delta_t:    $(f.contact_parallel_displacement)\n")
 
     println("  granular_stress:   $(f.granular_stress) Pa")
@@ -655,6 +658,7 @@ function compareGrains(if1::GrainCylindrical, if2::GrainCylindrical)
     Base.Test.@test if1.n_contacts == if2.n_contacts
     Base.Test.@test if1.ocean_grid_pos == if2.ocean_grid_pos
     Base.Test.@test if1.contacts == if2.contacts
+    Base.Test.@test if1.position_vector == if2.position_vector
     Base.Test.@test if1.contact_parallel_displacement == 
         if2.contact_parallel_displacement
     Base.Test.@test if1.contact_age ≈ if2.contact_age
diff --git a/src/interaction.jl b/src/interaction.jl
@@ -16,17 +16,7 @@ function interact!(simulation::Simulation)
                 continue
             end
 
-            #=if norm(simulation.grains[i].lin_pos - 
-                    simulation.grains[j].lin_pos) - 
-                    (simulation.grains[i].contact_radius + 
-                    simulation.grains[j].contact_radius) > 0.
-
-                simulation.grains[i].contacts[ic] = 0  # remove contact
-                simulation.grains[i].n_contacts -= 1
-                simulation.grains[j].n_contacts -= 1
-            else=#
             interactGrains!(simulation, i, j, ic)
-            #end
         end
     end
 
@@ -58,8 +48,9 @@ function interactGrains!(simulation::Simulation, i::Int, j::Int, ic::Int)
     force_n = 0.  # Contact-normal force
     force_t = 0.  # Contact-parallel (tangential) force
 
-    # Inter-position vector
-    p = simulation.grains[i].lin_pos - simulation.grains[j].lin_pos
+    # Inter-position vector, use stored value which is corrected for boundary
+    # periodicity
+    p = simulation.grains[i].position_vector[ic]
     dist = norm(p)
 
     r_i = simulation.grains[i].contact_radius
@@ -81,7 +72,7 @@ function interactGrains!(simulation::Simulation, i::Int, j::Int, ic::Int)
                                 simulation.grains[j].ang_vel)
 
     # Correct old tangential displacement for contact rotation and add new
-    δ_t0 =simulation.grains[i].contact_parallel_displacement[ic]
+    δ_t0 = simulation.grains[i].contact_parallel_displacement[ic]
     δ_t = dot(t, δ_t0 - (n*dot(n, δ_t0))) + vel_t*simulation.time_step
 
     # Effective radius