tadd ! to addIceFloeCylindrical name, store Nc_max in simulation object - 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 f51076d53050ade8a00bc1216257be807a5df907
 (DIR) parent 9bab04cfb39c08c4b92304968959fca205ec46e8
 (HTM) Author: Anders Damsgaard <andersd@riseup.net>
       Date:   Fri, 23 Jun 2017 10:20:00 -0400
       
       add ! to addIceFloeCylindrical name, store Nc_max in simulation object
       
       Diffstat:
         M examples/double_gyre.jl             |      18 +++++++++---------
         M examples/strait.jl                  |      39 ++++++++++++++++---------------
         M src/contact_search.jl               |       6 +++---
         M src/datatypes.jl                    |       2 ++
         M src/icefloe.jl                      |     104 +++++++++++++++----------------
         M src/interaction.jl                  |       2 +-
         M src/io.jl                           |       2 +-
         M src/simulation.jl                   |       8 ++++++--
         M test/atmosphere.jl                  |       2 +-
         M test/cohesion.jl                    |       8 ++++----
         M test/collision-2floes-normal.jl     |      16 ++++++++--------
         M test/collision-2floes-oblique.jl    |      36 ++++++++++++++++----------------
         M test/collision-5floes-normal.jl     |      40 ++++++++++++++++----------------
         M test/contact-search-and-geometry.jl |      44 ++++++++++++++++----------------
         M test/grid.jl                        |      60 ++++++++++++++++----------------
         M test/icefloe.jl                     |      26 +++++++++++++++-----------
         M test/jld.jl                         |       4 ++--
         M test/ocean.jl                       |       2 +-
         M test/profiling.jl                   |       4 ++--
         M test/temporal.jl                    |       2 +-
         M test/vtk.jl                         |       4 ++--
       
       21 files changed, 218 insertions(+), 211 deletions(-)
       ---
 (DIR) diff --git a/examples/double_gyre.jl b/examples/double_gyre.jl
       t@@ -35,18 +35,18 @@ h = 1.
        
        ## N-S wall segments
        for y in linspace(r, L[2]-r, Int(round((L[2] - 2.*r)/(r*2))))
       -    SeaIce.addIceFloeCylindrical(sim, [r, y], r, h, fixed=true,
       -                                 verbose=false)
       -    SeaIce.addIceFloeCylindrical(sim, [L[1]-r, y], r, h, fixed=true,
       -                                 verbose=false)
       +    SeaIce.addIceFloeCylindrical!(sim, [r, y], r, h, fixed=true,
       +                                  verbose=false)
       +    SeaIce.addIceFloeCylindrical!(sim, [L[1]-r, y], r, h, fixed=true,
       +                                  verbose=false)
        end
        
        ## E-W wall segments
        for x in linspace(3.*r, L[1]-3.*r, Int(round((L[1] - 6.*r)/(r*2))))
       -    SeaIce.addIceFloeCylindrical(sim, [x, r], r, h, fixed=true,
       -                                 verbose=false)
       -    SeaIce.addIceFloeCylindrical(sim, [x, L[2]-r], r, h, fixed=true,
       -                                 verbose=false)
       +    SeaIce.addIceFloeCylindrical!(sim, [x, r], r, h, fixed=true,
       +                                  verbose=false)
       +    SeaIce.addIceFloeCylindrical!(sim, [x, L[2]-r], r, h, fixed=true,
       +                                  verbose=false)
        end
        
        n_walls = length(sim.ice_floes)
       t@@ -70,7 +70,7 @@ for y in (4.*r + noise_amplitude):(2.*r + floe_padding):(L[2] - 4.*r -
                x_ = x + noise_amplitude*(0.5 - Base.Random.rand())
                y_ = y + noise_amplitude*(0.5 - Base.Random.rand())
        
       -        SeaIce.addIceFloeCylindrical(sim, [x_, y_], r, h, verbose=false)
       +        SeaIce.addIceFloeCylindrical!(sim, [x_, y_], r, h, verbose=false)
            end
        end
        n = length(sim.ice_floes) - n_walls
 (DIR) diff --git a/examples/strait.jl b/examples/strait.jl
       t@@ -25,15 +25,16 @@ r_walls = r_min
        for y in linspace((L[2] - Ly_constriction)/2.,
                          Ly_constriction + (L[2] - Ly_constriction)/2., 
                          Int(round(Ly_constriction/(r_walls*2))))
       -    SeaIce.addIceFloeCylindrical(sim, [(Lx - Lx_constriction)/2., y], r_walls, 
       -                                 h, fixed=true, verbose=false)
       +    SeaIce.addIceFloeCylindrical!(sim, [(Lx - Lx_constriction)/2., y], r_walls, 
       +                                  h, fixed=true, verbose=false)
        end
        for y in linspace((L[2] - Ly_constriction)/2.,
                          Ly_constriction + (L[2] - Ly_constriction)/2., 
                          Int(round(Ly_constriction/(r_walls*2))))
       -    SeaIce.addIceFloeCylindrical(sim,
       -                                 [Lx_constriction + (L[1] - Lx_constriction)/2., 
       -                                  y], r_walls, h, fixed=true, verbose=false)
       +    SeaIce.addIceFloeCylindrical!(sim,
       +                                  [Lx_constriction +
       +                                   (L[1] - Lx_constriction)/2., y], r_walls, h, 
       +                                   fixed=true, verbose=false)
        end
        
        dx = 2.*r_walls*sin(atan((Lx - Lx_constriction)/(L[2] - Ly_constriction)))
       t@@ -43,8 +44,8 @@ x = r_walls:dx:((Lx - Lx_constriction)/2.)
        y = linspace(L[2] - r_walls, (L[2] - Ly_constriction)/2. + Ly_constriction + 
                     r_walls, length(x))
        for i in 1:length(x)
       -    SeaIce.addIceFloeCylindrical(sim, [x[i], y[i]], r_walls, h, fixed=true, 
       -                                 verbose=false)
       +    SeaIce.addIceFloeCylindrical!(sim, [x[i], y[i]], r_walls, h, fixed=true, 
       +                                  verbose=false)
        end
        
        ## NE diagonal
       t@@ -52,24 +53,24 @@ x = (L[1] - r_walls):(-dx):((Lx - Lx_constriction)/2. + Lx_constriction)
        y = linspace(L[2] - r_walls, (L[2] - Ly_constriction)/2. + Ly_constriction + 
                     r_walls, length(x))
        for i in 1:length(x)
       -    SeaIce.addIceFloeCylindrical(sim, [x[i], y[i]], r_walls, h, fixed=true, 
       -                                 verbose=false)
       +    SeaIce.addIceFloeCylindrical!(sim, [x[i], y[i]], r_walls, h, fixed=true, 
       +                                  verbose=false)
        end
        
        ## SW diagonal
        x = r_walls:dx:((Lx - Lx_constriction)/2.)
        y = linspace(r, (L[2] - Ly_constriction)/2. - r_walls, length(x))
        for i in 1:length(x)
       -    SeaIce.addIceFloeCylindrical(sim, [x[i], y[i]], r_walls, h, fixed=true, 
       -                                 verbose=false)
       +    SeaIce.addIceFloeCylindrical!(sim, [x[i], y[i]], r_walls, h, fixed=true, 
       +                                  verbose=false)
        end
        
        ## SE diagonal
        x = (L[1] - r_walls):(-dx):((Lx - Lx_constriction)/2. + Lx_constriction)
        y = linspace(r_walls, (L[2] - Ly_constriction)/2. - r_walls, length(x))
        for i in 1:length(x)
       -    SeaIce.addIceFloeCylindrical(sim, [x[i], y[i]], r_walls, h, fixed=true, 
       -                                 verbose=false)
       +    SeaIce.addIceFloeCylindrical!(sim, [x[i], y[i]], r_walls, h, fixed=true, 
       +                                  verbose=false)
        end
        
        n_walls = length(sim.ice_floes)
       t@@ -102,7 +103,7 @@ for y in (L[2] - r - noise_amplitude):(-(2.*r + floe_padding)):((L[2] -
                end
                    
                r_rand = r_min + Base.Random.rand()*(r - r_min)
       -        SeaIce.addIceFloeCylindrical(sim, [x_, y_], r_rand, h, verbose=false)
       +        SeaIce.addIceFloeCylindrical!(sim, [x_, y_], r_rand, h, verbose=false)
            end
            iy += 1
        end
       t@@ -144,7 +145,7 @@ while sim.time < sim.time_total
        
                    # Enable for high mass flux
                    r_rand = r_min + Base.Random.rand()*(r - r_min)
       -            SeaIce.addIceFloeCylindrical(sim, [x-r, y-r], r_rand, h, 
       +            SeaIce.addIceFloeCylindrical!(sim, [x-r, y-r], r_rand, h, 
                            verbose=false,
                            contact_stiffness_normal=k_n,
                            contact_stiffness_tangential=k_t,
       t@@ -152,7 +153,7 @@ while sim.time < sim.time_total
                            contact_dynamic_friction = mu_d,
                            rotating=rotating)
                    r_rand = r_min + Base.Random.rand()*(r - r_min)
       -            SeaIce.addIceFloeCylindrical(sim, [x+r, y-r], r_rand, h, 
       +            SeaIce.addIceFloeCylindrical!(sim, [x+r, y-r], r_rand, h, 
                            verbose=false,
                            contact_stiffness_normal=k_n,
                            contact_stiffness_tangential=k_t,
       t@@ -160,7 +161,7 @@ while sim.time < sim.time_total
                            contact_dynamic_friction = mu_d,
                            rotating=rotating)
                    r_rand = r_min + Base.Random.rand()*(r - r_min)
       -            SeaIce.addIceFloeCylindrical(sim, [x+r, y+r], r_rand, h, 
       +            SeaIce.addIceFloeCylindrical!(sim, [x+r, y+r], r_rand, h, 
                            verbose=false,
                            contact_stiffness_normal=k_n,
                            contact_stiffness_tangential=k_t,
       t@@ -168,7 +169,7 @@ while sim.time < sim.time_total
                            contact_dynamic_friction = mu_d,
                            rotating=rotating)
                    r_rand = r_min + Base.Random.rand()*(r - r_min)
       -            SeaIce.addIceFloeCylindrical(sim, [x-r, y+r], r_rand, h, 
       +            SeaIce.addIceFloeCylindrical!(sim, [x-r, y+r], r_rand, h, 
                            verbose=false,
                            contact_stiffness_normal=k_n,
                            contact_stiffness_tangential=k_t,
       t@@ -179,7 +180,7 @@ while sim.time < sim.time_total
                    # Enable for low mass flux
                    #x += noise_amplitude*(0.5 - Base.Random.rand())
                    #y += noise_amplitude*(0.5 - Base.Random.rand())
       -            #SeaIce.addIceFloeCylindrical(sim, [x, y], r, h, verbose=false)
       +            #SeaIce.addIceFloeCylindrical!(sim, [x, y], r, h, verbose=false)
                end
            end
        end
 (DIR) diff --git a/src/contact_search.jl b/src/contact_search.jl
       t@@ -145,10 +145,10 @@ function checkAndAddContact!(sim::Simulation, i::Int, j::Int)
        
                # Check if grains overlap (overlap when negative)
                if overlap_ij < 0.
       -            for ic=1:(Nc_max + 1)
       -                if ic == (Nc_max + 1)
       +            for ic=1:(sim.Nc_max + 1)
       +                if ic == (sim.Nc_max + 1)
                            error("contact $i-$j exceeds max. number of contacts " *
       -                          "(Nc_max = $Nc_max) for ice floe $i")
       +                          "(sim.Nc_max = $(sim.Nc_max)) for ice floe $i")
        
                        else
                            if sim.ice_floes[i].contacts[ic] == j
 (DIR) diff --git a/src/datatypes.jl b/src/datatypes.jl
       t@@ -268,4 +268,6 @@ type Simulation
        
            ocean::Ocean
            atmosphere::Atmosphere
       +
       +    Nc_max::Int
        end
 (DIR) diff --git a/src/icefloe.jl b/src/icefloe.jl
       t@@ -1,63 +1,55 @@
        ## Manage icefloes in the model
        
       -const Nc_max = 32  # max. no. of contacts per ice floe
       -
       -export addIceFloeCylindrical
       +export addIceFloeCylindrical!
        """
        Adds a grain to the simulation. Example:
        
            SeaIce.addIceFloeCylindrical([1.0, 2.0, 3.0], 1.0)
        """
       -function addIceFloeCylindrical(simulation::Simulation,
       -                               lin_pos::vector,
       -                               contact_radius::float,
       -                               thickness::float;
       -                               areal_radius = false,
       -                               lin_vel::vector = [0., 0.],
       -                               lin_acc::vector = [0., 0.],
       -                               force::vector = [0., 0.],
       -                               ang_pos::float = 0.,
       -                               ang_vel::float = 0.,
       -                               ang_acc::float = 0.,
       -                               torque::float = 0.,
       -                               density::float = 934.,
       -                               contact_stiffness_normal::float = 1e7,
       -                               contact_stiffness_tangential::float = 0.,
       -                               contact_viscosity_normal::float = 0.,
       -                               contact_viscosity_tangential::float = 0.,
       -                               contact_static_friction::float = 0.4,
       -                               contact_dynamic_friction::float = 0.4,
       -                               youngs_modulus::float = 2e7,
       -                               #youngs_modulus::float = 2e9,  # Hopkins 2004
       -                               poissons_ratio::float = 0.185,  # Hopkins 2004
       -                               #tensile_strength::float = 500e3,  # Hopkins 2004
       -                               tensile_strength::float = 0.,
       -                               tensile_heal_rate::float = 0.,
       -                               compressive_strength_prefactor::float = 1285e3,  
       -                                   # Hopkins 2004
       -                               ocean_drag_coeff_vert::float = 0.85, # H&C 2011
       -                               ocean_drag_coeff_horiz::float = 5e-4, # H&C 2011
       -                               atmosphere_drag_coeff_vert::float = 0.4,
       -                                   # H&C 2011
       -                               atmosphere_drag_coeff_horiz::float = 2.5e-4,
       -                                   # H&C2011
       -                               pressure::float = 0.,
       -                               fixed::Bool = false,
       -                               rotating::Bool = true,
       -                               enabled::Bool = true,
       -                               verbose::Bool = true,
       -                               ocean_grid_pos::Array{Int, 1} = [0, 0],
       -                               atmosphere_grid_pos::Array{Int, 1} = [0, 0],
       -                               n_contacts::Int = 0,
       -                               contacts::Array{Int, 1} = zeros(Int, Nc_max),
       -                               contact_parallel_displacement::
       -                                   Array{Array{Float64, 1}, 1} =
       -                                   Array{Array{Float64, 1}, 1}(Nc_max),
       -                               contact_age::Array{Float64, 1} =
       -                                   zeros(Float64, Nc_max),
       -                               granular_stress::vector = [0., 0.],
       -                               ocean_stress::vector = [0., 0.],
       -                               atmosphere_stress::vector = [0., 0.])
       +function addIceFloeCylindrical!(simulation::Simulation,
       +                                lin_pos::vector,
       +                                contact_radius::float,
       +                                thickness::float;
       +                                areal_radius = false,
       +                                lin_vel::vector = [0., 0.],
       +                                lin_acc::vector = [0., 0.],
       +                                force::vector = [0., 0.],
       +                                ang_pos::float = 0.,
       +                                ang_vel::float = 0.,
       +                                ang_acc::float = 0.,
       +                                torque::float = 0.,
       +                                density::float = 934.,
       +                                contact_stiffness_normal::float = 1e7,
       +                                contact_stiffness_tangential::float = 0.,
       +                                contact_viscosity_normal::float = 0.,
       +                                contact_viscosity_tangential::float = 0.,
       +                                contact_static_friction::float = 0.4,
       +                                contact_dynamic_friction::float = 0.4,
       +                                youngs_modulus::float = 2e7,
       +                                #youngs_modulus::float = 2e9,  # Hopkins 2004
       +                                poissons_ratio::float = 0.185,  # Hopkins 2004
       +                                #tensile_strength::float = 500e3,  # Hopkins2004
       +                                tensile_strength::float = 0.,
       +                                tensile_heal_rate::float = 0.,
       +                                compressive_strength_prefactor::float = 1285e3,  
       +                                    # Hopkins 2004
       +                                ocean_drag_coeff_vert::float = 0.85, # H&C 2011
       +                                ocean_drag_coeff_horiz::float = 5e-4, # H&C 2011
       +                                atmosphere_drag_coeff_vert::float = 0.4,
       +                                    # H&C 2011
       +                                atmosphere_drag_coeff_horiz::float = 2.5e-4,
       +                                    # H&C2011
       +                                pressure::float = 0.,
       +                                fixed::Bool = false,
       +                                rotating::Bool = true,
       +                                enabled::Bool = true,
       +                                verbose::Bool = true,
       +                                ocean_grid_pos::Array{Int, 1} = [0, 0],
       +                                atmosphere_grid_pos::Array{Int, 1} = [0, 0],
       +                                n_contacts::Int = 0,
       +                                granular_stress::vector = [0., 0.],
       +                                ocean_stress::vector = [0., 0.],
       +                                atmosphere_stress::vector = [0., 0.])
        
            # Check input values
            if length(lin_pos) != 2
       t@@ -83,7 +75,11 @@ function addIceFloeCylindrical(simulation::Simulation,
                areal_radius = contact_radius
            end
        
       -    for i=1:Nc_max
       +    contacts::Array{Int, 1} = zeros(Int, simulation.Nc_max)
       +    contact_parallel_displacement =
       +        Array{Array{Float64, 1}, 1}(simulation.Nc_max)
       +    contact_age::Array{Float64, 1} = zeros(Float64, simulation.Nc_max)
       +    for i=1:simulation.Nc_max
                contact_parallel_displacement[i] = zeros(2)
            end
        
 (DIR) diff --git a/src/interaction.jl b/src/interaction.jl
       t@@ -8,7 +8,7 @@ Resolve mechanical interaction between all particle pairs.
        """
        function interact!(simulation::Simulation)
            for i=1:length(simulation.ice_floes)
       -        for ic=1:Nc_max
       +        for ic=1:simulation.Nc_max
        
                    j = simulation.ice_floes[i].contacts[ic]
        
 (DIR) diff --git a/src/io.jl b/src/io.jl
       t@@ -368,7 +368,7 @@ function writeIceFloeInteractionVTK(simulation::Simulation,
            shear_displacement = Array{float, 1}[]
            contact_age = float[]
            for i=1:length(simulation.ice_floes)
       -        for ic=1:Nc_max
       +        for ic=1:simulation.Nc_max
                    if simulation.ice_floes[i].contacts[ic] > 0
                        j = simulation.ice_floes[i].contacts[ic]
        
 (DIR) diff --git a/src/simulation.jl b/src/simulation.jl
       t@@ -29,7 +29,8 @@ function createSimulation(;id::String="unnamed",
                                  file_time_since_output_file::Float64=0.,
                                  ice_floes=Array{IceFloeCylindrical, 1}[],
                                  ocean::Ocean=createEmptyOcean(),
       -                          atmosphere::Atmosphere=createEmptyAtmosphere())
       +                          atmosphere::Atmosphere=createEmptyAtmosphere(),
       +                          Nc_max::Int=16)
        
            return Simulation(id,
                              time_iteration,
       t@@ -41,7 +42,8 @@ function createSimulation(;id::String="unnamed",
                              file_time_since_output_file,
                              ice_floes,
                              ocean,
       -                      atmosphere)
       +                      atmosphere,
       +                      Nc_max)
        end
        
        export run!
       t@@ -261,4 +263,6 @@ function compareSimulations(sim1::Simulation, sim2::Simulation)
            end
            compareOceans(sim1.ocean, sim2.ocean)
            compareAtmospheres(sim1.atmosphere, sim2.atmosphere)
       +
       +    Base.Test.@test sim1.Nc_max == sim2.Nc_max
        end
 (DIR) diff --git a/test/atmosphere.jl b/test/atmosphere.jl
       t@@ -24,7 +24,7 @@ sim.atmosphere = SeaIce.createRegularAtmosphereGrid([6, 6, 6], [1., 1., 1.])
        @test sim.atmosphere.v ≈ zeros(7, 7, 6, 1)
        
        info("Testing velocity drag interaction (static atmosphere)")
       -SeaIce.addIceFloeCylindrical(sim, [.5, .5], .25, .1)
       +SeaIce.addIceFloeCylindrical!(sim, [.5, .5], .25, .1)
        SeaIce.setTotalTime!(sim, 5.)
        SeaIce.setTimeStep!(sim)
        sim_init = deepcopy(sim)
 (DIR) diff --git a/test/cohesion.jl b/test/cohesion.jl
       t@@ -8,8 +8,8 @@ info("#### $(basename(@__FILE__)) ####")
        verbose=false
        
        sim_init = SeaIce.createSimulation()
       -SeaIce.addIceFloeCylindrical(sim_init, [0., 0.], 10., 1.)
       -SeaIce.addIceFloeCylindrical(sim_init, [18., 0.], 10., 1.)
       +SeaIce.addIceFloeCylindrical!(sim_init, [0., 0.], 10., 1.)
       +SeaIce.addIceFloeCylindrical!(sim_init, [18., 0.], 10., 1.)
        sim_init.ice_floes[1].youngs_modulus = 1e-5  # repulsion is negligible
        sim_init.ice_floes[2].youngs_modulus = 1e-5  # repulsion is negligible
        SeaIce.setTimeStep!(sim_init, verbose=verbose)
       t@@ -23,8 +23,8 @@ SeaIce.run!(sim, verbose=verbose)
        
        info("# Check if bonds add tensile strength")
        sim = SeaIce.createSimulation(id="cohesion")
       -SeaIce.addIceFloeCylindrical(sim, [0., 0.], 10., 1., tensile_strength=500e3)
       -SeaIce.addIceFloeCylindrical(sim, [20.1, 0.], 10., 1., tensile_strength=500e3)
       +SeaIce.addIceFloeCylindrical!(sim, [0., 0.], 10., 1., tensile_strength=500e3)
       +SeaIce.addIceFloeCylindrical!(sim, [20.1, 0.], 10., 1., tensile_strength=500e3)
        sim.ice_floes[1].lin_vel[1] = 0.1
        SeaIce.setTimeStep!(sim)
        SeaIce.setTotalTime!(sim, 10.)
 (DIR) diff --git a/test/collision-2floes-normal.jl b/test/collision-2floes-normal.jl
       t@@ -9,8 +9,8 @@ verbose=false
        
        info("# One ice floe fixed")
        sim = SeaIce.createSimulation(id="test")
       -SeaIce.addIceFloeCylindrical(sim, [0., 0.], 10., 1., verbose=verbose)
       -SeaIce.addIceFloeCylindrical(sim, [20.05, 0.], 10., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [0., 0.], 10., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [20.05, 0.], 10., 1., verbose=verbose)
        sim.ice_floes[1].lin_vel[1] = 0.1
        sim.ice_floes[2].fixed = true
        
       t@@ -65,8 +65,8 @@ E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
        info("# Ice floes free to move")
        
        sim = SeaIce.createSimulation(id="test")
       -SeaIce.addIceFloeCylindrical(sim, [0., 0.], 10., 1., verbose=verbose)
       -SeaIce.addIceFloeCylindrical(sim, [20.05, 0.], 10., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [0., 0.], 10., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [20.05, 0.], 10., 1., verbose=verbose)
        sim.ice_floes[1].lin_vel[1] = 0.1
        
        E_kin_lin_init = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
       t@@ -120,8 +120,8 @@ E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
        info("# Adding contact-normal viscosity")
        info("# One ice floe fixed")
        sim = SeaIce.createSimulation(id="test")
       -SeaIce.addIceFloeCylindrical(sim, [0., 0.], 10., 1., verbose=verbose)
       -SeaIce.addIceFloeCylindrical(sim, [20.05, 0.], 10., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [0., 0.], 10., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [20.05, 0.], 10., 1., verbose=verbose)
        sim.ice_floes[1].lin_vel[1] = 0.1
        sim.ice_floes[1].contact_viscosity_normal = 1e4
        sim.ice_floes[2].contact_viscosity_normal = 1e4
       t@@ -167,8 +167,8 @@ E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
        info("# Ice floes free to move")
        
        sim = SeaIce.createSimulation(id="test")
       -SeaIce.addIceFloeCylindrical(sim, [0., 0.], 10., 1., verbose=verbose)
       -SeaIce.addIceFloeCylindrical(sim, [20.05, 0.], 10., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [0., 0.], 10., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [20.05, 0.], 10., 1., verbose=verbose)
        sim.ice_floes[1].lin_vel[1] = 0.1
        sim.ice_floes[1].contact_viscosity_normal = 1e4
        sim.ice_floes[2].contact_viscosity_normal = 1e4
 (DIR) diff --git a/test/collision-2floes-oblique.jl b/test/collision-2floes-oblique.jl
       t@@ -10,8 +10,8 @@ verbose=false
        info("## Contact-normal elasticity only")
        info("# One ice floe fixed")
        sim = SeaIce.createSimulation(id="test")
       -SeaIce.addIceFloeCylindrical(sim, [0., 10.], 10., 1., verbose=verbose)
       -SeaIce.addIceFloeCylindrical(sim, [19., 0.], 10., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [0., 10.], 10., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [19., 0.], 10., 1., verbose=verbose)
        sim.ice_floes[1].lin_vel[1] = 0.1
        sim.ice_floes[1].contact_dynamic_friction = 0.
        sim.ice_floes[2].contact_dynamic_friction = 0.
       t@@ -67,8 +67,8 @@ E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
        info("# Ice floes free to move")
        
        sim = SeaIce.createSimulation(id="test")
       -SeaIce.addIceFloeCylindrical(sim, [0., 10.], 10., 1., verbose=verbose)
       -SeaIce.addIceFloeCylindrical(sim, [19.0, 0.], 10., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [0., 10.], 10., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [19.0, 0.], 10., 1., verbose=verbose)
        sim.ice_floes[1].lin_vel[1] = 0.1
        sim.ice_floes[1].contact_dynamic_friction = 0.
        sim.ice_floes[2].contact_dynamic_friction = 0.
       t@@ -200,8 +200,8 @@ E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
        info("# Ice floes free to move")
        
        sim = SeaIce.createSimulation(id="test")
       -SeaIce.addIceFloeCylindrical(sim, [0., 10.], 10., 1., verbose=verbose)
       -SeaIce.addIceFloeCylindrical(sim, [19.0, 0.], 10., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [0., 10.], 10., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [19.0, 0.], 10., 1., verbose=verbose)
        sim.ice_floes[1].lin_vel[1] = 0.1
        sim.ice_floes[1].contact_viscosity_tangential = 1e4
        sim.ice_floes[2].contact_viscosity_tangential = 1e4
       t@@ -263,8 +263,8 @@ E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
        info("# Ice floes free to move, mirrored")
        
        sim = SeaIce.createSimulation(id="test")
       -SeaIce.addIceFloeCylindrical(sim, [0., 0.], 10., 1., verbose=verbose)
       -SeaIce.addIceFloeCylindrical(sim, [19.0, 10.], 10., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [0., 0.], 10., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [19.0, 10.], 10., 1., verbose=verbose)
        sim.ice_floes[2].lin_vel[1] = -0.1
        sim.ice_floes[1].contact_viscosity_tangential = 1e4
        sim.ice_floes[2].contact_viscosity_tangential = 1e4
       t@@ -326,8 +326,8 @@ E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
        info("# Ice floes free to move, mirrored #2")
        
        sim = SeaIce.createSimulation(id="test")
       -SeaIce.addIceFloeCylindrical(sim, [0., 0.], 10., 1., verbose=verbose)
       -SeaIce.addIceFloeCylindrical(sim, [19.0, -10.], 10., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [0., 0.], 10., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [19.0, -10.], 10., 1., verbose=verbose)
        sim.ice_floes[2].lin_vel[1] = -0.1
        
        E_kin_lin_init = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
       t@@ -387,8 +387,8 @@ E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
        info("# Tangential elasticity, no tangential viscosity, no Coulomb slip")
        
        sim = SeaIce.createSimulation(id="test")
       -SeaIce.addIceFloeCylindrical(sim, [0., 0.], 10., 1., verbose=verbose)
       -SeaIce.addIceFloeCylindrical(sim, [19.0, -10.], 10., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [0., 0.], 10., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [19.0, -10.], 10., 1., verbose=verbose)
        sim.ice_floes[2].lin_vel[1] = -0.1
        sim.ice_floes[1].contact_dynamic_friction = 1e3  # disable Coulomb slip
        sim.ice_floes[2].contact_dynamic_friction = 1e3  # disable Coulomb slip
       t@@ -456,8 +456,8 @@ E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
        info("# Tangential elasticity, no tangential viscosity, Coulomb slip")
        
        sim = SeaIce.createSimulation(id="test")
       -SeaIce.addIceFloeCylindrical(sim, [0., 0.], 10., 1., verbose=verbose)
       -SeaIce.addIceFloeCylindrical(sim, [19.0, -10.], 10., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [0., 0.], 10., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [19.0, -10.], 10., 1., verbose=verbose)
        sim.ice_floes[2].lin_vel[1] = -0.1
        sim.ice_floes[1].contact_dynamic_friction = 0.1  # enable Coulomb slip
        sim.ice_floes[2].contact_dynamic_friction = 0.1  # enable Coulomb slip
       t@@ -509,8 +509,8 @@ E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
        info("# Tangential elasticity, tangential viscosity, no Coulomb slip")
        
        sim = SeaIce.createSimulation(id="test")
       -SeaIce.addIceFloeCylindrical(sim, [0., 0.], 10., 1., verbose=verbose)
       -SeaIce.addIceFloeCylindrical(sim, [19.0, -10.], 10., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [0., 0.], 10., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [19.0, -10.], 10., 1., verbose=verbose)
        sim.ice_floes[2].lin_vel[1] = -0.1
        sim.ice_floes[1].contact_dynamic_friction = 1e3  # disable Coulomb slip
        sim.ice_floes[2].contact_dynamic_friction = 1e3  # disable Coulomb slip
       t@@ -562,8 +562,8 @@ E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
        info("# Tangential elasticity, tangential viscosity, Coulomb slip")
        
        sim = SeaIce.createSimulation(id="test")
       -SeaIce.addIceFloeCylindrical(sim, [0., 0.], 10., 1., verbose=verbose)
       -SeaIce.addIceFloeCylindrical(sim, [19.0, -10.], 10., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [0., 0.], 10., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [19.0, -10.], 10., 1., verbose=verbose)
        sim.ice_floes[2].lin_vel[1] = -0.1
        sim.ice_floes[1].contact_dynamic_friction = 0.1  # enable Coulomb slip
        sim.ice_floes[2].contact_dynamic_friction = 0.1  # enable Coulomb slip
 (DIR) diff --git a/test/collision-5floes-normal.jl b/test/collision-5floes-normal.jl
       t@@ -9,11 +9,11 @@ verbose=false
        
        info("# One ice floe fixed")
        sim = SeaIce.createSimulation(id="test")
       -SeaIce.addIceFloeCylindrical(sim, [0., 0.], 10., 1., verbose=verbose)
       -SeaIce.addIceFloeCylindrical(sim, [20.05, 0.], 10., 1., verbose=verbose)
       -SeaIce.addIceFloeCylindrical(sim, [40.05, 0.], 10., 1., verbose=verbose)
       -SeaIce.addIceFloeCylindrical(sim, [60.05, 0.], 10., 1., verbose=verbose)
       -SeaIce.addIceFloeCylindrical(sim, [80.05, 0.], 10., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [0., 0.], 10., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [20.05, 0.], 10., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [40.05, 0.], 10., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [60.05, 0.], 10., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [80.05, 0.], 10., 1., verbose=verbose)
        sim.ice_floes[1].lin_vel[1] = 0.1
        sim.ice_floes[2].fixed = true
        sim.ice_floes[3].fixed = true
       t@@ -86,11 +86,11 @@ end
        info("# Ice floes free to move")
        
        sim = SeaIce.createSimulation(id="test")
       -SeaIce.addIceFloeCylindrical(sim, [0., 0.], 10., 1., verbose=verbose)
       -SeaIce.addIceFloeCylindrical(sim, [20.05, 0.], 10., 1., verbose=verbose)
       -SeaIce.addIceFloeCylindrical(sim, [40.05, 0.], 10., 1., verbose=verbose)
       -SeaIce.addIceFloeCylindrical(sim, [60.05, 0.], 10., 1., verbose=verbose)
       -SeaIce.addIceFloeCylindrical(sim, [80.05, 0.], 10., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [0., 0.], 10., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [20.05, 0.], 10., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [40.05, 0.], 10., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [60.05, 0.], 10., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [80.05, 0.], 10., 1., verbose=verbose)
        sim.ice_floes[1].lin_vel[1] = 0.1
        
        E_kin_lin_init = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
       t@@ -156,11 +156,11 @@ end
        info("# Adding contact-normal viscosity")
        info("# One ice floe fixed")
        sim = SeaIce.createSimulation(id="test")
       -SeaIce.addIceFloeCylindrical(sim, [0., 0.], 10., 1., verbose=verbose)
       -SeaIce.addIceFloeCylindrical(sim, [20.05, 0.], 10., 1., verbose=verbose)
       -SeaIce.addIceFloeCylindrical(sim, [40.05, 0.], 10., 1., verbose=verbose)
       -SeaIce.addIceFloeCylindrical(sim, [60.05, 0.], 10., 1., verbose=verbose)
       -SeaIce.addIceFloeCylindrical(sim, [80.05, 0.], 10., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [0., 0.], 10., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [20.05, 0.], 10., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [40.05, 0.], 10., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [60.05, 0.], 10., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [80.05, 0.], 10., 1., verbose=verbose)
        sim.ice_floes[1].lin_vel[1] = 0.1
        sim.ice_floes[1].contact_viscosity_normal = 1e4
        sim.ice_floes[2].contact_viscosity_normal = 1e4
       t@@ -222,11 +222,11 @@ end
        info("# Ice floes free to move")
        
        sim = SeaIce.createSimulation(id="test")
       -SeaIce.addIceFloeCylindrical(sim, [0., 0.], 10., 1., verbose=verbose)
       -SeaIce.addIceFloeCylindrical(sim, [20.05, 0.], 10., 1., verbose=verbose)
       -SeaIce.addIceFloeCylindrical(sim, [40.05, 0.], 10., 1., verbose=verbose)
       -SeaIce.addIceFloeCylindrical(sim, [60.05, 0.], 10., 1., verbose=verbose)
       -SeaIce.addIceFloeCylindrical(sim, [80.05, 0.], 10., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [0., 0.], 10., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [20.05, 0.], 10., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [40.05, 0.], 10., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [60.05, 0.], 10., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [80.05, 0.], 10., 1., verbose=verbose)
        sim.ice_floes[1].lin_vel[1] = 0.1
        sim.ice_floes[1].contact_viscosity_normal = 1e4
        sim.ice_floes[2].contact_viscosity_normal = 1e4
 (DIR) diff --git a/test/contact-search-and-geometry.jl b/test/contact-search-and-geometry.jl
       t@@ -6,8 +6,8 @@ info("#### $(basename(@__FILE__)) ####")
        
        info("Testing interIceFloePositionVector(...) and findOverlap(...)")
        sim = SeaIce.createSimulation(id="test")
       -SeaIce.addIceFloeCylindrical(sim, [ 0., 0.], 10., 1., verbose=false)
       -SeaIce.addIceFloeCylindrical(sim, [18., 0.], 10., 1., verbose=false)
       +SeaIce.addIceFloeCylindrical!(sim, [ 0., 0.], 10., 1., verbose=false)
       +SeaIce.addIceFloeCylindrical!(sim, [18., 0.], 10., 1., verbose=false)
        
        position_ij = SeaIce.interIceFloePositionVector(sim, 1, 2)
        overlap_ij = SeaIce.findOverlap(sim, 1, 2, position_ij)
       t@@ -28,11 +28,11 @@ SeaIce.findContacts!(sim)
        sim.ice_floes[1].fixed = true
        # The contact should be registered in ice floe 1, but not ice floe 2
        @test 2 == sim.ice_floes[1].contacts[1]
       -for ic=2:32
       +for ic=2:sim.Nc_max
            @test 0 == sim.ice_floes[1].contacts[ic]
            @test [0., 0.] ≈ sim.ice_floes[1].contact_parallel_displacement[ic]
        end
       -for ic=1:32
       +for ic=1:sim.Nc_max
            @test 0 == sim.ice_floes[2].contacts[ic]
            @test [0., 0.] ≈ sim.ice_floes[2].contact_parallel_displacement[ic]
        end
       t@@ -44,11 +44,11 @@ sim = deepcopy(sim_copy)
        SeaIce.findContacts!(sim)
        
        @test 2 == sim.ice_floes[1].contacts[1]
       -for ic=2:32
       +for ic=2:sim.Nc_max
            @test 0 == sim.ice_floes[1].contacts[ic]
            @test [0., 0.] ≈ sim.ice_floes[1].contact_parallel_displacement[ic]
        end
       -for ic=1:32
       +for ic=1:sim.Nc_max
            @test 0 == sim.ice_floes[2].contacts[ic]
            @test [0., 0.] ≈ sim.ice_floes[2].contact_parallel_displacement[ic]
        end
       t@@ -61,11 +61,11 @@ sim = deepcopy(sim_copy)
        sim.ice_floes[1].fixed = true
        sim.ice_floes[2].fixed = true
        SeaIce.findContacts!(sim)
       -for ic=1:32
       +for ic=1:sim.Nc_max
            @test 0 == sim.ice_floes[1].contacts[ic]
            @test [0., 0.] ≈ sim.ice_floes[1].contact_parallel_displacement[ic]
        end
       -for ic=1:32
       +for ic=1:sim.Nc_max
            @test 0 == sim.ice_floes[2].contacts[ic]
            @test [0., 0.] ≈ sim.ice_floes[2].contact_parallel_displacement[ic]
        end
       t@@ -76,11 +76,11 @@ end
        sim = deepcopy(sim_copy)
        SeaIce.disableIceFloe!(sim, 1)
        SeaIce.findContacts!(sim)
       -for ic=1:32
       +for ic=1:sim.Nc_max
            @test 0 == sim.ice_floes[1].contacts[ic]
            @test [0., 0.] ≈ sim.ice_floes[1].contact_parallel_displacement[ic]
        end
       -for ic=1:32
       +for ic=1:sim.Nc_max
            @test 0 == sim.ice_floes[2].contacts[ic]
            @test [0., 0.] ≈ sim.ice_floes[2].contact_parallel_displacement[ic]
        end
       t@@ -92,11 +92,11 @@ sim = deepcopy(sim_copy)
        SeaIce.disableIceFloe!(sim, 1)
        SeaIce.disableIceFloe!(sim, 2)
        SeaIce.findContacts!(sim)
       -for ic=1:32
       +for ic=1:sim.Nc_max
            @test 0 == sim.ice_floes[1].contacts[ic]
            @test [0., 0.] ≈ sim.ice_floes[1].contact_parallel_displacement[ic]
        end
       -for ic=1:32
       +for ic=1:sim.Nc_max
            @test 0 == sim.ice_floes[2].contacts[ic]
            @test [0., 0.] ≈ sim.ice_floes[2].contact_parallel_displacement[ic]
        end
       t@@ -110,11 +110,11 @@ SeaIce.interact!(sim)
        SeaIce.findContacts!(sim)
        
        @test 2 == sim.ice_floes[1].contacts[1]
       -for ic=2:32
       +for ic=2:sim.Nc_max
            @test 0 == sim.ice_floes[1].contacts[ic]
            @test [0., 0.] ≈ sim.ice_floes[1].contact_parallel_displacement[ic]
        end
       -for ic=1:32
       +for ic=1:sim.Nc_max
            @test 0 == sim.ice_floes[2].contacts[ic]
            @test [0., 0.] ≈ sim.ice_floes[2].contact_parallel_displacement[ic]
        end
       t@@ -129,11 +129,11 @@ SeaIce.sortIceFloesInGrid!(sim, sim.ocean)
        SeaIce.findContactsInGrid!(sim, sim.ocean)
        
        @test 2 == sim.ice_floes[1].contacts[1]
       -for ic=2:32
       +for ic=2:sim.Nc_max
            @test 0 == sim.ice_floes[1].contacts[ic]
            @test [0., 0.] ≈ sim.ice_floes[1].contact_parallel_displacement[ic]
        end
       -for ic=1:32
       +for ic=1:sim.Nc_max
            @test 0 == sim.ice_floes[2].contacts[ic]
            @test [0., 0.] ≈ sim.ice_floes[2].contact_parallel_displacement[ic]
        end
       t@@ -148,11 +148,11 @@ SeaIce.sortIceFloesInGrid!(sim, sim.ocean)
        SeaIce.findContactsInGrid!(sim, sim.ocean)
        
        @test 2 == sim.ice_floes[1].contacts[1]
       -for ic=2:32
       +for ic=2:sim.Nc_max
            @test 0 == sim.ice_floes[1].contacts[ic]
            @test [0., 0.] ≈ sim.ice_floes[1].contact_parallel_displacement[ic]
        end
       -for ic=1:32
       +for ic=1:sim.Nc_max
            @test 0 == sim.ice_floes[2].contacts[ic]
            @test [0., 0.] ≈ sim.ice_floes[2].contact_parallel_displacement[ic]
        end
       t@@ -167,11 +167,11 @@ sim.ice_floes[2].fixed = true
        SeaIce.sortIceFloesInGrid!(sim, sim.ocean)
        SeaIce.findContactsInGrid!(sim, sim.ocean)
        
       -for ic=1:32
       +for ic=1:sim.Nc_max
            @test 0 == sim.ice_floes[1].contacts[ic]
            @test [0., 0.] ≈ sim.ice_floes[1].contact_parallel_displacement[ic]
        end
       -for ic=1:32
       +for ic=1:sim.Nc_max
            @test 0 == sim.ice_floes[2].contacts[ic]
            @test [0., 0.] ≈ sim.ice_floes[2].contact_parallel_displacement[ic]
        end
       t@@ -185,11 +185,11 @@ SeaIce.sortIceFloesInGrid!(sim, sim.ocean)
        SeaIce.findContacts!(sim)
        
        @test 2 == sim.ice_floes[1].contacts[1]
       -for ic=2:32
       +for ic=2:sim.Nc_max
            @test 0 == sim.ice_floes[1].contacts[ic]
            @test [0., 0.] ≈ sim.ice_floes[1].contact_parallel_displacement[ic]
        end
       -for ic=1:32
       +for ic=1:sim.Nc_max
            @test 0 == sim.ice_floes[2].contacts[ic]
            @test [0., 0.] ≈ sim.ice_floes[2].contact_parallel_displacement[ic]
        end
 (DIR) diff --git a/test/grid.jl b/test/grid.jl
       t@@ -103,9 +103,9 @@ info("Testing cell binning - Conformal mapping")
        sim = SeaIce.createSimulation()
        sim.ocean = SeaIce.readOceanNetCDF("Baltic/00010101.ocean_month.nc",
                                           "Baltic/ocean_hgrid.nc")
       -SeaIce.addIceFloeCylindrical(sim, [6.5, 53.5], 10., 1., verbose=verbose)
       -SeaIce.addIceFloeCylindrical(sim, [6.6, 53.5], 10., 1., verbose=verbose)
       -SeaIce.addIceFloeCylindrical(sim, [7.5, 53.5], 10., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [6.5, 53.5], 10., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [6.6, 53.5], 10., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [7.5, 53.5], 10., 1., verbose=verbose)
        SeaIce.sortIceFloesInGrid!(sim, sim.ocean, verbose=verbose)
        @test sim.ice_floes[1].ocean_grid_pos == [1, 1]
        @test sim.ice_floes[2].ocean_grid_pos == [1, 1]
       t@@ -117,8 +117,8 @@ info("Testing ocean drag")
        sim = SeaIce.createSimulation()
        sim.ocean = SeaIce.createRegularOceanGrid([4, 4, 2], [4., 4., 2.])
        sim.ocean.u[:,:,1,1] = 5.
       -SeaIce.addIceFloeCylindrical(sim, [2.5, 3.5], 1., 1., verbose=verbose)
       -SeaIce.addIceFloeCylindrical(sim, [2.6, 2.5], 1., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [2.5, 3.5], 1., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [2.6, 2.5], 1., 1., verbose=verbose)
        SeaIce.sortIceFloesInGrid!(sim, sim.ocean, verbose=verbose)
        sim.time = ocean.time[1]
        SeaIce.addOceanDrag!(sim)
       t@@ -128,8 +128,8 @@ SeaIce.addOceanDrag!(sim)
        @test sim.ice_floes[2].force[2] ≈ 0.
        sim.ocean.u[:,:,1,1] = -5.
        sim.ocean.v[:,:,1,1] = 5.
       -SeaIce.addIceFloeCylindrical(sim, [2.5, 3.5], 1., 1., verbose=verbose)
       -SeaIce.addIceFloeCylindrical(sim, [2.6, 2.5], 1., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [2.5, 3.5], 1., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [2.6, 2.5], 1., 1., verbose=verbose)
        SeaIce.sortIceFloesInGrid!(sim, sim.ocean, verbose=verbose)
        sim.time = ocean.time[1]
        SeaIce.addOceanDrag!(sim)
       t@@ -158,8 +158,8 @@ sim = SeaIce.createSimulation()
        sim.atmosphere = SeaIce.createRegularAtmosphereGrid([4, 4, 2], [4., 4., 2.])
        atmosphere = SeaIce.createRegularAtmosphereGrid([4, 4, 2], [4., 4., 2.])
        sim.atmosphere.u[:,:,1,1] = 5.
       -SeaIce.addIceFloeCylindrical(sim, [2.5, 3.5], 1., 1., verbose=verbose)
       -SeaIce.addIceFloeCylindrical(sim, [2.6, 2.5], 1., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [2.5, 3.5], 1., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [2.6, 2.5], 1., 1., verbose=verbose)
        SeaIce.sortIceFloesInGrid!(sim, sim.atmosphere, verbose=verbose)
        sim.time = ocean.time[1]
        SeaIce.addAtmosphereDrag!(sim)
       t@@ -169,8 +169,8 @@ SeaIce.addAtmosphereDrag!(sim)
        @test sim.ice_floes[2].force[2] ≈ 0.
        sim.atmosphere.u[:,:,1,1] = -5.
        sim.atmosphere.v[:,:,1,1] = 5.
       -SeaIce.addIceFloeCylindrical(sim, [2.5, 3.5], 1., 1., verbose=verbose)
       -SeaIce.addIceFloeCylindrical(sim, [2.6, 2.5], 1., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [2.5, 3.5], 1., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [2.6, 2.5], 1., 1., verbose=verbose)
        SeaIce.sortIceFloesInGrid!(sim, sim.atmosphere, verbose=verbose)
        sim.time = ocean.time[1]
        SeaIce.addAtmosphereDrag!(sim)
       t@@ -208,7 +208,7 @@ pos = SeaIce.findEmptyPositionInGridCell(sim, sim.ocean, 1, 1, 0.5,
        info("# Insert into cell with one other ice floe")
        sim = SeaIce.createSimulation()
        sim.ocean = SeaIce.createRegularOceanGrid([4, 4, 2], [4., 4., 2.])
       -SeaIce.addIceFloeCylindrical(sim, [.25, .25], .25, 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [.25, .25], .25, 1., verbose=verbose)
        SeaIce.sortIceFloesInGrid!(sim, sim.ocean, verbose=verbose)
        pos = SeaIce.findEmptyPositionInGridCell(sim, sim.ocean, 1, 1, .25, 
                                                 verbose=true)
       t@@ -218,8 +218,8 @@ pos = SeaIce.findEmptyPositionInGridCell(sim, sim.ocean, 1, 1, .25,
        info("# Insert into cell with two other ice floes")
        sim = SeaIce.createSimulation()
        sim.ocean = SeaIce.createRegularOceanGrid([4, 4, 2], [4., 4., 2.])
       -SeaIce.addIceFloeCylindrical(sim, [.25, .25], .25, 1., verbose=verbose)
       -SeaIce.addIceFloeCylindrical(sim, [.75, .75], .25, 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [.25, .25], .25, 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [.75, .75], .25, 1., verbose=verbose)
        SeaIce.sortIceFloesInGrid!(sim, sim.ocean, verbose=verbose)
        pos = SeaIce.findEmptyPositionInGridCell(sim, sim.ocean, 1, 1, .25, 
                                                 verbose=true)
       t@@ -229,10 +229,10 @@ pos = SeaIce.findEmptyPositionInGridCell(sim, sim.ocean, 1, 1, .25,
        info("# Insert into full cell")
        sim = SeaIce.createSimulation()
        sim.ocean = SeaIce.createRegularOceanGrid([4, 4, 2], [4., 4., 2.])
       -SeaIce.addIceFloeCylindrical(sim, [.5, .5], 1., 1., verbose=verbose)
       -SeaIce.addIceFloeCylindrical(sim, [.75, .5], 1., 1., verbose=verbose)
       -SeaIce.addIceFloeCylindrical(sim, [.5, .75], 1., 1., verbose=verbose)
       -SeaIce.addIceFloeCylindrical(sim, [.75, .75], 1., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [.5, .5], 1., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [.75, .5], 1., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [.5, .75], 1., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [.75, .75], 1., 1., verbose=verbose)
        SeaIce.sortIceFloesInGrid!(sim, sim.ocean, verbose=verbose)
        pos = SeaIce.findEmptyPositionInGridCell(sim, sim.ocean, 1, 1, 0.5,
                                                 verbose=false)
       t@@ -250,10 +250,10 @@ pos = SeaIce.findEmptyPositionInGridCell(sim, sim.ocean, 2, 2, 0.5,
        info("# Insert into full cell")
        sim = SeaIce.createSimulation()
        sim.ocean = SeaIce.createRegularOceanGrid([4, 4, 2], [4., 4., 2.])
       -SeaIce.addIceFloeCylindrical(sim, [1.5, 1.5], 1., 1., verbose=verbose)
       -SeaIce.addIceFloeCylindrical(sim, [1.75, 1.5], 1., 1., verbose=verbose)
       -SeaIce.addIceFloeCylindrical(sim, [1.5, 1.75], 1., 1., verbose=verbose)
       -SeaIce.addIceFloeCylindrical(sim, [1.75, 1.75], 1., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [1.5, 1.5], 1., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [1.75, 1.5], 1., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [1.5, 1.75], 1., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [1.75, 1.75], 1., 1., verbose=verbose)
        SeaIce.sortIceFloesInGrid!(sim, sim.ocean, verbose=verbose)
        pos = SeaIce.findEmptyPositionInGridCell(sim, sim.ocean, 2, 2, 0.5,
                                                 verbose=false)
       t@@ -263,9 +263,9 @@ info("Test default sorting with ocean/atmosphere grids")
        sim = SeaIce.createSimulation()
        sim.ocean = SeaIce.createRegularOceanGrid([4, 4, 2], [4., 4., 2.])
        sim.atmosphere = SeaIce.createRegularAtmosphereGrid([4, 4, 2], [4., 4.000001, 2.])
       -SeaIce.addIceFloeCylindrical(sim, [0.5, 0.5], .1, 1., verbose=verbose)
       -SeaIce.addIceFloeCylindrical(sim, [0.7, 0.7], .1, 1., verbose=verbose)
       -SeaIce.addIceFloeCylindrical(sim, [2.6, 2.5], .1, 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [0.5, 0.5], .1, 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [0.7, 0.7], .1, 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [2.6, 2.5], .1, 1., verbose=verbose)
        SeaIce.sortIceFloesInGrid!(sim, sim.ocean, verbose=verbose)
        SeaIce.setTimeStep!(sim)
        SeaIce.setTotalTime!(sim, 1.0)
       t@@ -288,9 +288,9 @@ info("Test optimization when ocean/atmosphere grids are collocated")
        sim = SeaIce.createSimulation()
        sim.ocean = SeaIce.createRegularOceanGrid([4, 4, 2], [4., 4., 2.])
        sim.atmosphere = SeaIce.createRegularAtmosphereGrid([4, 4, 2], [4., 4., 2.])
       -SeaIce.addIceFloeCylindrical(sim, [0.5, 0.5], .1, 1., verbose=verbose)
       -SeaIce.addIceFloeCylindrical(sim, [0.7, 0.7], .1, 1., verbose=verbose)
       -SeaIce.addIceFloeCylindrical(sim, [2.6, 2.5], .1, 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [0.5, 0.5], .1, 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [0.7, 0.7], .1, 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [2.6, 2.5], .1, 1., verbose=verbose)
        SeaIce.sortIceFloesInGrid!(sim, sim.ocean, verbose=verbose)
        SeaIce.setTimeStep!(sim)
        SeaIce.setTotalTime!(sim, 1.0)
       t@@ -312,5 +312,5 @@ SeaIce.run!(sim, single_step=true, verbose=verbose)
        info("Testing ocean drag")
        sim = SeaIce.createSimulation()
        sim.ocean.u[:,:,1,1] = 5.
       -SeaIce.addIceFloeCylindrical(sim, [2.5, 3.5], 1., 1., verbose=verbose)
       -SeaIce.addIceFloeCylindrical(sim, [2.6, 2.5], 1., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [2.5, 3.5], 1., 1., verbose=verbose)
       +SeaIce.addIceFloeCylindrical!(sim, [2.6, 2.5], 1., 1., verbose=verbose)
 (DIR) diff --git a/test/icefloe.jl b/test/icefloe.jl
       t@@ -6,21 +6,25 @@ info("#### $(basename(@__FILE__)) ####")
        
        info("Writing simple simulation to VTK file")
        sim = SeaIce.createSimulation(id="test")
       -SeaIce.addIceFloeCylindrical(sim, [ 0., 0.], 10., 1., verbose=false)
       +SeaIce.addIceFloeCylindrical!(sim, [ 0., 0.], 10., 1., verbose=false)
        SeaIce.printIceFloeInfo(sim.ice_floes[1])
        
        
       -@test_throws ErrorException SeaIce.addIceFloeCylindrical(sim, [.1, .1, .1], 10., 1.)
       -@test_throws ErrorException SeaIce.addIceFloeCylindrical(sim, [.1, .1], 10., 1., 
       -    lin_vel=[.2,.2,.2])
       -@test_throws ErrorException SeaIce.addIceFloeCylindrical(sim, [.1, .1], 10., 1., 
       -    lin_acc=[.2,.2,.2])
       -@test_throws ErrorException SeaIce.addIceFloeCylindrical(sim, [.1, .1], 0., 1.)
       -@test_throws ErrorException SeaIce.addIceFloeCylindrical(sim, [.1, .1], 10., 1., 
       -    density=-2.)
       +@test_throws ErrorException SeaIce.addIceFloeCylindrical!(sim, [.1, .1, .1],
       +                                                          10., 1.)
       +@test_throws ErrorException SeaIce.addIceFloeCylindrical!(sim, [.1, .1],
       +                                                          10., 1., 
       +                                                          lin_vel=[.2,.2,.2])
       +@test_throws ErrorException SeaIce.addIceFloeCylindrical!(sim, [.1, .1],
       +                                                          10., 1., 
       +                                                          lin_acc=[.2,.2,.2])
       +@test_throws ErrorException SeaIce.addIceFloeCylindrical!(sim, [.1, .1],
       +                                                          0., 1.)
       +@test_throws ErrorException SeaIce.addIceFloeCylindrical!(sim, [.1, .1],
       +                                                          10., 1., density=-2.)
        @test_throws ErrorException SeaIce.disableIceFloe!(sim, 0)
        
        sim = SeaIce.createSimulation(id="test")
       -SeaIce.addIceFloeCylindrical(sim, [ 0., 0.], 10., 1., verbose=false)
       -SeaIce.addIceFloeCylindrical(sim, [ 0., 0.], 10., 1., verbose=false)
       +SeaIce.addIceFloeCylindrical!(sim, [ 0., 0.], 10., 1., verbose=false)
       +SeaIce.addIceFloeCylindrical!(sim, [ 0., 0.], 10., 1., verbose=false)
        SeaIce.compareIceFloes(sim.ice_floes[1], sim.ice_floes[2])
 (DIR) diff --git a/test/jld.jl b/test/jld.jl
       t@@ -8,8 +8,8 @@ if typeof(Pkg.installed("JLD")) == VersionNumber
        
            info("Writing simple simulation to JLD file")
            sim = SeaIce.createSimulation(id="test")
       -    SeaIce.addIceFloeCylindrical(sim, [ 0., 0.], 10., 1., verbose=false)
       -    SeaIce.addIceFloeCylindrical(sim, [18., 0.], 10., 1., verbose=false)
       +    SeaIce.addIceFloeCylindrical!(sim, [ 0., 0.], 10., 1., verbose=false)
       +    SeaIce.addIceFloeCylindrical!(sim, [18., 0.], 10., 1., verbose=false)
            sim.ocean = SeaIce.createRegularOceanGrid([10, 20, 5], [10., 25., 2.])  
            SeaIce.findContacts!(sim, method="all to all")
            SeaIce.writeVTK(sim, verbose=false)
 (DIR) diff --git a/test/ocean.jl b/test/ocean.jl
       t@@ -28,7 +28,7 @@ sim.ocean = SeaIce.createRegularOceanGrid([6, 6, 6], [1., 1., 1.])
        @test sim.ocean.e ≈ zeros(7, 7, 6, 1)
        
        info("Testing velocity drag interaction (static ocean)")
       -SeaIce.addIceFloeCylindrical(sim, [.5, .5], .25, .1)
       +SeaIce.addIceFloeCylindrical!(sim, [.5, .5], .25, .1)
        SeaIce.setTotalTime!(sim, 5.)
        SeaIce.setTimeStep!(sim)
        sim_init = deepcopy(sim)
 (DIR) diff --git a/test/profiling.jl b/test/profiling.jl
       t@@ -39,8 +39,8 @@ function timeSingleStepInDenseSimulation(nx::Int; verbose::Bool=true,
                    if ix == 1 || iy == 1 || ix == nx || iy == ny
                        fixed = true
                    end
       -            SeaIce.addIceFloeCylindrical(sim, [x, y], r*1.1, 1.,
       -                                         fixed=fixed, verbose=false)
       +            SeaIce.addIceFloeCylindrical!(sim, [x, y], r*1.1, 1.,
       +                                          fixed=fixed, verbose=false)
                end
            end
            print_with_color(:green, "number of ice floes: $(length(sim.ice_floes))\n")
 (DIR) diff --git a/test/temporal.jl b/test/temporal.jl
       t@@ -12,7 +12,7 @@ SeaIce.setCurrentTime!(sim, 1.)
        SeaIce.setOutputFileInterval!(sim, 0.)
        SeaIce.disableOutputFiles!(sim)
        @test_throws ErrorException SeaIce.checkTimeParameters(sim)
       -SeaIce.addIceFloeCylindrical(sim, [.1,.1], 2., 2.)
       +SeaIce.addIceFloeCylindrical!(sim, [.1,.1], 2., 2.)
        sim.ice_floes[1].mass = 0.
        @test_throws ErrorException SeaIce.setTimeStep!(sim)
        
 (DIR) diff --git a/test/vtk.jl b/test/vtk.jl
       t@@ -6,8 +6,8 @@ info("#### $(basename(@__FILE__)) ####")
        
        info("Writing simple simulation to VTK file")
        sim = SeaIce.createSimulation(id="test")
       -SeaIce.addIceFloeCylindrical(sim, [ 0., 0.], 10., 1., verbose=false)
       -SeaIce.addIceFloeCylindrical(sim, [18., 0.], 10., 1., verbose=false)
       +SeaIce.addIceFloeCylindrical!(sim, [ 0., 0.], 10., 1., verbose=false)
       +SeaIce.addIceFloeCylindrical!(sim, [18., 0.], 10., 1., verbose=false)
        sim.ocean = SeaIce.createRegularOceanGrid([10, 20, 5], [10., 25., 2.])  
        SeaIce.findContacts!(sim, method="all to all")
        SeaIce.writeVTK(sim, verbose=false)