tadd elastic interaction with walls - Granular.jl - Julia package for granular dynamics simulation
 (HTM) git clone git://src.adamsgaard.dk/Granular.jl
 (DIR) Log
 (DIR) Files
 (DIR) Refs
 (DIR) README
 (DIR) LICENSE
       ---
 (DIR) commit c480bcc2d14ba436fa4663aa46cd6d503c8ca905
 (DIR) parent 6e8adb0d8a14aa77aaa516f197c72bbe242a9fbb
 (HTM) Author: Anders Damsgaard <andersd@riseup.net>
       Date:   Wed, 15 Nov 2017 11:53:04 -0500
       
       add elastic interaction with walls
       
       Diffstat:
         M examples/image.jl                   |       6 +++---
         M examples/two-grains.jl              |       1 +
         M src/interaction.jl                  |      28 +++++++++++++++++-----------
         M src/simulation.jl                   |       1 +
         M src/wall.jl                         |       8 ++++++--
         M test/wall.jl                        |      72 ++++++++++++++++++++++++++++---
       
       6 files changed, 94 insertions(+), 22 deletions(-)
       ---
 (DIR) diff --git a/examples/image.jl b/examples/image.jl
       t@@ -6,7 +6,7 @@ import Colors
        
        const verbose = true
        
       -const img_file = "Cundall2008.png"
       +const img_file = "aadcroft.png"
        
        img = FileIO.load(img_file)
        
       t@@ -20,10 +20,10 @@ end
        
        const img_bw = Colors.Gray.(img)
        
       -#const forcing = "gyres"
       +const forcing = "gyres"
        #const forcing = "down"
        #const forcing = "convergent"
       -const forcing = "sandpile"
       +#const forcing = "sandpile"
        
        const dx = 1.
        const dy = dx
 (DIR) diff --git a/examples/two-grains.jl b/examples/two-grains.jl
       t@@ -39,3 +39,4 @@ info("Kinetic energy before: $E_kin_before J")
        info("Kinetic energy after:  $E_kin_after J")
        
        
       +Granular.render(sim, animation=true)
 (DIR) diff --git a/src/interaction.jl b/src/interaction.jl
       t@@ -54,18 +54,24 @@ function interactWalls!(sim::Simulation)
        
                    # get overlap distance by projecting grain position onto wall-normal
                    # vector
       -            δ_n = dot(sim.walls[iw].normal, sim.grains[i].lin_pos) -
       -                sim.walls[iw].pos - sim.grains[i].contact_radius
       -
       -            if sim.grains[i].youngs_modulus > 0.
       -                k_n = sim.grains[i].youngs_modulus *
       -                    sim.grains[i].thickness
       -            else
       -                k_n = sim.grains[i].contact_stiffness_normal
       +            δ_n = abs(dot(sim.walls[iw].normal, sim.grains[i].lin_pos) -
       +                      sim.walls[iw].pos) - sim.grains[i].contact_radius
       +
       +            #println(dot(sim.walls[iw].normal, sim.grains[i].lin_pos))
       +            #println(sim.walls[iw].pos)
       +            #println(sim.grains[i].contact_radius)
       +            #println(δ_n)
       +
       +            if δ_n < 0.
       +                if sim.grains[i].youngs_modulus > 0.
       +                    k_n = sim.grains[i].youngs_modulus * sim.grains[i].thickness
       +                else
       +                    k_n = sim.grains[i].contact_stiffness_normal
       +                end
       +
       +                sim.walls[iw].force += -k_n * δ_n
       +                sim.grains[i].force += k_n * abs(δ_n) * sim.walls[iw].normal
                    end
       -
       -            sim.walls[iw].force += k_n * abs(δ_n)
       -            sim.grains[i].force += k_n * abs(δ_n) * sim.grains[iw].normal
                end
            end
        end
 (DIR) diff --git a/src/simulation.jl b/src/simulation.jl
       t@@ -161,6 +161,7 @@ function run!(simulation::Simulation;
                end
        
                interact!(simulation)
       +        interactWalls!(simulation)
        
                if typeof(simulation.ocean.input_file) != Bool
                    addOceanDrag!(simulation)
 (DIR) diff --git a/src/wall.jl b/src/wall.jl
       t@@ -98,7 +98,9 @@ function addWallLinearFrictionless!(simulation::Simulation,
                for grain in simulation.grains
                    mass += grain.mass
                end
       -        info("Setting wall mass to total grain mass: $mass kg")
       +        if verbose
       +            info("Setting wall mass to total grain mass: $mass kg")
       +        end
            end
        
            # if not set, set wall thickness to equal largest grain thickness
       t@@ -113,7 +115,9 @@ function addWallLinearFrictionless!(simulation::Simulation,
                        thickness = grain.thickness
                    end
                end
       -        info("Setting wall thickness to largest grain thickness: $thickness m")
       +        if verbose
       +            info("Setting wall thickness to max grain thickness: $thickness m")
       +        end
            end
        
            # Create wall object
 (DIR) diff --git a/test/wall.jl b/test/wall.jl
       t@@ -5,7 +5,7 @@
        info("#### $(basename(@__FILE__)) ####")
        
        info("Testing argument value checks")
       -sim = Granular.createSimulation(id="test")
       +sim = Granular.createSimulation()
        Granular.addGrainCylindrical!(sim, [ 0., 0.], 10., 2., verbose=false)
        @test_throws ErrorException Granular.addWallLinearFrictionless!(sim,
                                                                        [.1, .1, .1],
       t@@ -16,13 +16,13 @@ Granular.addGrainCylindrical!(sim, [ 0., 0.], 10., 2., verbose=false)
        @test_throws ErrorException Granular.addWallLinearFrictionless!(sim,
                                                                        [.1, .1, .1],
                                                                        1.)
       -sim = Granular.createSimulation(id="test")
       +sim = Granular.createSimulation()
        @test_throws ErrorException Granular.addWallLinearFrictionless!(sim, [1., 0.],
                                                                        1.)
        
        
        info("Check that wall mass equals total grain mass and max. thickness")
       -sim = Granular.createSimulation(id="test")
       +sim = Granular.createSimulation()
        @test length(sim.walls) == 0
        Granular.addGrainCylindrical!(sim, [ 0., 0.], 10., 2., verbose=false)
        sim.grains[1].mass = 1.0
       t@@ -32,12 +32,12 @@ Granular.addWallLinearFrictionless!(sim, [1., 0.], 1., verbose=true)
        @test sim.walls[1].thickness ≈ 2.0
        
        info("Test wall surface area and defined normal stress")
       -sim = Granular.createSimulation(id="test")
       +sim = Granular.createSimulation()
        sim.ocean = Granular.createRegularOceanGrid([1, 1, 1], [10., 20., 1.0])
        Granular.addGrainCylindrical!(sim, [ 0., 0.], 10., 2., verbose=false)
        sim.grains[1].mass = 1.0
       -Granular.addWallLinearFrictionless!(sim, [1., 0.], 1., verbose=true)
       -Granular.addWallLinearFrictionless!(sim, [0., 1.], 1., verbose=true)
       +Granular.addWallLinearFrictionless!(sim, [1., 0.], 1., verbose=false)
       +Granular.addWallLinearFrictionless!(sim, [0., 1.], 1., verbose=false)
        @test length(sim.walls) == 2
        @test sim.walls[1].mass ≈ 1.0
        @test sim.walls[1].thickness ≈ 2.0
       t@@ -46,4 +46,64 @@ Granular.addWallLinearFrictionless!(sim, [0., 1.], 1., verbose=true)
        @test Granular.getWallSurfaceArea(sim, 1) ≈ 20.0*2.0
        @test Granular.getWallSurfaceArea(sim, 2) ≈ 10.0*2.0
        
       +info("Test wall-grain interaction")
       +
       +info("...wall present but no contact")
       +sim = Granular.createSimulation()
       +sim.ocean = Granular.createRegularOceanGrid([1, 1, 1], [10., 20., 1.0])
       +Granular.addGrainCylindrical!(sim, [ 0., 0.], 1., 2., verbose=false)
       +Granular.addWallLinearFrictionless!(sim, [1., 0.], -1.01, verbose=false)
       +Granular.interactWalls!(sim)
       +@test sim.walls[1].force ≈ 0.
       +@test sim.grains[1].force[1] ≈ 0.
       +@test sim.grains[1].force[2] ≈ 0.
       +
       +info("...wall present but no contact")
       +sim = Granular.createSimulation()
       +sim.ocean = Granular.createRegularOceanGrid([1, 1, 1], [10., 20., 1.0])
       +Granular.addGrainCylindrical!(sim, [ 0., 0.], 1., 2., verbose=false)
       +Granular.addWallLinearFrictionless!(sim, [-1., 0.], +2.01, verbose=false)
       +Granular.interactWalls!(sim)
       +@test sim.walls[1].force ≈ 0.
       +@test sim.grains[1].force[1] ≈ 0.
       +@test sim.grains[1].force[2] ≈ 0.
       +
       +info("...wall at -x")
       +sim = Granular.createSimulation()
       +sim.ocean = Granular.createRegularOceanGrid([1, 1, 1], [10., 20., 1.0])
       +Granular.addGrainCylindrical!(sim, [ 0., 0.], 1., 2., verbose=false)
       +Granular.addWallLinearFrictionless!(sim, [1., 0.], -1. + .01, verbose=false)
       +Granular.interactWalls!(sim)
       +@test sim.walls[1].force > 0.
       +@test sim.grains[1].force[1] > 0.
       +@test sim.grains[1].force[2] ≈ 0.
       +
       +info("...wall at +x")
       +sim = Granular.createSimulation()
       +sim.ocean = Granular.createRegularOceanGrid([1, 1, 1], [10., 20., 1.0])
       +Granular.addGrainCylindrical!(sim, [ 0., 0.], 1., 2., verbose=false)
       +Granular.addWallLinearFrictionless!(sim, [-1., 0.], +1. - .01, verbose=false)
       +Granular.interactWalls!(sim)
       +@test sim.walls[1].force > 0.
       +@test sim.grains[1].force[1] < 0.
       +@test sim.grains[1].force[2] ≈ 0.
        
       +info("...wall at -y")
       +sim = Granular.createSimulation()
       +sim.ocean = Granular.createRegularOceanGrid([1, 1, 1], [10., 20., 1.0])
       +Granular.addGrainCylindrical!(sim, [ 0., 0.], 1., 2., verbose=false)
       +Granular.addWallLinearFrictionless!(sim, [0., 1.], -1. + .01, verbose=false)
       +Granular.interactWalls!(sim)
       +@test sim.walls[1].force > 0.
       +@test sim.grains[1].force[1] ≈ 0.
       +@test sim.grains[1].force[2] > 0.
       +
       +info("...wall at +y")
       +sim = Granular.createSimulation()
       +sim.ocean = Granular.createRegularOceanGrid([1, 1, 1], [10., 20., 1.0])
       +Granular.addGrainCylindrical!(sim, [ 0., 0.], 1., 2., verbose=false)
       +Granular.addWallLinearFrictionless!(sim, [0., -1.], +1. - .01, verbose=false)
       +Granular.interactWalls!(sim)
       +@test sim.walls[1].force > 0.
       +@test sim.grains[1].force[1] ≈ 0.
       +@test sim.grains[1].force[2] < 0.