commit 76c46ba160aa80b209aa7e20e302b170c87ef454
parent 2379861ee82b96220b667e149c44ed11172f7dc4
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date: Mon, 9 Jul 2018 11:31:33 -0400
Merge pull request #6 from anders-dc/ridging
Implement compressive failure parameterization
Diffstat:
32 files changed, 885 insertions(+), 520 deletions(-)
diff --git a/src/atmosphere.jl b/src/atmosphere.jl
@@ -169,10 +169,10 @@ function applyAtmosphereDragToGrain!(grain::GrainCylindrical,
drag_force = ρ_a * π *
(2.0*grain.ocean_drag_coeff_vert*grain.areal_radius*.1*grain.thickness +
grain.atmosphere_drag_coeff_horiz*grain.areal_radius^2.0) *
- ([u, v] - grain.lin_vel)*norm([u, v] - grain.lin_vel)
+ ([u, v] - grain.lin_vel[1:2])*norm([u, v] - grain.lin_vel[1:2])
- grain.force += drag_force
- grain.atmosphere_stress = drag_force/grain.horizontal_surface_area
+ grain.force += vecTo3d(drag_force)
+ grain.atmosphere_stress = vecTo3d(drag_force/grain.horizontal_surface_area)
nothing
end
@@ -186,12 +186,12 @@ function applyAtmosphereVorticityToGrain!(grain::GrainCylindrical,
atmosphere_curl::Float64)
ρ_a = 1.2754 # atmosphere density
- grain.torque +=
+ grain.torque[3] +=
π * grain.areal_radius^4. * ρ_a *
(grain.areal_radius / 5. * grain.atmosphere_drag_coeff_horiz +
.1 * grain.thickness * grain.atmosphere_drag_coeff_vert) *
- abs(.5 * atmosphere_curl - grain.ang_vel) *
- (.5 * atmosphere_curl - grain.ang_vel)
+ norm(.5 * atmosphere_curl - grain.ang_vel[3]) *
+ (.5 * atmosphere_curl - grain.ang_vel[3])
nothing
end
diff --git a/src/contact_search.jl b/src/contact_search.jl
@@ -109,7 +109,7 @@ Perform an O(n*log(n)) cell-based contact search between all grains in the
"""
function findContactsInGrid!(simulation::Simulation, grid::Any)
- distance_modifier = [0., 0.]
+ distance_modifier = [0., 0., 0.]
i_corrected = 0
j_corrected = 0
@@ -173,7 +173,7 @@ function checkForContacts(simulation::Simulation,
r_candidate::Float64;
return_when_overlap_found::Bool=false)
- distance_modifier = zeros(2)
+ distance_modifier = zeros(3)
overlaps_found::Integer = 0
# Inter-grain position vector and grain overlap
@@ -225,7 +225,7 @@ function periodicBoundaryCorrection!(grid::Any, i::Integer, j::Integer)
# vector for correcting inter-particle distance in case of
# boundary periodicity
- distance_modifier = zeros(2)
+ distance_modifier = zeros(3)
# i and j are not corrected for periodic boundaries
i_corrected = i
@@ -275,7 +275,7 @@ written to `simulation.contact_parallel_displacement`.
boundaries.
"""
function checkAndAddContact!(sim::Simulation, i::Int, j::Int,
- distance_modifier::Vector{Float64} = [0., 0.])
+ distance_modifier::Vector{Float64} = [0., 0., 0.])
if i < j
@inbounds if !sim.grains[i].enabled || !sim.grains[j].enabled
@@ -327,8 +327,10 @@ function checkAndAddContact!(sim::Simulation, i::Int, j::Int,
position_ij
@inbounds fill!(sim.grains[i].
contact_parallel_displacement[ic] , 0.)
- @inbounds sim.grains[i].contact_rotation[ic] = 0.
+ @inbounds fill!(sim.grains[i].contact_rotation[ic] , 0.)
@inbounds sim.grains[i].contact_age[ic] = 0.
+ @inbounds sim.grains[i].contact_area[ic] = 0.
+ @inbounds sim.grains[i].compressive_failure[ic] = false
break
end
end
diff --git a/src/datatypes.jl b/src/datatypes.jl
@@ -24,15 +24,16 @@ mutable struct GrainCylindrical
lin_disp::Vector{Float64}
# Angular kinematic degrees of freedom for vertical rotation around center
- ang_pos::Float64
- ang_vel::Float64
- ang_acc::Float64
- torque::Float64
+ ang_pos::Vector{Float64}
+ ang_vel::Vector{Float64}
+ ang_acc::Vector{Float64}
+ torque::Vector{Float64}
# Kinematic constraint flags
fixed::Bool
allow_x_acc::Bool
allow_y_acc::Bool
+ allow_z_acc::Bool
rotating::Bool
enabled::Bool
@@ -65,8 +66,10 @@ mutable struct GrainCylindrical
contacts::Vector{Int}
position_vector::Vector{Vector{Float64}}
contact_parallel_displacement::Vector{Vector{Float64}}
- contact_rotation::Vector{Float64}
+ contact_rotation::Vector{Vector{Float64}}
contact_age::Vector{Float64}
+ contact_area::Vector{Float64}
+ compressive_failure::Vector{Bool}
granular_stress::Vector{Float64}
ocean_stress::Vector{Float64}
@@ -128,6 +131,7 @@ mutable struct GrainArrays
fixed::Vector{Int}
allow_x_acc::Vector{Int}
allow_y_acc::Vector{Int}
+ allow_z_acc::Vector{Int}
rotating::Vector{Int}
enabled::Vector{Int}
diff --git a/src/grain.jl b/src/grain.jl
@@ -24,7 +24,7 @@ export addGrainCylindrical!
atmosphere_drag_coeff_vert,
atmosphere_drag_coeff_horiz,
pressure, fixed,
- allow_x_acc, allow_y_acc,
+ allow_x_acc, allow_y_acc, allow_z_acc,
rotating, enabled, verbose,
ocean_grid_pos, atmosphere_grid_pos,
n_contact, granular_stress, ocean_stress,
@@ -39,21 +39,34 @@ are optional, and come with default values. The only required arguments are
# Arguments
* `simulation::Simulation`: the simulation object where the grain should be
added to.
-* `lin_pos::Vector{Float64}`: linear position of grain center [m].
+* `lin_pos::Vector{Float64}`: linear position of grain center [m]. If a
+ two-component vector is used, the values will be mapped to *x* and *y*, and
+ the *z* component will be set to zero.
* `contact_radius::Float64`: grain radius for granular interaction [m].
* `thickness::Float64`: grain thickness [m].
* `areal_radius = false`: grain radius for determining sea-ice concentration
[m].
-* `lin_vel::Vector{Float64} = [0., 0.]`: linear velocity [m/s].
-* `lin_acc::Vector{Float64} = [0., 0.]`: linear acceleration [m/s^2].
-* `force::Vector{Float64} = [0., 0.]`: linear force balance [N].
-* `ang_pos::Float64 = 0.`: angular position around its center vertical axis
- [rad].
-* `ang_vel::Float64 = 0.`: angular velocity around its center vertical axis
- [rad/s].
-* `ang_acc::Float64 = 0.`: angular acceleration around its center vertical axis
- [rad/s^2].
-* `torque::Float64 = 0.`: torque around its center vertical axis [N*m]
+* `lin_vel::Vector{Float64} = [0., 0., 0.]`: linear velocity [m/s]. If a
+ two-component vector is used, the values will be mapped to *x* and *y*, and
+ the *z* component will be set to zero.
+* `lin_acc::Vector{Float64} = [0., 0., 0.]`: linear acceleration [m/s^2]. If a
+ two-component vector is used, the values will be mapped to *x* and *y*, and
+ the *z* component will be set to zero.
+* `force::Vector{Float64} = [0., 0., 0.]`: linear force balance [N]. If a
+ two-component vector is used, the values will be mapped to *x* and *y*, and
+ the *z* component will be set to zero.
+* `ang_pos::Float64 = [0., 0., 0.]`: angular position around its center vertical
+ axis [rad]. If a scalar is used, the value will be mapped to *z*, and the
+ *x* and *y* components will be set to zero.
+* `ang_vel::Float64 = [0., 0., 0.]`: angular velocity around its center vertical
+ axis [rad/s]. If a scalar is used, the value will be mapped to *z*, and the
+ *x* and *y* components will be set to zero.
+* `ang_acc::Float64 = [0., 0., 0.]`: angular acceleration around its center
+ vertical axis [rad/s^2]. If a scalar is used, the value will be mapped to
+ *z*, and the *x* and *y* components will be set to zero.
+* `torque::Float64 = [0., 0., 0.]`: torque around its center vertical axis
+ [N*m]. If a scalar is used, the value will be mapped to *z*, and the *x* and
+ *y* components will be set to zero.
* `density::Float64 = 934.`: grain mean density [kg/m^3].
* `contact_stiffness_normal::Float64 = 1e7`: contact-normal stiffness [N/m];
overridden if `youngs_modulus` is set to a positive value.
@@ -75,8 +88,8 @@ are optional, and come with default values. The only required arguments are
* `shear_strength::Float64 = 0.`: shear strength of bonded contacts [Pa].
* `strength_heal_rate::Float64 = 0.`: rate at which contact bond
strength is obtained [Pa/s].
-* `fracture_toughness::Float64 = 0.`: maximum compressive
- strength on granular contact (not currently enforced) [m^{1/2}*Pa]. A value
+* `fracture_toughness::Float64 = 0.`: fracture toughness which influences the
+ maximum compressive strength on granular contact [m^{1/2}*Pa]. A value
of 1.285e3 m^{1/2}*Pa is used for sea ice by Hopkins 2004.
* `ocean_drag_coeff_vert::Float64 = 0.85`: vertical drag coefficient for ocean
against grain sides [-].
@@ -90,6 +103,7 @@ are optional, and come with default values. The only required arguments are
* `fixed::Bool = false`: grain is fixed to a constant velocity (e.g. zero).
* `allow_x_acc::Bool = false`: override `fixed` along `x`.
* `allow_y_acc::Bool = false`: override `fixed` along `y`.
+* `allow_z_acc::Bool = false`: override `fixed` along `z`.
* `rotating::Bool = true`: grain is allowed to rotate.
* `enabled::Bool = true`: grain interacts with other grains.
* `verbose::Bool = true`: display diagnostic information during the function
@@ -99,11 +113,11 @@ are optional, and come with default values. The only required arguments are
* `atmosphere_grid_pos::Array{Int, 1} = [0, 0]`: position of grain in the
atmosphere grid.
* `n_contacts::Int = 0`: number of contacts with other grains.
-* `granular_stress::Vector{Float64} = [0., 0.]`: resultant stress on grain
+* `granular_stress::Vector{Float64} = [0., 0., 0.]`: resultant stress on grain
from granular interactions [Pa].
-* `ocean_stress::Vector{Float64} = [0., 0.]`: resultant stress on grain from
+* `ocean_stress::Vector{Float64} = [0., 0., 0.]`: resultant stress on grain from
ocean drag [Pa].
-* `atmosphere_stress::Vector{Float64} = [0., 0.]`: resultant stress on grain
+* `atmosphere_stress::Vector{Float64} = [0., 0., 0.]`: resultant stress on grain
from atmosphere drag [Pa].
* `thermal_energy::Float64 = 0.0`: thermal energy of grain [J].
* `color::Int=0`: type number, usually used for associating a color to the grain
@@ -111,12 +125,13 @@ are optional, and come with default values. The only required arguments are
# Examples
The most basic example adds a new grain to the simulation `sim`, with a
-center at `[1., 2.]`, a radius of `1.` meter, and a thickness of `0.5`
+center at `[1., 2., 0.]`, a radius of `1.` meter, and a thickness of `0.5`
meter:
```julia
Granular.addGrainCylindrical!(sim, [1., 2.], 1., .5)
```
+Note that the *z* component is set to zero if a two-component vector is passed.
The following example will create a grain with tensile and shear strength, and a
velocity of 0.5 m/s towards -x:
@@ -142,13 +157,13 @@ function addGrainCylindrical!(simulation::Simulation,
contact_radius::Float64,
thickness::Float64;
areal_radius = false,
- lin_vel::Vector{Float64} = [0., 0.],
- lin_acc::Vector{Float64} = [0., 0.],
- force::Vector{Float64} = [0., 0.],
- ang_pos::Float64 = 0.,
- ang_vel::Float64 = 0.,
- ang_acc::Float64 = 0.,
- torque::Float64 = 0.,
+ lin_vel::Vector{Float64} = [0., 0., 0.],
+ lin_acc::Vector{Float64} = [0., 0., 0.],
+ force::Vector{Float64} = [0., 0., 0.],
+ ang_pos::Vector{Float64} = [0., 0., 0.],
+ ang_vel::Vector{Float64} = [0., 0., 0.],
+ ang_acc::Vector{Float64} = [0., 0., 0.],
+ torque::Vector{Float64} = [0., 0., 0.],
density::Float64 = 934.,
contact_stiffness_normal::Float64 = 1e7,
contact_stiffness_tangential::Float64 = 0.,
@@ -170,30 +185,52 @@ function addGrainCylindrical!(simulation::Simulation,
fixed::Bool = false,
allow_x_acc::Bool = false,
allow_y_acc::Bool = false,
+ allow_z_acc::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{Float64} = [0., 0.],
- ocean_stress::Vector{Float64} = [0., 0.],
- atmosphere_stress::Vector{Float64} = [0., 0.],
+ granular_stress::Vector{Float64} = [0., 0., 0.],
+ ocean_stress::Vector{Float64} = [0., 0., 0.],
+ atmosphere_stress::Vector{Float64} = [0., 0., 0.],
thermal_energy::Float64 = 0.,
color::Int = 0)
# Check input values
- if length(lin_pos) != 2
- error("Linear position must be a two-element array " *
- "(lin_pos = $lin_pos)")
+ if length(lin_pos) != 3
+ lin_pos = vecTo3d(lin_pos)
end
- if length(lin_vel) != 2
- error("Linear velocity must be a two-element array " *
- "(lin_vel = $lin_vel)")
+ if length(lin_vel) != 3
+ lin_vel = vecTo3d(lin_vel)
end
- if length(lin_acc) != 2
- error("Linear acceleration must be a two-element array " *
- "(lin_acc = $lin_acc)")
+ if length(lin_acc) != 3
+ lin_acc = vecTo3d(lin_acc)
+ end
+ if length(force) != 3
+ force = vecTo3d(force)
+ end
+ if length(ang_pos) != 3
+ ang_pos = vecTo3d(ang_pos)
+ end
+ if length(ang_vel) != 3
+ ang_vel = vecTo3d(ang_vel)
+ end
+ if length(ang_acc) != 3
+ ang_acc = vecTo3d(ang_acc)
+ end
+ if length(torque) != 3
+ torque = vecTo3d(torque)
+ end
+ if length(granular_stress) != 3
+ granular_stress = vecTo3d(granular_stress)
+ end
+ if length(ocean_stress) != 3
+ ocean_stress = vecTo3d(ocean_stress)
+ end
+ if length(atmosphere_stress) != 3
+ atmosphere_stress = vecTo3d(atmosphere_stress)
end
if contact_radius <= 0.0
error("Radius must be greater than 0.0 " *
@@ -212,81 +249,87 @@ function addGrainCylindrical!(simulation::Simulation,
position_vector = Vector{Vector{Float64}}(undef, simulation.Nc_max)
contact_parallel_displacement =
Vector{Vector{Float64}}(undef, simulation.Nc_max)
- contact_rotation::Vector{Float64} = zeros(Float64, simulation.Nc_max)
+ contact_rotation = Vector{Vector{Float64}}(undef, simulation.Nc_max)
contact_age::Vector{Float64} = zeros(Float64, simulation.Nc_max)
+ contact_area::Vector{Float64} = zeros(Float64, simulation.Nc_max)
+ compressive_failure::Vector{Bool} = zeros(Bool, simulation.Nc_max)
for i=1:simulation.Nc_max
- position_vector[i] = zeros(2)
- contact_parallel_displacement[i] = zeros(2)
+ position_vector[i] = zeros(3)
+ contact_rotation[i] = zeros(3)
+ contact_parallel_displacement[i] = zeros(3)
end
# Create grain object with placeholder values for surface area, volume,
# mass, and moment of inertia.
grain = GrainCylindrical(density,
- thickness,
- contact_radius,
- areal_radius,
- 1.0, # circumreference
- 1.0, # horizontal_surface_area
- 1.0, # side_surface_area
- 1.0, # volume
- 1.0, # mass
- 1.0, # moment_of_inertia
- lin_pos,
- lin_vel,
- lin_acc,
- force,
- [0., 0.], # external_body_force
- [0., 0.], # lin_disp
-
- ang_pos,
- ang_vel,
- ang_acc,
- torque,
-
- fixed,
- allow_x_acc,
- allow_y_acc,
- rotating,
- enabled,
-
- contact_stiffness_normal,
- contact_stiffness_tangential,
- contact_viscosity_normal,
- contact_viscosity_tangential,
- contact_static_friction,
- contact_dynamic_friction,
-
- youngs_modulus,
- poissons_ratio,
- tensile_strength,
- shear_strength,
- strength_heal_rate,
- fracture_toughness,
-
- ocean_drag_coeff_vert,
- ocean_drag_coeff_horiz,
- atmosphere_drag_coeff_vert,
- atmosphere_drag_coeff_horiz,
-
- pressure,
- n_contacts,
- ocean_grid_pos,
- atmosphere_grid_pos,
- contacts,
- position_vector,
- contact_parallel_displacement,
- contact_rotation,
- contact_age,
-
- granular_stress,
- ocean_stress,
- atmosphere_stress,
-
- thermal_energy,
-
- color
- )
+ thickness,
+ contact_radius,
+ areal_radius,
+ 1.0, # circumreference
+ 1.0, # horizontal_surface_area
+ 1.0, # side_surface_area
+ 1.0, # volume
+ 1.0, # mass
+ 1.0, # moment_of_inertia
+ lin_pos,
+ lin_vel,
+ lin_acc,
+ force,
+ [0., 0., 0.], # external_body_force
+ [0., 0., 0.], # lin_disp
+
+ ang_pos,
+ ang_vel,
+ ang_acc,
+ torque,
+
+ fixed,
+ allow_x_acc,
+ allow_y_acc,
+ allow_z_acc,
+ rotating,
+ enabled,
+
+ contact_stiffness_normal,
+ contact_stiffness_tangential,
+ contact_viscosity_normal,
+ contact_viscosity_tangential,
+ contact_static_friction,
+ contact_dynamic_friction,
+
+ youngs_modulus,
+ poissons_ratio,
+ tensile_strength,
+ shear_strength,
+ strength_heal_rate,
+ fracture_toughness,
+
+ ocean_drag_coeff_vert,
+ ocean_drag_coeff_horiz,
+ atmosphere_drag_coeff_vert,
+ atmosphere_drag_coeff_horiz,
+
+ pressure,
+ n_contacts,
+ ocean_grid_pos,
+ atmosphere_grid_pos,
+ contacts,
+ position_vector,
+ contact_parallel_displacement,
+ contact_rotation,
+ contact_age,
+ contact_area,
+ compressive_failure,
+
+ granular_stress,
+ ocean_stress,
+ atmosphere_stress,
+
+ thermal_energy,
+
+ color
+ )
# Overwrite previous placeholder values
grain.circumreference = grainCircumreference(grain)
@@ -400,6 +443,8 @@ function convertGrainDataToArrays(simulation::Simulation)
Array{Int}(undef, length(simulation.grains)),
## allow_y_acc
Array{Int}(undef, length(simulation.grains)),
+ ## allow_z_acc
+ Array{Int}(undef, length(simulation.grains)),
## rotating
Array{Int}(undef, length(simulation.grains)),
## enabled
@@ -479,22 +524,23 @@ function convertGrainDataToArrays(simulation::Simulation)
ifarr.mass[i] = simulation.grains[i].mass
ifarr.moment_of_inertia[i] = simulation.grains[i].moment_of_inertia
- ifarr.lin_pos[1:2, i] = simulation.grains[i].lin_pos
- ifarr.lin_vel[1:2, i] = simulation.grains[i].lin_vel
- ifarr.lin_acc[1:2, i] = simulation.grains[i].lin_acc
- ifarr.force[1:2, i] = simulation.grains[i].force
- ifarr.external_body_force[1:2, i] =
+ ifarr.lin_pos[1:3, i] = simulation.grains[i].lin_pos
+ ifarr.lin_vel[1:3, i] = simulation.grains[i].lin_vel
+ ifarr.lin_acc[1:3, i] = simulation.grains[i].lin_acc
+ ifarr.force[1:3, i] = simulation.grains[i].force
+ ifarr.external_body_force[1:3, i] =
simulation.grains[i].external_body_force
- ifarr.lin_disp[1:2, i] = simulation.grains[i].lin_disp
+ ifarr.lin_disp[1:3, i] = simulation.grains[i].lin_disp
- ifarr.ang_pos[3, i] = simulation.grains[i].ang_pos
- ifarr.ang_vel[3, i] = simulation.grains[i].ang_vel
- ifarr.ang_acc[3, i] = simulation.grains[i].ang_acc
- ifarr.torque[3, i] = simulation.grains[i].torque
+ ifarr.ang_pos[1:3, i] = simulation.grains[i].ang_pos
+ ifarr.ang_vel[1:3, i] = simulation.grains[i].ang_vel
+ ifarr.ang_acc[1:3, i] = simulation.grains[i].ang_acc
+ ifarr.torque[1:3, i] = simulation.grains[i].torque
ifarr.fixed[i] = Int(simulation.grains[i].fixed)
ifarr.allow_x_acc[i] = Int(simulation.grains[i].allow_x_acc)
ifarr.allow_y_acc[i] = Int(simulation.grains[i].allow_y_acc)
+ ifarr.allow_z_acc[i] = Int(simulation.grains[i].allow_z_acc)
ifarr.rotating[i] = Int(simulation.grains[i].rotating)
ifarr.enabled[i] = Int(simulation.grains[i].enabled)
@@ -531,10 +577,9 @@ function convertGrainDataToArrays(simulation::Simulation)
ifarr.pressure[i] = simulation.grains[i].pressure
ifarr.n_contacts[i] = simulation.grains[i].n_contacts
- ifarr.granular_stress[1:2, i] = simulation.grains[i].granular_stress
- ifarr.ocean_stress[1:2, i] = simulation.grains[i].ocean_stress
- ifarr.atmosphere_stress[1:2, i] =
- simulation.grains[i].atmosphere_stress
+ ifarr.granular_stress[1:3, i] = simulation.grains[i].granular_stress
+ ifarr.ocean_stress[1:3, i] = simulation.grains[i].ocean_stress
+ ifarr.atmosphere_stress[1:3, i] = simulation.grains[i].atmosphere_stress
ifarr.thermal_energy[i] = simulation.grains[i].thermal_energy
@@ -576,6 +621,7 @@ function deleteGrainArrays!(ifarr::GrainArrays)
ifarr.fixed = i1
ifarr.allow_x_acc = i1
ifarr.allow_y_acc = i1
+ ifarr.allow_z_acc = i1
ifarr.rotating = i1
ifarr.enabled = i1
@@ -643,22 +689,23 @@ function printGrainInfo(f::GrainCylindrical)
println(" fixed: $(f.fixed)")
println(" allow_x_acc: $(f.allow_x_acc)")
println(" allow_y_acc: $(f.allow_y_acc)")
+ println(" allow_z_acc: $(f.allow_z_acc)")
println(" rotating: $(f.rotating)")
println(" enabled: $(f.enabled)\n")
- println(" k_n: $(f.contact_stiffness_normal) N/m")
- println(" k_t: $(f.contact_stiffness_tangential) N/m")
+ println(" k_n: $(f.contact_stiffness_normal) N/m")
+ println(" k_t: $(f.contact_stiffness_tangential) N/m")
println(" γ_n: $(f.contact_viscosity_normal) N/(m/s)")
println(" γ_t: $(f.contact_viscosity_tangential) N/(m/s)")
- println(" μ_s: $(f.contact_static_friction)")
- println(" μ_d: $(f.contact_dynamic_friction)\n")
+ println(" μ_s: $(f.contact_static_friction)")
+ println(" μ_d: $(f.contact_dynamic_friction)\n")
- println(" E: $(f.youngs_modulus) Pa")
- println(" ν: $(f.poissons_ratio)")
- println(" tensile_strength: $(f.tensile_strength) Pa")
- println(" shear_strength: $(f.shear_strength) Pa")
- println(" strength_heal_rate: $(f.strength_heal_rate) Pa/s")
- println(" fracture_toughness: $(f.fracture_toughness) m^0.5 Pa\n")
+ println(" E: $(f.youngs_modulus) Pa")
+ println(" ν: $(f.poissons_ratio)")
+ println(" tensile_strength: $(f.tensile_strength) Pa")
+ println(" shear_strength: $(f.shear_strength) Pa")
+ println(" strength_heal_rate: $(f.strength_heal_rate) Pa/s")
+ println(" fracture_toughness: $(f.fracture_toughness) m^0.5 Pa\n")
println(" c_o_v: $(f.ocean_drag_coeff_vert)")
println(" c_o_h: $(f.ocean_drag_coeff_horiz)")
@@ -754,7 +801,7 @@ Add to the value of the external body force on a grain.
* `force::Vector{Float64}`: a vector of force [N]
"""
function addBodyForce!(grain::GrainCylindrical, force::Vector{Float64})
- grain.external_body_force += force
+ grain.external_body_force += vecTo3d(force)
end
export setBodyForce!
@@ -768,7 +815,7 @@ Set the value of the external body force on a grain.
* `force::Vector{Float64}`: a vector of force [N]
"""
function setBodyForce!(grain::GrainCylindrical, force::Vector{Float64})
- grain.external_body_force = force
+ grain.external_body_force = vecTo3d(force)
end
export compareGrains
@@ -781,11 +828,9 @@ function compareGrains(if1::GrainCylindrical, if2::GrainCylindrical)
@test if1.density ≈ if2.density
@test if1.thickness ≈ if2.thickness
- @test if1.contact_radius ≈
- if2.contact_radius
+ @test if1.contact_radius ≈ if2.contact_radius
@test if1.areal_radius ≈ if2.areal_radius
- @test if1.circumreference ≈
- if2.circumreference
+ @test if1.circumreference ≈ if2.circumreference
@test if1.horizontal_surface_area ≈ if2.horizontal_surface_area
@test if1.side_surface_area ≈ if2.side_surface_area
@test if1.volume ≈ if2.volume
@@ -809,11 +854,9 @@ function compareGrains(if1::GrainCylindrical, if2::GrainCylindrical)
@test if1.enabled == if2.enabled
@test if1.contact_stiffness_normal ≈ if2.contact_stiffness_normal
- @test if1.contact_stiffness_tangential ≈
- if2.contact_stiffness_tangential
+ @test if1.contact_stiffness_tangential ≈ if2.contact_stiffness_tangential
@test if1.contact_viscosity_normal ≈ if2.contact_viscosity_normal
- @test if1.contact_viscosity_tangential ≈
- if2.contact_viscosity_tangential
+ @test if1.contact_viscosity_tangential ≈ if2.contact_viscosity_tangential
@test if1.contact_static_friction ≈ if2.contact_static_friction
@test if1.contact_dynamic_friction ≈ if2.contact_dynamic_friction
@@ -822,25 +865,23 @@ function compareGrains(if1::GrainCylindrical, if2::GrainCylindrical)
@test if1.tensile_strength ≈ if2.tensile_strength
@test if1.shear_strength ≈ if2.shear_strength
@test if1.strength_heal_rate ≈ if2.strength_heal_rate
- @test if1.fracture_toughness ≈
- if2.fracture_toughness
+ @test if1.fracture_toughness ≈ if2.fracture_toughness
@test if1.ocean_drag_coeff_vert ≈ if2.ocean_drag_coeff_vert
@test if1.ocean_drag_coeff_horiz ≈ if2.ocean_drag_coeff_horiz
- @test if1.atmosphere_drag_coeff_vert ≈
- if2.atmosphere_drag_coeff_vert
- @test if1.atmosphere_drag_coeff_horiz ≈
- if2.atmosphere_drag_coeff_horiz
+ @test if1.atmosphere_drag_coeff_vert ≈ if2.atmosphere_drag_coeff_vert
+ @test if1.atmosphere_drag_coeff_horiz ≈ if2.atmosphere_drag_coeff_horiz
@test if1.pressure ≈ if2.pressure
@test if1.n_contacts == if2.n_contacts
@test if1.ocean_grid_pos == if2.ocean_grid_pos
@test if1.contacts == if2.contacts
@test if1.position_vector == if2.position_vector
- @test if1.contact_parallel_displacement ==
- if2.contact_parallel_displacement
+ @test if1.contact_parallel_displacement == if2.contact_parallel_displacement
@test if1.contact_rotation == if2.contact_rotation
@test if1.contact_age ≈ if2.contact_age
+ @test if1.contact_area ≈ if2.contact_area
+ @test if1.compressive_failure ≈ if2.compressive_failure
@test if1.granular_stress ≈ if2.granular_stress
@test if1.ocean_stress ≈ if2.ocean_stress
@@ -905,12 +946,12 @@ rotational) to zero in order to get rid of all kinetic energy.
"""
function zeroKinematics!(sim::Simulation)
for grain in sim.grains
- grain.lin_vel .= zeros(2)
- grain.lin_acc .= zeros(2)
- grain.force .= zeros(2)
- grain.lin_disp .= zeros(2)
- grain.ang_vel = 0.
- grain.ang_acc = 0.
- grain.torque = 0.
+ grain.lin_vel .= zeros(3)
+ grain.lin_acc .= zeros(3)
+ grain.force .= zeros(3)
+ grain.lin_disp .= zeros(3)
+ grain.ang_vel .= zeros(3)
+ grain.ang_acc .= zeros(3)
+ grain.torque .= zeros(3)
end
end
diff --git a/src/grid.jl b/src/grid.jl
@@ -161,14 +161,15 @@ function sortGrainsInGrid!(simulation::Simulation, grid::Any; verbose=true)
!grid.regular_grid &&
i_old > 0 && j_old > 0 &&
isPointInCell(grid, i_old, j_old,
- simulation.grains[idx].lin_pos, sw, se, ne, nw)
+ simulation.grains[idx].lin_pos[1:2], sw, se, ne, nw)
i = i_old
j = j_old
else
if grid.regular_grid
- i, j = Int.(floor.((simulation.grains[idx].lin_pos - grid.origo)
+ i, j = Int.(floor.((simulation.grains[idx].lin_pos[1:2]
+ - grid.origo)
./ grid.dx[1:2])) + [1,1]
else
@@ -185,8 +186,9 @@ function sortGrainsInGrid!(simulation::Simulation, grid::Any; verbose=true)
j_t = max(min(j_old + j_rel, ny), 1)
@inbounds if isPointInCell(grid, i_t, j_t,
- simulation.grains[idx].lin_pos,
- sw, se, ne, nw)
+ simulation.grains[idx].
+ lin_pos[1:2],
+ sw, se, ne, nw)
i = i_t
j = j_t
found = true
@@ -201,7 +203,7 @@ function sortGrainsInGrid!(simulation::Simulation, grid::Any; verbose=true)
if !found
i, j = findCellContainingPoint(grid,
simulation.grains[idx].
- lin_pos,
+ lin_pos[1:2],
sw, se, ne, nw)
end
end
@@ -285,10 +287,10 @@ This function is a wrapper for `getCellCornerCoordinates()` and
function getNonDimensionalCellCoordinates(grid::Any, i::Int, j::Int,
point::Vector{Float64})
if grid.regular_grid
- return (point - Float64.([i-1,j-1]).*grid.dx[1:2])./grid.dx[1:2]
+ return (point[1:2] - Float64.([i-1,j-1]).*grid.dx[1:2])./grid.dx[1:2]
else
sw, se, ne, nw = getCellCornerCoordinates(grid.xq, grid.yq, i, j)
- return conformalQuadrilateralCoordinates(sw, se, ne, nw, point)
+ return conformalQuadrilateralCoordinates(sw, se, ne, nw, point[1:2])
end
end
@@ -308,7 +310,7 @@ function isPointInCell(grid::Any, i::Int, j::Int,
method::String="Conformal")
if grid.regular_grid
- if [i,j] == Int.(floor.((point - grid.origo) ./ grid.dx[1:2])) + [1,1]
+ if [i,j] == Int.(floor.((point[1:2] - grid.origo) ./ grid.dx[1:2])) + [1,1]
return true
else
return false
@@ -620,8 +622,8 @@ function findEmptyPositionInGridCell(simulation::Simulation,
# traverse list of grains in the target cell and check
# for overlaps
for grain_idx in grid.grain_list[it, jt]
- overlap = norm(simulation.grains[grain_idx].lin_pos -
- pos) -
+ overlap = norm(simulation.grains[grain_idx].
+ lin_pos[1:2] - pos) -
(simulation.grains[grain_idx].contact_radius + r)
if overlap < 0.
diff --git a/src/interaction.jl b/src/interaction.jl
@@ -63,10 +63,11 @@ function interactWalls!(sim::Simulation)
# get overlap distance by projecting grain position onto wall-normal
# vector. Overlap when δ_n < 0.0
- δ_n = abs(dot(sim.walls[iw].normal, sim.grains[i].lin_pos) -
- sim.walls[iw].pos) - sim.grains[i].contact_radius
+ δ_n = norm(dot(sim.walls[iw].normal[1:2],
+ sim.grains[i].lin_pos[1:2]) -
+ sim.walls[iw].pos) - sim.grains[i].contact_radius
- vel_n = dot(sim.walls[iw].normal, sim.grains[i].lin_vel)
+ vel_n = dot(sim.walls[iw].normal[1:2], sim.grains[i].lin_vel[1:2])
if δ_n < 0.
if sim.grains[i].youngs_modulus > 0.
@@ -98,7 +99,7 @@ function interactWalls!(sim::Simulation)
end
sim.grains[i].force += -force_n .* sim.walls[iw].normal .*
- orientation
+ orientation
sim.walls[iw].force += force_n * orientation
end
end
@@ -127,7 +128,7 @@ function interactGrains!(simulation::Simulation, i::Int, j::Int, ic::Int)
# Inter-position vector, use stored value which is corrected for boundary
# periodicity
- p = simulation.grains[i].position_vector[ic]
+ p = simulation.grains[i].position_vector[ic][1:2]
dist = norm(p)
r_i = simulation.grains[i].contact_radius
@@ -137,11 +138,12 @@ function interactGrains!(simulation::Simulation, i::Int, j::Int, ic::Int)
δ_n = dist - (r_i + r_j)
# Local axes
- n = p/dist
+ n = p/max(dist, 1e-12)
t = [-n[2], n[1]]
- # Contact kinematics
- vel_lin = simulation.grains[i].lin_vel - simulation.grains[j].lin_vel
+ # Contact kinematics (2d)
+ vel_lin = simulation.grains[i].lin_vel[1:2] -
+ simulation.grains[j].lin_vel[1:2]
vel_n = dot(vel_lin, n)
if !simulation.grains[i].rotating && !simulation.grains[j].rotating
@@ -152,25 +154,38 @@ function interactGrains!(simulation::Simulation, i::Int, j::Int, ic::Int)
if rotation
vel_t = dot(vel_lin, t) -
- harmonicMean(r_i, r_j)*(simulation.grains[i].ang_vel +
- simulation.grains[j].ang_vel)
+ harmonicMean(r_i, r_j)*(simulation.grains[i].ang_vel[3] +
+ simulation.grains[j].ang_vel[3])
# Correct old tangential displacement for contact rotation and add new
- δ_t0 = simulation.grains[i].contact_parallel_displacement[ic]
- δ_t = dot(t, δ_t0 - (n*dot(n, δ_t0))) + vel_t*simulation.time_step
+ δ_t0 = simulation.grains[i].contact_parallel_displacement[ic][1:2]
+ if !simulation.grains[i].compressive_failure[ic]
+ δ_t = dot(t, δ_t0 - (n*dot(n, δ_t0))) + vel_t*simulation.time_step
+ else
+ # Use δ_t as total slip-distance vector
+ δ_t = δ_t0 .+ (vel_lin .+ vel_t.*t) .* simulation.time_step
+ end
- # Determine the contact rotation
- θ_t = simulation.grains[i].contact_rotation[ic] +
- (simulation.grains[j].ang_vel - simulation.grains[i].ang_vel) *
- simulation.time_step
+ # Determine the contact rotation (2d)
+ θ_t = simulation.grains[i].contact_rotation[ic][3] +
+ (simulation.grains[j].ang_vel[3] -
+ simulation.grains[i].ang_vel[3]) * simulation.time_step
end
# Effective radius
R_ij = harmonicMean(r_i, r_j)
+ # Determine which ice floe is the thinnest
+ h_min, idx_thinnest = findmin([simulation.grains[i].thickness,
+ simulation.grains[j].thickness])
+ idx_thinnest == 1 ? idx_thinnest = i : idx_thinnest = j
+
+ # Contact length along z
+ Lz_ij = h_min
+
# Contact area
- A_ij = R_ij*min(simulation.grains[i].thickness,
- simulation.grains[j].thickness)
+ A_ij = R_ij*Lz_ij
+ simulation.grains[i].contact_area[ic] = A_ij
# Contact mechanical parameters
if simulation.grains[i].youngs_modulus > 0. &&
@@ -190,7 +205,7 @@ function interactGrains!(simulation::Simulation, i::Int, j::Int, ic::Int)
else # Micromechanical parameterization: k_n and k_t set explicitly
k_n = harmonicMean(simulation.grains[i].contact_stiffness_normal,
- simulation.grains[j].contact_stiffness_normal)
+ simulation.grains[j].contact_stiffness_normal)
if rotation
k_t = harmonicMean(simulation.grains[i].contact_stiffness_tangential,
@@ -223,26 +238,20 @@ function interactGrains!(simulation::Simulation, i::Int, j::Int, ic::Int)
end
else
- error("unknown contact_normal_rheology (k_n = $k_n, γ_n = $γ_n")
+ error("unknown contact_normal_rheology (k_n = $k_n, γ_n = $γ_n,
+ E = $E, A_ij = $A_ij, R_ij = $R_ij)
+ ")
end
- # Determine which grain is the weakest
- compressive_strength = simulation.grains[i].fracture_toughness *
- sqrt(simulation.grains[i].thickness)
- compressive_strength_j = simulation.grains[j].fracture_toughness *
- sqrt(simulation.grains[j].thickness)
- idx_weakest = i
- if compressive_strength_j < compressive_strength
- compressive_strength = compressive_strength_j
- idx_weakest = j
- end
+ # Determine the compressive strength in Pa by the contact thickness and the
+ # minimum compressive strength [N]
+ compressive_strength = min(simulation.grains[i].fracture_toughness,
+ simulation.grains[j].fracture_toughness) *
+ Lz_ij^1.5
# Grain-pair moment of inertia [m^4]
- if rotation
- I_ij = 2.0/3.0*R_ij^3*min(simulation.grains[i].thickness,
- simulation.grains[j].thickness)
- end
-
+ I_ij = 2.0/3.0*R_ij^3*min(simulation.grains[i].thickness,
+ simulation.grains[j].thickness)
# Contact tensile strength increases linearly with contact age until
# tensile stress exceeds tensile strength.
@@ -264,7 +273,7 @@ function interactGrains!(simulation::Simulation, i::Int, j::Int, ic::Int)
# Reset contact age (breaking bond) if bond strength is exceeded
if rotation
- if δ_n >= 0.0 && abs(force_n)/A_ij + abs(M_t)*R_ij/I_ij > tensile_strength
+ if δ_n >= 0.0 && norm(force_n)/A_ij + norm(M_t)*R_ij/I_ij > tensile_strength
force_n = 0.
force_t = 0.
simulation.grains[i].contacts[ic] = 0 # remove contact
@@ -272,7 +281,7 @@ function interactGrains!(simulation::Simulation, i::Int, j::Int, ic::Int)
simulation.grains[j].n_contacts -= 1
end
else
- if δ_n >= 0.0 && abs(force_n)/A_ij > tensile_strength
+ if δ_n >= 0.0 && norm(force_n)/A_ij > tensile_strength
force_n = 0.
force_t = 0.
simulation.grains[i].contacts[ic] = 0 # remove contact
@@ -281,36 +290,20 @@ function interactGrains!(simulation::Simulation, i::Int, j::Int, ic::Int)
end
end
- # Limit compressive stress if the prefactor is set to a positive value
- if δ_n <= 0.0 && compressive_strength > 0. &&
- abs(force_n) >= compressive_strength
-
- # Determine the overlap distance where yeild stress is reached
- δ_n_yield = -compressive_strength*A_ij/k_n
-
- # Determine the excess overlap distance relative to yield
- Δr = abs(δ_n) - abs(δ_n_yield)
-
- # Adjust radius and thickness of the weakest grain
- simulation.grains[idx_weakest].contact_radius -= Δr
- simulation.grains[idx_weakest].areal_radius -= Δr
- simulation.grains[idx_weakest].thickness += 1.0/(π*Δr)
- end
-
- if rotation
+ if rotation && !simulation.grains[i].compressive_failure[ic]
if k_t ≈ 0. && γ_t ≈ 0.
# do nothing
elseif k_t ≈ 0. && γ_t > 0.
- force_t = abs(γ_t * vel_t)
+ force_t = norm(γ_t * vel_t)
# Coulomb slip
- if force_t > μ_d_minimum*abs(force_n)
- force_t = μ_d_minimum*abs(force_n)
+ if force_t > μ_d_minimum*norm(force_n)
+ force_t = μ_d_minimum*norm(force_n)
# Nguyen et al 2009 "Thermomechanical modelling of friction effects
# in granular flows using the discrete element method"
- E_shear = abs(force_t)*abs(vel_t)*simulation.time_step
+ E_shear = norm(force_t)*norm(vel_t)*simulation.time_step
# Assume equal thermal properties
simulation.grains[i].thermal_energy += 0.5*E_shear
@@ -325,13 +318,13 @@ function interactGrains!(simulation::Simulation, i::Int, j::Int, ic::Int)
force_t = -k_t*δ_t - γ_t*vel_t
# Coulomb slip
- if abs(force_t) > μ_d_minimum*abs(force_n)
- force_t = μ_d_minimum*abs(force_n)*force_t/abs(force_t)
+ if norm(force_t) > μ_d_minimum*norm(force_n)
+ force_t = μ_d_minimum*norm(force_n)*force_t/norm(force_t)
δ_t = (-force_t - γ_t*vel_t)/k_t
- # Nguyen et al 2009 "Thermomechanical modelling of friction effects
- # in granular flows using the discrete element method"
- E_shear = abs(force_t)*abs(vel_t)*simulation.time_step
+ # Nguyen et al 2009 "Thermomechanical modelling of friction
+ # effects in granular flows using the discrete element method"
+ E_shear = norm(force_t)*norm(vel_t)*simulation.time_step
# Assume equal thermal properties
simulation.grains[i].thermal_energy += 0.5*E_shear
@@ -343,28 +336,104 @@ function interactGrains!(simulation::Simulation, i::Int, j::Int, ic::Int)
end
end
+ just_failed = false
+ # Detect compressive failure if compressive force exceeds compressive
+ # strength
+ if δ_n <= 0.0 && compressive_strength > 0. &&
+ !simulation.grains[i].compressive_failure[ic] &&
+ norm(force_n.*n .+ force_t.*t) >= compressive_strength
+
+ # Register that compressive failure has occurred for this contact
+ simulation.grains[i].compressive_failure[ic] = true
+ just_failed = true
+ end
+
+ # Overwrite previous force results if compressive failure occurred
+ if simulation.grains[i].compressive_failure[ic] && δ_n < 0.0
+
+ # Normal stress on horizontal interface, assuming bending stresses are
+ # negligible (ocean density and gravity are hardcoded for now)
+ σ_n = (1e3 - simulation.grains[i].density) * 9.81 *
+ (simulation.grains[i].thickness + simulation.grains[j].thickness)
+
+ # Horizontal overlap area (two overlapping circles)
+ if 2.0*min(r_i, r_j) <= abs(δ_n)
+
+ # Complete overlap
+ A_h = π*min(r_i, r_j)^2
+
+ else
+ # Partial overlap, eq 14 from
+ # http://mathworld.wolfram.com/Circle-CircleIntersection.html
+ A_h = r_i^2.0*acos((dist^2.0 + r_i^2.0 - r_j^2.0)/(2.0*dist*r_i)) +
+ r_j^2.0*acos((dist^2.0 + r_j^2.0 - r_i^2.0)/(2.0*dist*r_j)) -
+ 0.5*sqrt((-dist + r_i + r_j)*
+ ( dist + r_i - r_j)*
+ ( dist - r_i + r_j)*
+ ( dist + r_i + r_j))
+ end
+ simulation.grains[i].contact_area[ic] = A_h
+
+ if just_failed
+ # Use δ_t as travel distance on horizontal slip surface, and set the
+ # original displacement to Coulomb-frictional limit in the direction
+ # of the pre-failure force
+ F = force_n.*n .+ force_t.*t
+ σ_t = μ_d_minimum*norm(σ_n).*F/norm(F)
+ δ_t = σ_t*A_h/(-k_t)
+
+ # Save energy loss in E_shear
+ E_shear = norm(F)*norm(vel_lin)*simulation.time_step
+ simulation.grains[i].thermal_energy += 0.5*E_shear
+ simulation.grains[j].thermal_energy += 0.5*E_shear
+
+ else
+
+ # Find horizontal stress on slip interface
+ σ_t = -k_t*δ_t/A_h
+
+ # Check for Coulomb-failure on the slip interface
+ if norm(σ_t) > μ_d_minimum*norm(σ_n)
+
+ σ_t = μ_d_minimum*norm(σ_n)*σ_t/norm(σ_t)
+ δ_t = σ_t*A_h/(-k_t)
+
+ E_shear = norm(σ_t*A_h)*norm(vel_lin)*simulation.time_step
+ simulation.grains[i].thermal_energy += 0.5*E_shear
+ simulation.grains[j].thermal_energy += 0.5*E_shear
+ end
+ end
+
+ force_n = dot(σ_t, n)*A_h
+ force_t = dot(σ_t, t)*A_h
+ end
+
# Break bond under extension through bending failure
if δ_n < 0.0 && tensile_strength > 0.0 && rotation &&
- shear_strength > 0.0 && abs(force_t)/A_ij > shear_strength
+ shear_strength > 0.0 && norm(force_t)/A_ij > shear_strength
force_n = 0.
force_t = 0.
simulation.grains[i].contacts[ic] = 0 # remove contact
simulation.grains[i].n_contacts -= 1
simulation.grains[j].n_contacts -= 1
+ else
+ simulation.grains[i].contact_age[ic] += simulation.time_step
end
- simulation.grains[i].contact_age[ic] += simulation.time_step
-
if rotation
- simulation.grains[i].contact_parallel_displacement[ic] = δ_t*t
- simulation.grains[i].contact_rotation[ic] = θ_t
+ if simulation.grains[i].compressive_failure[ic]
+ simulation.grains[i].contact_parallel_displacement[ic] = vecTo3d(δ_t)
+ else
+ simulation.grains[i].contact_parallel_displacement[ic] = vecTo3d(δ_t.*t)
+ end
+ simulation.grains[i].contact_rotation[ic] = [0., 0., θ_t]
- simulation.grains[i].force += force_n*n + force_t*t;
- simulation.grains[j].force -= force_n*n + force_t*t;
+ simulation.grains[i].force += vecTo3d(force_n.*n + force_t.*t);
+ simulation.grains[j].force -= vecTo3d(force_n.*n + force_t.*t);
- simulation.grains[i].torque += -force_t*R_ij + M_t
- simulation.grains[j].torque += -force_t*R_ij - M_t
+ simulation.grains[i].torque[3] += -force_t*R_ij + M_t
+ simulation.grains[j].torque[3] += -force_t*R_ij - M_t
else
simulation.grains[i].force += force_n*n;
simulation.grains[j].force -= force_n*n;
diff --git a/src/io.jl b/src/io.jl
@@ -393,6 +393,8 @@ function writeGrainVTK(simulation::Simulation,
"Fixed but allow (x) acceleration [-]")
WriteVTK.vtk_point_data(vtkfile, ifarr.allow_y_acc,
"Fixed but allow (y) acceleration [-]")
+ WriteVTK.vtk_point_data(vtkfile, ifarr.allow_z_acc,
+ "Fixed but allow (z) acceleration [-]")
WriteVTK.vtk_point_data(vtkfile, ifarr.rotating, "Free to rotate [-]")
WriteVTK.vtk_point_data(vtkfile, ifarr.enabled, "Enabled [-]")
@@ -481,8 +483,10 @@ function writeGrainInteractionVTK(simulation::Simulation,
contact_stiffness = Float64[]
tensile_stress = Float64[]
shear_displacement = Vector{Float64}[]
- contact_rotation = Float64[]
+ contact_rotation = Vector{Float64}[]
contact_age = Float64[]
+ contact_area = Float64[]
+ compressive_failure = Int[]
for i=1:length(simulation.grains)
for ic=1:simulation.Nc_max
@@ -549,6 +553,9 @@ function writeGrainInteractionVTK(simulation::Simulation,
simulation.grains[i].contact_rotation[ic])
push!(contact_age, simulation.grains[i].contact_age[ic])
+ push!(contact_area, simulation.grains[i].contact_area[ic])
+ push!(compressive_failure,
+ Int(simulation.grains[i].compressive_failure[ic]))
end
end
end
@@ -577,7 +584,7 @@ function writeGrainInteractionVTK(simulation::Simulation,
for i=1:length(i1)
write(f, "$(inter_particle_vector[i][1]) ")
write(f, "$(inter_particle_vector[i][2]) ")
- write(f, "0.0 ")
+ write(f, "$(inter_particle_vector[i][3]) ")
end
write(f, "\n")
write(f, " </DataArray>\n")
@@ -587,7 +594,7 @@ function writeGrainInteractionVTK(simulation::Simulation,
for i=1:length(i1)
@inbounds write(f, "$(shear_displacement[i][1]) ")
@inbounds write(f, "$(shear_displacement[i][2]) ")
- write(f, "0.0 ")
+ @inbounds write(f, "$(shear_displacement[i][3]) ")
end
write(f, "\n")
write(f, " </DataArray>\n")
@@ -637,10 +644,12 @@ function writeGrainInteractionVTK(simulation::Simulation,
write(f, " </DataArray>\n")
write(f, " <DataArray type=\"Float32\" " *
- "Name=\"Contact rotation [rad]\" NumberOfComponents=\"1\"
+ "Name=\"Contact rotation [rad]\" NumberOfComponents=\"3\"
format=\"ascii\">\n")
for i=1:length(i1)
- @inbounds write(f, "$(contact_rotation[i]) ")
+ @inbounds write(f, "$(contact_rotation[i][1]) ")
+ @inbounds write(f, "$(contact_rotation[i][2]) ")
+ @inbounds write(f, "$(contact_rotation[i][3]) ")
end
write(f, "\n")
write(f, " </DataArray>\n")
@@ -654,6 +663,24 @@ function writeGrainInteractionVTK(simulation::Simulation,
write(f, "\n")
write(f, " </DataArray>\n")
+ write(f, " <DataArray type=\"Float32\" " *
+ "Name=\"Contact area [m*m]\" NumberOfComponents=\"1\"
+ format=\"ascii\">\n")
+ for i=1:length(i1)
+ @inbounds write(f, "$(contact_area[i]) ")
+ end
+ write(f, "\n")
+ write(f, " </DataArray>\n")
+
+ write(f, " <DataArray type=\"Float32\" " *
+ "Name=\"Compressive failure [-]\" NumberOfComponents=\"1\"
+ format=\"ascii\">\n")
+ for i=1:length(i1)
+ @inbounds write(f, "$(Int(compressive_failure[i])) ")
+ end
+ write(f, "\n")
+ write(f, " </DataArray>\n")
+
write(f, " </CellData>\n")
write(f, " <Points>\n")
@@ -662,7 +689,7 @@ function writeGrainInteractionVTK(simulation::Simulation,
write(f, " <DataArray type=\"Float32\" Name=\"Points\" " *
"NumberOfComponents=\"3\" format=\"ascii\">\n")
for i in simulation.grains
- @inbounds write(f, "$(i.lin_pos[1]) $(i.lin_pos[2]) 0.0 ")
+ @inbounds write(f, "$(i.lin_pos[1]) $(i.lin_pos[2]) $(i.lin_pos[3]) ")
end
write(f, "\n")
write(f, " </DataArray>\n")
@@ -868,6 +895,7 @@ imagegrains.PointArrayStatus = [
'Fixed in space [-]',
'Fixed but allow (x) acceleration [-]',
'Fixed but allow (y) acceleration [-]',
+'Fixed but allow (z) acceleration [-]',
'Free to rotate [-]',
'Enabled [-]',
'Contact stiffness (normal) [N m^-1]',
@@ -1103,6 +1131,12 @@ function render(simulation::Simulation; pvpython::String="pvpython",
if isfile("$(simulation.id)/$(simulation.id).mp4")
rm("$(simulation.id)/$(simulation.id).avi")
end
+ catch return_signal
+ if isa(return_signal, Base.UVError)
+ Compat.@warn "Could not run external ffmpeg command, " *
+ "skipping conversion from " *
+ "$(simulation.id)/$(simulation.id).avi to mp4."
+ end
end
end
@@ -1131,7 +1165,7 @@ function render(simulation::Simulation; pvpython::String="pvpython",
end
catch return_signal
if isa(return_signal, Base.UVError)
- Compat.@info "Skipping gif merge since `$convert` " *
+ Compat.@warn "Skipping gif merge since `$convert` " *
"was not found."
end
end
@@ -1321,8 +1355,9 @@ function plotGrains(sim::Simulation;
contact_stiffness = Float64[]
tensile_stress = Float64[]
shear_displacement = Vector{Float64}[]
- contact_rotation = Float64[]
+ contact_rotation = Vector{Float64}[]
contact_age = Float64[]
+ compressive_failure = Int[]
for i=1:length(sim.grains)
for ic=1:sim.Nc_max
if sim.grains[i].contacts[ic] > 0
@@ -1374,6 +1409,8 @@ function plotGrains(sim::Simulation;
sim.grains[i].contact_rotation[ic])
push!(contact_age, sim.grains[i].contact_age[ic])
+ push!(compressive_failure,
+ sim.grains[i].compressive_failure[ic])
end
end
end
diff --git a/src/ocean.jl b/src/ocean.jl
@@ -388,10 +388,10 @@ function applyOceanDragToGrain!(grain::GrainCylindrical,
drag_force = ρ_o * π *
(2.0*grain.ocean_drag_coeff_vert*grain.areal_radius*draft +
grain.ocean_drag_coeff_horiz*grain.areal_radius^2.0) *
- ([u, v] - grain.lin_vel)*norm([u, v] - grain.lin_vel)
+ ([u, v] - grain.lin_vel[1:2])*norm([u, v] - grain.lin_vel[1:2])
- grain.force += drag_force
- grain.ocean_stress = drag_force/grain.horizontal_surface_area
+ grain.force += vecTo3d(drag_force)
+ grain.ocean_stress = vecTo3d(drag_force/grain.horizontal_surface_area)
nothing
end
@@ -407,11 +407,12 @@ function applyOceanVorticityToGrain!(grain::GrainCylindrical,
ρ_o = 1000. # ocean density
draft = grain.thickness - freeboard # height of submerged thickness
- grain.torque +=
+ grain.torque[3] +=
π * grain.areal_radius^4. * ρ_o *
(grain.areal_radius/5. * grain.ocean_drag_coeff_horiz +
draft * grain.ocean_drag_coeff_vert) *
- abs(.5 * ocean_curl - grain.ang_vel) * (.5 * ocean_curl - grain.ang_vel)
+ norm(.5 * ocean_curl - grain.ang_vel[3]) *
+ (.5 * ocean_curl - grain.ang_vel[3])
nothing
end
diff --git a/src/packing.jl b/src/packing.jl
@@ -124,7 +124,7 @@ function generateNeighboringPoint(x_i::Vector, r_i::Real,
#R = rand() * (r_i + r_j) * max_padding_factor + 2. * (r_i + r_j)
R = r_i + r_j + padding
T = rand() * 2. * pi
- return [x_i[1] + R * sin(T), x_i[2] + R * cos(T)], r_j
+ return [x_i[1] + R * sin(T), x_i[2] + R * cos(T), x_i[3]], r_j
end
function generateRandomDirection()
@@ -132,7 +132,7 @@ function generateRandomDirection()
end
function getPositionDistancedFromPoint(T::Real, x_i::Vector, dist::Real)
- return [x_i[1] + dist * sin(T), x_i[2] + dist * cos(T)]
+ return [x_i[1] + dist * sin(T), x_i[2] + dist * cos(T), x_i[3]]
end
export irregularPacking!
@@ -224,7 +224,7 @@ function irregularPacking!(simulation::Simulation;
# and add it to the active list. If after `sample_limit` attempts no such
# point is found, instead remove `i` from the active list.
j = 0;
- x_active = zeros(2); x_candidate = zeros(2);
+ x_active = zeros(3); x_candidate = zeros(3);
r_active = 0.; r_candidate = 0.; T = 0.
n = 0
neighbor_found = false
@@ -456,9 +456,9 @@ function rasterMap(sim::Simulation, dx::Real)
#i, j = Int.(floor.((grain.lin_pos - origo) ./ dx)) + [1,1]
# Find corner indexes for box spanning the grain
- min_i, min_j = Int.(floor.((grain.lin_pos - origo .-
+ min_i, min_j = Int.(floor.((grain.lin_pos[1:2] - origo .-
grain.contact_radius) ./ dx)) .+ [1,1]
- max_i, max_j = Int.(floor.((grain.lin_pos - origo .+
+ max_i, max_j = Int.(floor.((grain.lin_pos[1:2] - origo .+
grain.contact_radius) ./ dx)) .+ [1,1]
# For each cell in box, check if the grain is contained
@@ -466,7 +466,7 @@ function rasterMap(sim::Simulation, dx::Real)
for j = min_j:max_j
cell_mid_point = dx .* Vector{Float64}([i,j]) .- 0.5 * dx
- if (norm(cell_mid_point - grain.lin_pos) -
+ if (norm(cell_mid_point - grain.lin_pos[1:2]) -
grain.contact_radius < dist)
occupied[i,j] = true
end
diff --git a/src/simulation.jl b/src/simulation.jl
@@ -259,7 +259,7 @@ export zeroForcesAndTorques!
function zeroForcesAndTorques!(simulation::Simulation)
for grain in simulation.grains
grain.force .= grain.external_body_force
- grain.torque = 0.
+ grain.torque = zeros(3)
grain.pressure = 0.
end
for wall in simulation.walls
diff --git a/src/temporal_integration.jl b/src/temporal_integration.jl
@@ -51,7 +51,7 @@ Use a two-term Taylor expansion for integrating the kinematic degrees of freedom
for a `grain`.
"""
function updateGrainKinematicsTwoTermTaylor!(grain::GrainCylindrical,
- simulation::Simulation)
+ simulation::Simulation)
grain.lin_acc = grain.force/grain.mass
grain.ang_acc = grain.torque/grain.moment_of_inertia
@@ -62,9 +62,12 @@ function updateGrainKinematicsTwoTermTaylor!(grain::GrainCylindrical,
if !grain.allow_y_acc
grain.lin_acc[2] = 0.
end
- grain.ang_acc = 0.
+ if !grain.allow_z_acc
+ grain.lin_acc[3] = 0.
+ end
+ grain.ang_acc = zeros(3)
elseif !grain.rotating
- grain.ang_acc = 0.
+ grain.ang_acc = zeros(3)
end
grain.lin_pos +=
@@ -90,11 +93,11 @@ Use a three-term Taylor expansion for integrating the kinematic degrees of
freedom for a `grain`.
"""
function updateGrainKinematicsThreeTermTaylor!(grain::GrainCylindrical,
- simulation::Simulation)
+ simulation::Simulation)
if simulation.time_iteration == 0
- lin_acc_0 = zeros(2)
- ang_acc_0 = 0.
+ lin_acc_0 = zeros(3)
+ ang_acc_0 = zeros(3)
else
lin_acc_0 = grain.lin_acc
ang_acc_0 = grain.ang_acc
@@ -110,9 +113,12 @@ function updateGrainKinematicsThreeTermTaylor!(grain::GrainCylindrical,
if !grain.allow_y_acc
grain.lin_acc[2] = 0.
end
- grain.ang_acc = 0.
+ if !grain.allow_z_acc
+ grain.lin_acc[3] = 0.
+ end
+ grain.ang_acc = zeros(3)
elseif !grain.rotating
- grain.ang_acc = 0.
+ grain.ang_acc = zeros(3)
end
# Temporal gradient in acceleration using backwards differences
diff --git a/src/util.jl b/src/util.jl
@@ -42,3 +42,31 @@ function harmonicMean(a::Number, b::Number)::Number
return 2. * a * b / (a + b)
end
end
+
+export vecTo3d
+"""
+ function vecTo3d(input, fill)
+
+Convert a scalar or 2d vector to 3d by filling the missing component with the
+value `fill`. The returned 3-component vector is a Vector (or 1d Array) of
+the same type as the input.
+
+# Arguments
+* `input`: a scalar or two-component vector.
+* `fill::Real`: value to use for third
+"""
+function vecTo3d(input::Any; fill::Real = 0.0)
+ if length(input) > 3
+ error("vecTo3d requires a scalar or input vector to of length 3 or " *
+ "less, but input is $input")
+ elseif length(input) == 3
+ return input
+ elseif length(input) == 2
+ return [input[1], input[2], typeof(input[1])(fill)]
+ elseif length(input) == 1
+ return [input[1], typeof(input[1])(fill), typeof(input[1])(fill)]
+ else
+ error("input not understood: $input")
+ end
+end
+
diff --git a/src/wall.jl b/src/wall.jl
@@ -16,9 +16,10 @@ The only required arguments are
# Arguments
* `simulation::Simulation`: the simulation object where the wall should be
added to.
-* `normal::Vector{Float64}`: 2d vector denoting the normal to the wall [m]. The
+* `normal::Vector{Float64}`: 3d vector denoting the normal to the wall [m]. The
wall will only interact in the opposite direction of this vector, so the
- normal vector should point in the direction of the grains.
+ normal vector should point in the direction of the grains. If a 2d vector is
+ passed, the third (z) component is set to zero.
* `pos::Float64`: position along axis parallel to the normal vector [m].
* `bc::String="fixed"`: boundary condition, possible values are `"fixed"`
(default), `"normal stress"`, or `"velocity"`.
@@ -50,13 +51,13 @@ wall-face normal of `[1., 0.]` (wall along *y* and normal to *x*), a position of
`1.5` meter:
```julia
-Granular.addWallLinearFrictionless!(sim, [1., 0.], 1.5)
+Granular.addWallLinearFrictionless!(sim, [1., 0., 0.], 1.5)
```
The following example creates a wall with a velocity of 0.5 m/s towards *-y*:
```julia
-Granular.addWallLinearFrictionless!(sim, [0., 1.], 1.5,
+Granular.addWallLinearFrictionless!(sim, [0., 1., 0.], 1.5,
bc="velocity",
vel=-0.5)
```
@@ -65,7 +66,7 @@ To create a wall parallel to the *y* axis pushing downwards with a constant
normal stress of 100 kPa, starting at a position of y = 3.5 m:
```julia
-Granular.addWallLinearFrictionless!(sim, [0., 1.], 3.5,
+Granular.addWallLinearFrictionless!(sim, [0., 1., 0.], 3.5,
bc="normal stress",
normal_stress=100e3)
```
@@ -85,20 +86,19 @@ function addWallLinearFrictionless!(simulation::Simulation,
verbose::Bool=true)
# Check input values
- if length(normal) != 2
- error("Wall normal must be a two-element array (normal = " *
- "$normal)")
+ if length(normal) != 3
+ normal = vecTo3d(normal)
end
if bc != "fixed" && bc != "velocity" && bc != "normal stress"
error("Wall BC must be 'fixed', 'velocity', or 'normal stress'.")
end
- if !(normal ≈ [1., 0.]) && !(normal ≈ [0., 1.])
+ if !(normal ≈ [1., 0., 0.]) && !(normal ≈ [0., 1., 0.])
error("Currently only walls with normals orthogonal to the " *
"coordinate system are allowed, i.e. normals parallel to the " *
"x or y axes. Accepted values for `normal` " *
- "are [1., 0.] and [0., 1.]. The passed normal was $normal")
+ "are [1., 0., 0.] and [0., 1., 0.]. The passed normal was $normal")
end
# if not set, set wall mass to equal the mass of all grains.
@@ -171,10 +171,10 @@ Returns the surface area of the wall given the grid size and its index.
"""
function getWallSurfaceArea(sim::Simulation, wall_index::Integer)
- if sim.walls[wall_index].normal ≈ [1., 0.]
+ if sim.walls[wall_index].normal ≈ [1., 0., 0.]
return (sim.ocean.yq[end,end] - sim.ocean.yq[1,1]) *
sim.walls[wall_index].thickness
- elseif sim.walls[wall_index].normal ≈ [0., 1.]
+ elseif sim.walls[wall_index].normal ≈ [0., 1., 0.]
return (sim.ocean.xq[end,end] - sim.ocean.xq[1,1]) *
sim.walls[wall_index].thickness
else
@@ -185,9 +185,9 @@ end
function getWallSurfaceArea(sim::Simulation, normal::Vector{Float64},
thickness::Float64)
- if normal ≈ [1., 0.]
+ if length(normal) == 3 && normal ≈ [1., 0., 0.]
return (sim.ocean.yq[end,end] - sim.ocean.yq[1,1]) * thickness
- elseif normal ≈ [0., 1.]
+ elseif length(normal) == 3 && normal ≈ [0., 1., 0.]
return (sim.ocean.xq[end,end] - sim.ocean.xq[1,1]) * thickness
else
error("Wall normal not understood")
diff --git a/test/atmosphere.jl b/test/atmosphere.jl
@@ -3,8 +3,6 @@
# Check if atmosphere-specific functions and grid operations are functioning
# correctly
-Compat.@info "#### $(basename(@__FILE__)) ####"
-
Compat.@info "Testing regular grid generation"
sim = Granular.createSimulation()
sim.atmosphere = Granular.createRegularAtmosphereGrid([6, 6, 6], [1., 1., 1.])
@@ -54,15 +52,15 @@ E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
Compat.@info "Testing vortex interaction (static atmosphere)"
sim = deepcopy(sim_init)
-sim.grains[1].ang_vel = 0.1
+sim.grains[1].ang_vel[3] = 0.1
E_kin_lin_init = Granular.totalGrainKineticTranslationalEnergy(sim)
E_kin_rot_init = Granular.totalGrainKineticRotationalEnergy(sim)
Granular.run!(sim, verbose=false)
E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
@test E_kin_rot_init > E_kin_rot_final # energy lost to atmosphere
-@test sim.grains[1].ang_vel > 0. # check angular velocity orientation
-@test sim.grains[1].ang_pos > 0. # check angular position orientation
+@test sim.grains[1].ang_vel[3] > 0. # check angular velocity orientation
+@test sim.grains[1].ang_pos[3] > 0. # check angular position orientation
@test E_kin_lin_init ≈ E_kin_lin_final # no linear velocity gained
Compat.@info "Testing vortex interaction (static ice floe)"
@@ -77,8 +75,8 @@ E_kin_rot_init = Granular.totalGrainKineticRotationalEnergy(sim)
Granular.run!(sim, verbose=false)
E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
-@test sim.grains[1].ang_vel > 0. # check angular velocity orientation
-@test sim.grains[1].ang_pos > 0. # check angular position orientation
+@test sim.grains[1].ang_vel[3] > 0. # check angular velocity orientation
+@test sim.grains[1].ang_pos[3] > 0. # check angular position orientation
@test E_kin_rot_init < E_kin_rot_final # rotation after due to atm vortex
@test E_kin_lin_init ≈ E_kin_lin_final # no linear velocity gained
@@ -93,8 +91,8 @@ E_kin_rot_init = Granular.totalGrainKineticRotationalEnergy(sim)
Granular.run!(sim, verbose=false)
E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
-@test sim.grains[1].ang_vel < 0. # check angular velocity orientation
-@test sim.grains[1].ang_pos < 0. # check angular position orientation
+@test sim.grains[1].ang_vel[3] < 0. # check angular velocity orientation
+@test sim.grains[1].ang_pos[3] < 0. # check angular position orientation
@test E_kin_rot_init < E_kin_rot_final # rotation after due to atm vortex
@test E_kin_lin_init ≈ E_kin_lin_final # no linear velocity gained
@@ -109,8 +107,8 @@ E_kin_rot_init = Granular.totalGrainKineticRotationalEnergy(sim)
Granular.run!(sim, verbose=false)
E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
-@test sim.grains[1].ang_vel < 0. # check angular velocity orientation
-@test sim.grains[1].ang_pos < 0. # check angular position orientation
+@test sim.grains[1].ang_vel[3] < 0. # check angular velocity orientation
+@test sim.grains[1].ang_pos[3] < 0. # check angular position orientation
@test E_kin_rot_init < E_kin_rot_final # rotation after due to atm vortex
@test E_kin_lin_init ≈ E_kin_lin_final # no linear velocity gained
@@ -125,8 +123,8 @@ E_kin_rot_init = Granular.totalGrainKineticRotationalEnergy(sim)
Granular.run!(sim, verbose=false)
E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
-@test sim.grains[1].ang_vel > 0. # check angular velocity orientation
-@test sim.grains[1].ang_pos > 0. # check angular position orientation
+@test sim.grains[1].ang_vel[3] > 0. # check angular velocity orientation
+@test sim.grains[1].ang_pos[3] > 0. # check angular position orientation
@test E_kin_rot_init < E_kin_rot_final # rotation after due to atm vortex
@test E_kin_lin_init ≈ E_kin_lin_final # no linear velocity gained
diff --git a/test/cohesion.jl b/test/cohesion.jl
@@ -5,8 +5,6 @@ import Granular
# Check for conservation of kinetic energy (=momentum) during a normal collision
# between two ice cylindrical grains
-Compat.@info "#### $(basename(@__FILE__)) ####"
-
verbose=false
sim_init = Granular.createSimulation()
@@ -21,6 +19,7 @@ sim = deepcopy(sim_init)
Granular.setTotalTime!(sim, 10.)
sim.time_step = 1.
Granular.run!(sim, verbose=verbose)
+Granular.removeSimulationFiles(sim)
@test sim.grains[1].contact_age[1] ≈ sim.time
Compat.@info "# Check if bonds add tensile strength"
@@ -31,12 +30,13 @@ sim.grains[1].lin_vel[1] = 0.1
Granular.setTimeStep!(sim)
Granular.setTotalTime!(sim, 10.)
Granular.run!(sim, verbose=verbose)
+Granular.removeSimulationFiles(sim)
@test sim.grains[1].lin_vel[1] > 0.
@test sim.grains[1].lin_vel[2] ≈ 0.
@test sim.grains[2].lin_vel[1] > 0.
@test sim.grains[2].lin_vel[2] ≈ 0.
-@test sim.grains[1].ang_vel ≈ 0.
-@test sim.grains[2].ang_vel ≈ 0.
+@test sim.grains[1].ang_vel ≈ zeros(3)
+@test sim.grains[2].ang_vel ≈ zeros(3)
Compat.@info "# Add shear strength and test bending resistance (one grain rotating)"
sim = Granular.createSimulation(id="cohesion")
@@ -44,7 +44,7 @@ Granular.addGrainCylindrical!(sim, [0., 0.], 10.1, 1., tensile_strength=500e3,
shear_strength=500e3)
Granular.addGrainCylindrical!(sim, [20., 0.], 10., 1., tensile_strength=500e3,
shear_strength=500e3)
-sim.grains[1].ang_vel = 0.1
+sim.grains[1].ang_vel[3] = 0.1
Granular.findContacts!(sim) # make sure contact is registered
sim.grains[1].contact_radius=10.0 # decrease radius so there isn't compression
E_kin_lin_init = Granular.totalGrainKineticTranslationalEnergy(sim)
@@ -53,12 +53,13 @@ Granular.setTimeStep!(sim)
Granular.setTotalTime!(sim, 5.)
#Granular.setOutputFileInterval!(sim, 0.1)
Granular.run!(sim, verbose=verbose)
+Granular.removeSimulationFiles(sim)
@test sim.grains[1].lin_vel[1] ≈ 0.
@test sim.grains[1].lin_vel[2] ≈ 0.
@test sim.grains[2].lin_vel[1] ≈ 0.
@test sim.grains[2].lin_vel[2] ≈ 0.
-@test sim.grains[1].ang_vel != 0.
-@test sim.grains[2].ang_vel != 0.
+@test sim.grains[1].ang_vel[3] != 0.
+@test sim.grains[2].ang_vel[3] != 0.
E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
E_therm_final = Granular.totalGrainThermalEnergy(sim)
@@ -71,7 +72,7 @@ Granular.addGrainCylindrical!(sim, [0., 0.], 10.1, 1., tensile_strength=500e3,
shear_strength=500e3)
Granular.addGrainCylindrical!(sim, [20., 0.], 10., 1., tensile_strength=500e3,
shear_strength=500e3)
-sim.grains[2].ang_vel = 0.1
+sim.grains[2].ang_vel[3] = 0.1
Granular.findContacts!(sim) # make sure contact is registered
sim.grains[1].contact_radius=10.0 # decrease radius so there isn't compression
E_kin_lin_init = Granular.totalGrainKineticTranslationalEnergy(sim)
@@ -80,12 +81,13 @@ Granular.setTimeStep!(sim)
Granular.setTotalTime!(sim, 5.)
#Granular.setOutputFileInterval!(sim, 0.1)
Granular.run!(sim, verbose=verbose)
+Granular.removeSimulationFiles(sim)
@test sim.grains[1].lin_vel[1] ≈ 0.
@test sim.grains[1].lin_vel[2] ≈ 0.
@test sim.grains[2].lin_vel[1] ≈ 0.
@test sim.grains[2].lin_vel[2] ≈ 0.
-@test sim.grains[1].ang_vel != 0.
-@test sim.grains[2].ang_vel != 0.
+@test sim.grains[1].ang_vel[3] != 0.
+@test sim.grains[2].ang_vel[3] != 0.
E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
E_therm_final = Granular.totalGrainThermalEnergy(sim)
@@ -98,8 +100,8 @@ Granular.addGrainCylindrical!(sim, [0., 0.], 10.0000001, 1., tensile_strength=50
shear_strength=500e3)
Granular.addGrainCylindrical!(sim, [20., 0.], 10., 1., tensile_strength=500e3,
shear_strength=500e3)
-sim.grains[1].ang_vel = 0.1
-sim.grains[2].ang_vel = -0.1
+sim.grains[1].ang_vel[3] = 0.1
+sim.grains[2].ang_vel[3] = -0.1
Granular.findContacts!(sim) # make sure contact is registered
sim.grains[1].contact_radius=10.0 # decrease radius so there isn't compression
E_kin_lin_init = Granular.totalGrainKineticTranslationalEnergy(sim)
@@ -108,12 +110,13 @@ Granular.setTimeStep!(sim)
Granular.setTotalTime!(sim, 5.)
#Granular.setOutputFileInterval!(sim, 0.1)
Granular.run!(sim, verbose=verbose)
+Granular.removeSimulationFiles(sim)
@test sim.grains[1].lin_vel[1] ≈ 0.
@test sim.grains[1].lin_vel[2] ≈ 0.
@test sim.grains[2].lin_vel[1] ≈ 0.
@test sim.grains[2].lin_vel[2] ≈ 0.
-@test sim.grains[1].ang_vel != 0.
-@test sim.grains[2].ang_vel != 0.
+@test sim.grains[1].ang_vel[3] != 0.
+@test sim.grains[2].ang_vel[3] != 0.
E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
E_therm_final = Granular.totalGrainThermalEnergy(sim)
@@ -126,8 +129,8 @@ Granular.addGrainCylindrical!(sim, [0., 0.], 10.0000001, 1., tensile_strength=50
shear_strength=500e3)
Granular.addGrainCylindrical!(sim, [20., 0.], 10., 1., tensile_strength=500e3,
shear_strength=500e3)
-sim.grains[1].ang_vel = 100
-sim.grains[2].ang_vel = -100
+sim.grains[1].ang_vel[3] = 100
+sim.grains[2].ang_vel[3] = -100
Granular.findContacts!(sim) # make sure contact is registered
sim.grains[1].contact_radius=10.0 # decrease radius so there isn't compression
E_kin_lin_init = Granular.totalGrainKineticTranslationalEnergy(sim)
@@ -136,12 +139,13 @@ Granular.setTimeStep!(sim)
Granular.setTotalTime!(sim, 5.)
#Granular.setOutputFileInterval!(sim, 0.1)
Granular.run!(sim, verbose=verbose)
+Granular.removeSimulationFiles(sim)
@test sim.grains[1].lin_vel[1] ≈ 0.
@test sim.grains[1].lin_vel[2] ≈ 0.
@test sim.grains[2].lin_vel[1] ≈ 0.
@test sim.grains[2].lin_vel[2] ≈ 0.
-@test sim.grains[1].ang_vel != 0.
-@test sim.grains[2].ang_vel != 0.
+@test sim.grains[1].ang_vel[3] != 0.
+@test sim.grains[2].ang_vel[3] != 0.
E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
E_therm_final = Granular.totalGrainThermalEnergy(sim)
@@ -155,8 +159,8 @@ Granular.addGrainCylindrical!(sim, [0., 0.], 10.1, 1., tensile_strength=500e3,
shear_strength=50e3)
Granular.addGrainCylindrical!(sim, [20., 0.], 10., 1., tensile_strength=500e3,
shear_strength=50e3)
-sim.grains[1].ang_vel = 100
-sim.grains[2].ang_vel = 0.0
+sim.grains[1].ang_vel[3] = 100
+sim.grains[2].ang_vel[3] = 0.0
Granular.findContacts!(sim) # make sure contact is registered
sim.grains[1].contact_radius=10.0 # decrease radius so there isn't compression
E_kin_lin_init = Granular.totalGrainKineticTranslationalEnergy(sim)
@@ -165,12 +169,13 @@ Granular.setTimeStep!(sim)
Granular.setTotalTime!(sim, 5.)
#Granular.setOutputFileInterval!(sim, 0.1)
Granular.run!(sim, verbose=verbose)
+Granular.removeSimulationFiles(sim)
@test sim.grains[1].lin_vel[1] ≈ 0.
@test sim.grains[1].lin_vel[2] ≈ 0.
@test sim.grains[2].lin_vel[1] ≈ 0.
@test sim.grains[2].lin_vel[2] ≈ 0.
-@test sim.grains[1].ang_vel != 0.
-@test sim.grains[2].ang_vel != 0.
+@test sim.grains[1].ang_vel[3] != 0.
+@test sim.grains[2].ang_vel[3] != 0.
E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
E_therm_final = Granular.totalGrainThermalEnergy(sim)
diff --git a/test/collision-2floes-normal.jl b/test/collision-2floes-normal.jl
@@ -3,8 +3,6 @@
# Check for conservation of kinetic energy (=momentum) during a normal collision
# between two ice cylindrical grains
-Compat.@info "#### $(basename(@__FILE__)) ####"
-
verbose=false
Compat.@info "# One ice floe fixed"
diff --git a/test/collision-2floes-oblique.jl b/test/collision-2floes-oblique.jl
@@ -3,8 +3,6 @@
# Check for conservation of kinetic energy (=momentum) during a normal collision
# between two ice cylindrical grains
-Compat.@info "#### $(basename(@__FILE__)) ####"
-
verbose=false
Compat.@info "## Contact-normal elasticity only"
@@ -36,9 +34,6 @@ Granular.run!(sim, temporal_integration_method="Two-term Taylor", verbose=verbos
E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
E_thermal_final = Granular.totalGrainThermalEnergy(sim)
-println(E_kin_lin_init)
-println(E_kin_lin_final)
-println(E_thermal_final)
@test E_kin_lin_init ≈ E_kin_lin_final+E_thermal_final atol=E_kin_lin_init*tol
@test E_kin_rot_init ≈ E_kin_rot_final
@@ -136,19 +131,25 @@ sim_init.grains[2].contact_dynamic_friction = 1e2 # no Coulomb sliding
sim_init.grains[2].fixed = true
sim = deepcopy(sim_init)
+
Compat.@info "Testing kinetic energy conservation with Two-term Taylor scheme"
Granular.setTimeStep!(sim, epsilon=0.07)
tol = 0.1
Compat.@info "Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)"
+Granular.setOutputFileInterval!(sim, 1.0)
Granular.run!(sim, temporal_integration_method="Two-term Taylor",
verbose=verbose)
-@test sim.grains[1].ang_pos < 0.
-@test sim.grains[1].ang_vel < 0.
-@test sim.grains[2].ang_pos ≈ 0.
-@test sim.grains[2].ang_vel ≈ 0.
+@test sim.grains[1].ang_pos[3] < 0.
+@test sim.grains[1].ang_vel[3] < 0.
+@test sim.grains[2].ang_pos[3] ≈ 0.
+@test sim.grains[2].ang_vel[3] ≈ 0.
E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
+println(E_kin_lin_init)
+println(E_kin_lin_final)
+println(E_kin_rot_init)
+println(E_kin_rot_final)
@test E_kin_lin_init+E_kin_rot_init ≈ E_kin_lin_final+E_kin_rot_final atol=E_kin_lin_init*tol
Compat.@info "mu_d = 0."
@@ -162,9 +163,9 @@ E_kin_rot_init = Granular.totalGrainKineticRotationalEnergy(sim)
Granular.run!(sim, temporal_integration_method="Three-term Taylor",
verbose=verbose)
-@test sim.grains[1].ang_pos ≈ 0.
-@test sim.grains[1].ang_vel ≈ 0.
-@test sim.grains[2].ang_pos ≈ 0.
+@test sim.grains[1].ang_pos[3] ≈ 0.
+@test sim.grains[1].ang_vel[3] ≈ 0.
+@test sim.grains[2].ang_pos[3] ≈ 0.
E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
@test E_kin_lin_init ≈ E_kin_lin_final atol=E_kin_lin_init*tol
@@ -178,10 +179,10 @@ Compat.@info "Relative tolerance: $(tol*100.)%"
Granular.run!(sim, temporal_integration_method="Two-term Taylor",
verbose=verbose)
-@test sim.grains[1].ang_pos < 0.
-@test sim.grains[1].ang_vel < 0.
-@test sim.grains[2].ang_pos ≈ 0.
-@test sim.grains[2].ang_vel ≈ 0.
+@test sim.grains[1].ang_pos[3] < 0.
+@test sim.grains[1].ang_vel[3] < 0.
+@test sim.grains[2].ang_pos[3] ≈ 0.
+@test sim.grains[2].ang_vel[3] ≈ 0.
E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
@test E_kin_lin_init+E_kin_rot_init ≈ E_kin_lin_final+E_kin_rot_final atol=E_kin_lin_init*tol
@@ -195,10 +196,10 @@ Compat.@info "Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)"
Granular.run!(sim, temporal_integration_method="Three-term Taylor",
verbose=verbose)
-@test sim.grains[1].ang_pos < 0.
-@test sim.grains[1].ang_vel < 0.
-@test sim.grains[2].ang_pos ≈ 0.
-@test sim.grains[2].ang_vel ≈ 0.
+@test sim.grains[1].ang_pos[3] < 0.
+@test sim.grains[1].ang_vel[3] < 0.
+@test sim.grains[2].ang_pos[3] ≈ 0.
+@test sim.grains[2].ang_vel[3] ≈ 0.
E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
@test E_kin_lin_init+E_kin_rot_init ≈ E_kin_lin_final+E_kin_rot_final atol=E_kin_lin_init*tol
@@ -228,10 +229,10 @@ Compat.@info "Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)"
Granular.run!(sim, temporal_integration_method="Two-term Taylor",
verbose=verbose)
-@test sim.grains[1].ang_pos < 0.
-@test sim.grains[1].ang_vel < 0.
-@test sim.grains[2].ang_pos < 0.
-@test sim.grains[2].ang_vel < 0.
+@test sim.grains[1].ang_pos[3] < 0.
+@test sim.grains[1].ang_vel[3] < 0.
+@test sim.grains[2].ang_pos[3] < 0.
+@test sim.grains[2].ang_vel[3] < 0.
E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
@test E_kin_lin_init+E_kin_rot_init ≈ E_kin_lin_final+E_kin_rot_final atol=E_kin_lin_init*tol
@@ -257,10 +258,10 @@ Compat.@info "Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)"
Granular.run!(sim, temporal_integration_method="Three-term Taylor",
verbose=verbose)
-@test sim.grains[1].ang_pos < 0.
-@test sim.grains[1].ang_vel < 0.
-@test sim.grains[2].ang_pos < 0.
-@test sim.grains[2].ang_vel < 0.
+@test sim.grains[1].ang_pos[3] < 0.
+@test sim.grains[1].ang_vel[3] < 0.
+@test sim.grains[2].ang_pos[3] < 0.
+@test sim.grains[2].ang_vel[3] < 0.
E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
@test E_kin_lin_init+E_kin_rot_init ≈ E_kin_lin_final+E_kin_rot_final atol=E_kin_lin_init*tol
@@ -291,10 +292,10 @@ Compat.@info "Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)"
Granular.run!(sim, temporal_integration_method="Two-term Taylor",
verbose=verbose)
-@test sim.grains[1].ang_pos > 0.
-@test sim.grains[1].ang_vel > 0.
-@test sim.grains[2].ang_pos > 0.
-@test sim.grains[2].ang_vel > 0.
+@test sim.grains[1].ang_pos[3] > 0.
+@test sim.grains[1].ang_vel[3] > 0.
+@test sim.grains[2].ang_pos[3] > 0.
+@test sim.grains[2].ang_vel[3] > 0.
E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
@test E_kin_lin_init+E_kin_rot_init ≈ E_kin_lin_final+E_kin_rot_final atol=E_kin_lin_init*tol
@@ -320,10 +321,10 @@ Compat.@info "Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)"
Granular.run!(sim, temporal_integration_method="Three-term Taylor",
verbose=verbose)
-@test sim.grains[1].ang_pos > 0.
-@test sim.grains[1].ang_vel > 0.
-@test sim.grains[2].ang_pos > 0.
-@test sim.grains[2].ang_vel > 0.
+@test sim.grains[1].ang_pos[3] > 0.
+@test sim.grains[1].ang_vel[3] > 0.
+@test sim.grains[2].ang_pos[3] > 0.
+@test sim.grains[2].ang_vel[3] > 0.
E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
@test E_kin_lin_init+E_kin_rot_init ≈ E_kin_lin_final+E_kin_rot_final atol=E_kin_lin_init*tol
@@ -352,10 +353,10 @@ Compat.@info "Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)"
Granular.run!(sim, temporal_integration_method="Two-term Taylor",
verbose=verbose)
-@test sim.grains[1].ang_pos < 0.
-@test sim.grains[1].ang_vel < 0.
-@test sim.grains[2].ang_pos < 0.
-@test sim.grains[2].ang_vel < 0.
+@test sim.grains[1].ang_pos[3] < 0.
+@test sim.grains[1].ang_vel[3] < 0.
+@test sim.grains[2].ang_pos[3] < 0.
+@test sim.grains[2].ang_vel[3] < 0.
E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
@test E_kin_lin_init+E_kin_rot_init ≈ E_kin_lin_final+E_kin_rot_final atol=E_kin_lin_init*tol
@@ -381,10 +382,10 @@ Compat.@info "Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)"
Granular.run!(sim, temporal_integration_method="Three-term Taylor",
verbose=verbose)
-@test sim.grains[1].ang_pos < 0.
-@test sim.grains[1].ang_vel < 0.
-@test sim.grains[2].ang_pos < 0.
-@test sim.grains[2].ang_vel < 0.
+@test sim.grains[1].ang_pos[3] < 0.
+@test sim.grains[1].ang_vel[3] < 0.
+@test sim.grains[2].ang_pos[3] < 0.
+@test sim.grains[2].ang_vel[3] < 0.
E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
@test E_kin_lin_init+E_kin_rot_init ≈ E_kin_lin_final+E_kin_rot_final atol=E_kin_lin_init*tol
@@ -421,10 +422,10 @@ Compat.@info "Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)"
Granular.run!(sim, temporal_integration_method="Two-term Taylor",
verbose=verbose)
-@test sim.grains[1].ang_pos < 0.
-@test sim.grains[1].ang_vel < 0.
-@test sim.grains[2].ang_pos < 0.
-@test sim.grains[2].ang_vel < 0.
+@test sim.grains[1].ang_pos[3] < 0.
+@test sim.grains[1].ang_vel[3] < 0.
+@test sim.grains[2].ang_pos[3] < 0.
+@test sim.grains[2].ang_vel[3] < 0.
E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
@test E_kin_lin_init+E_kin_rot_init ≈ E_kin_lin_final+E_kin_rot_final atol=E_kin_lin_init*tol
@@ -450,10 +451,10 @@ Compat.@info "Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)"
Granular.run!(sim, temporal_integration_method="Three-term Taylor",
verbose=verbose)
-@test sim.grains[1].ang_pos < 0.
-@test sim.grains[1].ang_vel < 0.
-@test sim.grains[2].ang_pos < 0.
-@test sim.grains[2].ang_vel < 0.
+@test sim.grains[1].ang_pos[3] < 0.
+@test sim.grains[1].ang_vel[3] < 0.
+@test sim.grains[2].ang_pos[3] < 0.
+@test sim.grains[2].ang_vel[3] < 0.
E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
@test E_kin_lin_init+E_kin_rot_init ≈ E_kin_lin_final+E_kin_rot_final atol=E_kin_lin_init*tol
@@ -503,10 +504,10 @@ Compat.@info "Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)"
Granular.run!(sim, temporal_integration_method="Three-term Taylor",
verbose=verbose)
-@test sim.grains[1].ang_pos < 0.
-@test sim.grains[1].ang_vel < 0.
-@test sim.grains[2].ang_pos < 0.
-@test sim.grains[2].ang_vel < 0.
+@test sim.grains[1].ang_pos[3] < 0.
+@test sim.grains[1].ang_vel[3] < 0.
+@test sim.grains[2].ang_pos[3] < 0.
+@test sim.grains[2].ang_vel[3] < 0.
E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
@test E_kin_lin_init+E_kin_rot_init > E_kin_lin_final+E_kin_rot_final
@@ -556,10 +557,10 @@ Compat.@info "Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)"
Granular.run!(sim, temporal_integration_method="Three-term Taylor",
verbose=verbose)
-@test sim.grains[1].ang_pos < 0.
-@test sim.grains[1].ang_vel < 0.
-@test sim.grains[2].ang_pos < 0.
-@test sim.grains[2].ang_vel < 0.
+@test sim.grains[1].ang_pos[3] < 0.
+@test sim.grains[1].ang_vel[3] < 0.
+@test sim.grains[2].ang_pos[3] < 0.
+@test sim.grains[2].ang_vel[3] < 0.
E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
@test E_kin_lin_init+E_kin_rot_init > E_kin_lin_final+E_kin_rot_final
@@ -609,10 +610,10 @@ Compat.@info "Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)"
Granular.run!(sim, temporal_integration_method="Three-term Taylor",
verbose=verbose)
-@test sim.grains[1].ang_pos < 0.
-@test sim.grains[1].ang_vel < 0.
-@test sim.grains[2].ang_pos < 0.
-@test sim.grains[2].ang_vel < 0.
+@test sim.grains[1].ang_pos[3] < 0.
+@test sim.grains[1].ang_vel[3] < 0.
+@test sim.grains[2].ang_pos[3] < 0.
+@test sim.grains[2].ang_vel[3] < 0.
E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
@test E_kin_lin_init+E_kin_rot_init > E_kin_lin_final+E_kin_rot_final
diff --git a/test/collision-5floes-normal.jl b/test/collision-5floes-normal.jl
@@ -4,8 +4,6 @@ using Compat.LinearAlgebra
# Check for conservation of kinetic energy (=momentum) during a normal collision
# between two ice cylindrical grains
-Compat.@info "#### $(basename(@__FILE__)) ####"
-
verbose=false
Compat.@info "# One ice floe fixed"
diff --git a/test/compressive_failure.jl b/test/compressive_failure.jl
@@ -0,0 +1,184 @@
+#!/usr/bin/env julia
+using Compat.Test
+import Granular
+
+verbose = false
+debug = false
+if debug
+ import PyPlot
+end
+
+function plot_interaction(sim::Granular.Simulation, output::String)
+ time = Float64[]
+ force_n_1 = Float64[]
+ force_n_2 = Float64[]
+ force_t_1 = Float64[]
+ force_t_2 = Float64[]
+ torque_1 = Float64[]
+ torque_2 = Float64[]
+ compressive_failure = Bool[]
+ while sim.time < sim.time_total
+ Granular.run!(sim, verbose=verbose, single_step=true)
+ append!(time, sim.time)
+ append!(compressive_failure, sim.grains[1].compressive_failure[1])
+ append!(force_n_1, sim.grains[1].force[1])
+ append!(force_n_2, sim.grains[2].force[1])
+ append!(force_t_1, sim.grains[1].force[2])
+ append!(force_t_2, sim.grains[2].force[2])
+ append!(torque_1, sim.grains[1].torque[3])
+ append!(torque_2, sim.grains[2].torque[3])
+ end
+ PyPlot.clf()
+ PyPlot.subplot(3,1,1)
+ PyPlot.plot(time, force_n_1, "-b", label="1")
+ PyPlot.plot(time, force_n_2, "--y", label="2")
+ PyPlot.legend(loc="upper right")
+ PyPlot.ylabel("\$f_x\$ [N]")
+ PyPlot.subplot(3,1,2)
+ PyPlot.plot(time, force_t_1, "-b", label="1")
+ PyPlot.plot(time, force_t_2, "--y", label="2")
+ PyPlot.legend(loc="upper right")
+ PyPlot.ylabel("\$f_y\$ [N]")
+ PyPlot.subplot(3,1,3)
+ PyPlot.plot(time, torque_1, "-b", label="1")
+ PyPlot.plot(time, torque_2, "--y", label="2")
+ PyPlot.legend(loc="upper right")
+ PyPlot.ylabel("Torque [Nm]")
+ PyPlot.xlabel("Time [s]")
+ PyPlot.tight_layout()
+ PyPlot.savefig(output)
+end
+
+Compat.@info "Testing compressive failure: uniaxial compression"
+sim = Granular.createSimulation("compressive_failure_uniaxial")
+Granular.addGrainCylindrical!(sim, [0.,0.], 1., 0.5,
+ fracture_toughness=1285e3,
+ lin_vel=[1., 0.], fixed=true, verbose=verbose)
+Granular.addGrainCylindrical!(sim, [2.,0.], 1., 0.5,
+ fracture_toughness=1285e3,
+ fixed=true, verbose=verbose)
+@test count(x->x==true, sim.grains[1].compressive_failure) == 0
+Granular.setTimeStep!(sim, verbose=verbose)
+Granular.setTotalTime!(sim, 1.0)
+
+if debug
+ Granular.removeSimulationFiles(sim)
+ Granular.setOutputFileInterval!(sim, 0.01)
+ plot_interaction(sim, sim.id * ".pdf")
+else
+ Granular.run!(sim, verbose=verbose)
+end
+
+@test sim.grains[1].compressive_failure[1] == true
+@test count(x->x==true, sim.grains[1].compressive_failure) == 1
+@test sim.grains[1].force[1] < 0.0
+@test sim.grains[1].force[2] ≈ 0.0
+@test sim.grains[2].force[1] > 0.0
+@test sim.grains[2].force[2] ≈ 0.0
+@test sim.grains[1].torque ≈ zeros(3)
+@test sim.grains[2].torque ≈ zeros(3)
+
+Compat.@info "Testing compressive failure: shear"
+sim = Granular.createSimulation("compressive_failure_shear")
+Granular.addGrainCylindrical!(sim, [0.,0.], 1., 0.5,
+ fracture_toughness=1285e3,
+ lin_vel=[0., 1.], fixed=true, verbose=verbose)
+Granular.addGrainCylindrical!(sim, [1.5,1.5], 1., 0.5,
+ fracture_toughness=1285e3,
+ fixed=true, verbose=verbose)
+@test count(x->x==true, sim.grains[1].compressive_failure) == 0
+Granular.setTimeStep!(sim, verbose=verbose)
+Granular.setTotalTime!(sim, 1.0)
+
+if debug
+ Granular.removeSimulationFiles(sim)
+ Granular.setOutputFileInterval!(sim, 0.01)
+ plot_interaction(sim, sim.id * ".pdf")
+else
+ Granular.run!(sim, verbose=verbose)
+end
+
+@test sim.grains[1].compressive_failure[1] == true
+@test count(x->x==true, sim.grains[1].compressive_failure) == 1
+@test sim.grains[1].force[1] > 0.0
+@test sim.grains[1].force[2] < 0.0
+@test abs(sim.grains[1].force[1]) < abs(sim.grains[1].force[2])
+@test sim.grains[2].force[1] < 0.0
+@test sim.grains[2].force[2] > 0.0
+@test abs(sim.grains[2].force[1]) < abs(sim.grains[2].force[2])
+@test sim.grains[1].torque[1:2] ≈ zeros(2)
+@test sim.grains[1].torque[3] < 0.0
+@test sim.grains[2].torque[1:2] ≈ zeros(2)
+@test sim.grains[2].torque[3] < 0.0
+
+Compat.@info "Testing robustness of overlap calculations"
+sim = Granular.createSimulation("overlap")
+Granular.addGrainCylindrical!(sim, [0.,0.], 1., 0.5,
+ fracture_toughness=1285e3,
+ lin_vel=[0., 1.], fixed=true, verbose=verbose)
+Granular.addGrainCylindrical!(sim, [2.,0.], 1., 0.5,
+ fracture_toughness=1285e3,
+ fixed=true, verbose=verbose)
+@test count(x->x==true, sim.grains[1].compressive_failure) == 0
+Granular.setTimeStep!(sim, verbose=verbose)
+Granular.setTotalTime!(sim, 1.0)
+Granular.run!(sim, single_step=true, verbose=verbose)
+@test sim.grains[1].compressive_failure[1] == false
+@test sim.grains[1].contact_area[1] == 0.0
+@test count(x->x==true, sim.grains[1].compressive_failure) == 0
+
+sim = Granular.createSimulation("overlap")
+Granular.addGrainCylindrical!(sim, [0.,0.], 1., 0.5,
+ fracture_toughness=1.,
+ fixed=true, verbose=verbose)
+Granular.addGrainCylindrical!(sim, [0.0+1e-9,0.], 1., 0.5,
+ fracture_toughness=1.,
+ fixed=true, verbose=verbose)
+Granular.setTimeStep!(sim, verbose=verbose)
+Granular.setTotalTime!(sim, 1.0)
+Granular.run!(sim, single_step=true, verbose=verbose)
+@test sim.grains[1].compressive_failure[1] == true
+@test sim.grains[1].contact_area[1] ≈ π*1.0^2
+
+sim = Granular.createSimulation("overlap")
+Granular.addGrainCylindrical!(sim, [0.,0.], 1., 0.5,
+ fracture_toughness=1.,
+ fixed=true, verbose=verbose)
+Granular.addGrainCylindrical!(sim, [0.1,0.], 1., 0.5,
+ fracture_toughness=1.,
+ fixed=true, verbose=verbose)
+Granular.setTimeStep!(sim, verbose=verbose)
+Granular.setTotalTime!(sim, 1.0)
+Granular.run!(sim, single_step=true, verbose=verbose)
+@test sim.grains[1].compressive_failure[1] == true
+@test sim.grains[1].contact_area[1] < π*1.0^2
+@test sim.grains[1].contact_area[1] > 0.
+
+sim = Granular.createSimulation("overlap")
+Granular.addGrainCylindrical!(sim, [0.,0.], 1., 0.5,
+ fracture_toughness=1.,
+ fixed=true, verbose=verbose)
+Granular.addGrainCylindrical!(sim, [0.+1e-9,0.], 0.1, 0.5,
+ fracture_toughness=1.,
+ fixed=true, verbose=verbose)
+@test count(x->x==true, sim.grains[1].compressive_failure) == 0
+Granular.setTimeStep!(sim, verbose=verbose)
+Granular.setTotalTime!(sim, 1.0)
+Granular.run!(sim, single_step=true, verbose=verbose)
+@test sim.grains[1].position_vector[1] ≈ [-1e-9, 0., 0.]
+@test sim.grains[1].compressive_failure[1] == true
+@test sim.grains[1].contact_area[1] ≈ π*0.1^2
+
+sim = Granular.createSimulation("overlap")
+Granular.addGrainCylindrical!(sim, [0.,0.], 1., 0.5,
+ fracture_toughness=1.,
+ fixed=true, verbose=verbose)
+Granular.addGrainCylindrical!(sim, [0.3,0.4], 0.1, 0.5,
+ fracture_toughness=1.,
+ fixed=true, verbose=verbose)
+@test count(x->x==true, sim.grains[1].compressive_failure) == 0
+Granular.setTimeStep!(sim, verbose=verbose)
+Granular.setTotalTime!(sim, 1.0)
+Granular.run!(sim, single_step=true, verbose=verbose)
+@test sim.grains[1].compressive_failure[1] == true
+@test sim.grains[1].contact_area[1] ≈ π*0.1^2
diff --git a/test/contact-search-and-geometry.jl b/test/contact-search-and-geometry.jl
@@ -4,8 +4,6 @@ import Granular
# Check the contact search and geometry of a two-particle interaction
-Compat.@info "#### $(basename(@__FILE__)) ####"
-
Compat.@info "Testing interGrainPositionVector(...) and findOverlap(...)"
sim = Granular.createSimulation("test")
sim = Granular.createSimulation(id="test")
@@ -15,7 +13,7 @@ Granular.addGrainCylindrical!(sim, [18.01, 0.01], 10., 1., verbose=false)
position_ij = Granular.interGrainPositionVector(sim, 1, 2)
overlap_ij = Granular.findOverlap(sim, 1, 2, position_ij)
-@test [-18., 0.] ≈ position_ij
+@test [-18., 0., 0.] ≈ position_ij
@test -2. ≈ overlap_ij
@@ -31,16 +29,16 @@ Granular.findContacts!(sim)
sim.grains[1].enabled = false
# The contact should be registered in ice floe 1, but not ice floe 2
@test 2 == sim.grains[1].contacts[1]
-@test [-18., 0.] ≈ sim.grains[1].position_vector[1]
+@test [-18., 0., 0.] ≈ sim.grains[1].position_vector[1]
for ic=2:sim.Nc_max
@test 0 == sim.grains[1].contacts[ic]
- @test [0., 0.] ≈ sim.grains[1].position_vector[ic]
- @test [0., 0.] ≈ sim.grains[1].contact_parallel_displacement[ic]
+ @test [0., 0., 0.] ≈ sim.grains[1].position_vector[ic]
+ @test [0., 0., 0.] ≈ sim.grains[1].contact_parallel_displacement[ic]
end
for ic=1:sim.Nc_max
@test 0 == sim.grains[2].contacts[ic]
- @test [0., 0.] ≈ sim.grains[2].position_vector[ic]
- @test [0., 0.] ≈ sim.grains[2].contact_parallel_displacement[ic]
+ @test [0., 0., 0.] ≈ sim.grains[2].position_vector[ic]
+ @test [0., 0., 0.] ≈ sim.grains[2].contact_parallel_displacement[ic]
end
@test 1 == sim.grains[1].n_contacts
@test 1 == sim.grains[2].n_contacts
@@ -50,16 +48,16 @@ sim = deepcopy(sim_copy)
Granular.findContacts!(sim)
@test 2 == sim.grains[1].contacts[1]
-@test [-18., 0.] ≈ sim.grains[1].position_vector[1]
+@test [-18., 0., 0.] ≈ sim.grains[1].position_vector[1]
for ic=2:sim.Nc_max
@test 0 == sim.grains[1].contacts[ic]
- @test [0., 0.] ≈ sim.grains[1].position_vector[ic]
- @test [0., 0.] ≈ sim.grains[1].contact_parallel_displacement[ic]
+ @test [0., 0., 0.] ≈ sim.grains[1].position_vector[ic]
+ @test [0., 0., 0.] ≈ sim.grains[1].contact_parallel_displacement[ic]
end
for ic=1:sim.Nc_max
@test 0 == sim.grains[2].contacts[ic]
- @test [0., 0.] ≈ sim.grains[2].position_vector[ic]
- @test [0., 0.] ≈ sim.grains[2].contact_parallel_displacement[ic]
+ @test [0., 0., 0.] ≈ sim.grains[2].position_vector[ic]
+ @test [0., 0., 0.] ≈ sim.grains[2].contact_parallel_displacement[ic]
end
@test 1 == sim.grains[1].n_contacts
@test 1 == sim.grains[2].n_contacts
@@ -72,13 +70,13 @@ sim.grains[2].enabled = false
Granular.findContacts!(sim)
for ic=1:sim.Nc_max
@test 0 == sim.grains[1].contacts[ic]
- @test [0., 0.] ≈ sim.grains[1].position_vector[ic]
- @test [0., 0.] ≈ sim.grains[1].contact_parallel_displacement[ic]
+ @test [0., 0., 0.] ≈ sim.grains[1].position_vector[ic]
+ @test [0., 0., 0.] ≈ sim.grains[1].contact_parallel_displacement[ic]
end
for ic=1:sim.Nc_max
@test 0 == sim.grains[2].contacts[ic]
- @test [0., 0.] ≈ sim.grains[2].position_vector[ic]
- @test [0., 0.] ≈ sim.grains[2].contact_parallel_displacement[ic]
+ @test [0., 0., 0.] ≈ sim.grains[2].position_vector[ic]
+ @test [0., 0., 0.] ≈ sim.grains[2].contact_parallel_displacement[ic]
end
@test 0 == sim.grains[1].n_contacts
@test 0 == sim.grains[2].n_contacts
@@ -89,13 +87,13 @@ Granular.disableGrain!(sim, 1)
Granular.findContacts!(sim)
for ic=1:sim.Nc_max
@test 0 == sim.grains[1].contacts[ic]
- @test [0., 0.] ≈ sim.grains[1].position_vector[ic]
- @test [0., 0.] ≈ sim.grains[1].contact_parallel_displacement[ic]
+ @test [0., 0., 0.] ≈ sim.grains[1].position_vector[ic]
+ @test [0., 0., 0.] ≈ sim.grains[1].contact_parallel_displacement[ic]
end
for ic=1:sim.Nc_max
@test 0 == sim.grains[2].contacts[ic]
- @test [0., 0.] ≈ sim.grains[2].position_vector[ic]
- @test [0., 0.] ≈ sim.grains[2].contact_parallel_displacement[ic]
+ @test [0., 0., 0.] ≈ sim.grains[2].position_vector[ic]
+ @test [0., 0., 0.] ≈ sim.grains[2].contact_parallel_displacement[ic]
end
@test 0 == sim.grains[1].n_contacts
@test 0 == sim.grains[2].n_contacts
@@ -107,13 +105,13 @@ Granular.disableGrain!(sim, 2)
Granular.findContacts!(sim)
for ic=1:sim.Nc_max
@test 0 == sim.grains[1].contacts[ic]
- @test [0., 0.] ≈ sim.grains[1].position_vector[ic]
- @test [0., 0.] ≈ sim.grains[1].contact_parallel_displacement[ic]
+ @test [0., 0., 0.] ≈ sim.grains[1].position_vector[ic]
+ @test [0., 0., 0.] ≈ sim.grains[1].contact_parallel_displacement[ic]
end
for ic=1:sim.Nc_max
@test 0 == sim.grains[2].contacts[ic]
- @test [0., 0.] ≈ sim.grains[2].position_vector[ic]
- @test [0., 0.] ≈ sim.grains[2].contact_parallel_displacement[ic]
+ @test [0., 0., 0.] ≈ sim.grains[2].position_vector[ic]
+ @test [0., 0., 0.] ≈ sim.grains[2].contact_parallel_displacement[ic]
end
@test 0 == sim.grains[1].n_contacts
@test 0 == sim.grains[2].n_contacts
@@ -125,16 +123,16 @@ Granular.interact!(sim)
Granular.findContacts!(sim)
@test 2 == sim.grains[1].contacts[1]
-@test [-18., 0.] ≈ sim.grains[1].position_vector[1]
+@test [-18., 0., 0.] ≈ sim.grains[1].position_vector[1]
for ic=2:sim.Nc_max
@test 0 == sim.grains[1].contacts[ic]
- @test [0., 0.] ≈ sim.grains[1].position_vector[ic]
- @test [0., 0.] ≈ sim.grains[1].contact_parallel_displacement[ic]
+ @test [0., 0., 0.] ≈ sim.grains[1].position_vector[ic]
+ @test [0., 0., 0.] ≈ sim.grains[1].contact_parallel_displacement[ic]
end
for ic=1:sim.Nc_max
@test 0 == sim.grains[2].contacts[ic]
- @test [0., 0.] ≈ sim.grains[2].position_vector[ic]
- @test [0., 0.] ≈ sim.grains[2].contact_parallel_displacement[ic]
+ @test [0., 0., 0.] ≈ sim.grains[2].position_vector[ic]
+ @test [0., 0., 0.] ≈ sim.grains[2].contact_parallel_displacement[ic]
end
@test 1 == sim.grains[1].n_contacts
@test 1 == sim.grains[2].n_contacts
@@ -149,13 +147,13 @@ Granular.findContactsInGrid!(sim, sim.ocean)
@test 2 == sim.grains[1].contacts[1]
for ic=2:sim.Nc_max
@test 0 == sim.grains[1].contacts[ic]
- @test [0., 0.] ≈ sim.grains[1].position_vector[ic]
- @test [0., 0.] ≈ sim.grains[1].contact_parallel_displacement[ic]
+ @test [0., 0., 0.] ≈ sim.grains[1].position_vector[ic]
+ @test [0., 0., 0.] ≈ sim.grains[1].contact_parallel_displacement[ic]
end
for ic=1:sim.Nc_max
@test 0 == sim.grains[2].contacts[ic]
- @test [0., 0.] ≈ sim.grains[2].position_vector[ic]
- @test [0., 0.] ≈ sim.grains[2].contact_parallel_displacement[ic]
+ @test [0., 0., 0.] ≈ sim.grains[2].position_vector[ic]
+ @test [0., 0., 0.] ≈ sim.grains[2].contact_parallel_displacement[ic]
end
@test 1 == sim.grains[1].n_contacts
@test 1 == sim.grains[2].n_contacts
@@ -170,13 +168,13 @@ Granular.findContactsInGrid!(sim, sim.ocean)
@test 2 == sim.grains[1].contacts[1]
for ic=2:sim.Nc_max
@test 0 == sim.grains[1].contacts[ic]
- @test [0., 0.] ≈ sim.grains[1].position_vector[ic]
- @test [0., 0.] ≈ sim.grains[1].contact_parallel_displacement[ic]
+ @test [0., 0., 0.] ≈ sim.grains[1].position_vector[ic]
+ @test [0., 0., 0.] ≈ sim.grains[1].contact_parallel_displacement[ic]
end
for ic=1:sim.Nc_max
@test 0 == sim.grains[2].contacts[ic]
- @test [0., 0.] ≈ sim.grains[2].position_vector[ic]
- @test [0., 0.] ≈ sim.grains[2].contact_parallel_displacement[ic]
+ @test [0., 0., 0.] ≈ sim.grains[2].position_vector[ic]
+ @test [0., 0., 0.] ≈ sim.grains[2].contact_parallel_displacement[ic]
end
@test 1 == sim.grains[1].n_contacts
@test 1 == sim.grains[2].n_contacts
@@ -191,13 +189,13 @@ Granular.findContactsInGrid!(sim, sim.ocean)
for ic=1:sim.Nc_max
@test 0 == sim.grains[1].contacts[ic]
- @test [0., 0.] ≈ sim.grains[1].position_vector[ic]
- @test [0., 0.] ≈ sim.grains[1].contact_parallel_displacement[ic]
+ @test [0., 0., 0.] ≈ sim.grains[1].position_vector[ic]
+ @test [0., 0., 0.] ≈ sim.grains[1].contact_parallel_displacement[ic]
end
for ic=1:sim.Nc_max
@test 0 == sim.grains[2].contacts[ic]
- @test [0., 0.] ≈ sim.grains[2].position_vector[ic]
- @test [0., 0.] ≈ sim.grains[2].contact_parallel_displacement[ic]
+ @test [0., 0., 0.] ≈ sim.grains[2].position_vector[ic]
+ @test [0., 0., 0.] ≈ sim.grains[2].contact_parallel_displacement[ic]
end
@test 0 == sim.grains[1].n_contacts
@test 0 == sim.grains[2].n_contacts
@@ -211,13 +209,13 @@ Granular.findContacts!(sim)
@test 2 == sim.grains[1].contacts[1]
for ic=2:sim.Nc_max
@test 0 == sim.grains[1].contacts[ic]
- @test [0., 0.] ≈ sim.grains[1].position_vector[ic]
- @test [0., 0.] ≈ sim.grains[1].contact_parallel_displacement[ic]
+ @test [0., 0., 0.] ≈ sim.grains[1].position_vector[ic]
+ @test [0., 0., 0.] ≈ sim.grains[1].contact_parallel_displacement[ic]
end
for ic=1:sim.Nc_max
@test 0 == sim.grains[2].contacts[ic]
- @test [0., 0.] ≈ sim.grains[2].position_vector[ic]
- @test [0., 0.] ≈ sim.grains[2].contact_parallel_displacement[ic]
+ @test [0., 0., 0.] ≈ sim.grains[2].position_vector[ic]
+ @test [0., 0., 0.] ≈ sim.grains[2].contact_parallel_displacement[ic]
end
@test 1 == sim.grains[1].n_contacts
@test 1 == sim.grains[2].n_contacts
@@ -302,5 +300,3 @@ for j=2:(ny-1)
@test sim.grains[idx].n_contacts == 4
end
end
-
-
diff --git a/test/grain.jl b/test/grain.jl
@@ -2,22 +2,20 @@
# Check the basic icefloe functionality
-Compat.@info "#### $(basename(@__FILE__)) ####"
-
Compat.@info "Writing simple simulation to VTK file"
sim = Granular.createSimulation(id="test")
Granular.addGrainCylindrical!(sim, [ 0., 0.], 10., 1., verbose=false)
Granular.printGrainInfo(sim.grains[1])
Compat.@info "Testing grain value checks "
-@test_throws ErrorException Granular.addGrainCylindrical!(sim, [.1, .1, .1],
+@test_throws ErrorException Granular.addGrainCylindrical!(sim, [.1, .1, .1, .1],
10., 1.)
@test_throws ErrorException Granular.addGrainCylindrical!(sim, [.1, .1],
10., 1.,
- lin_vel=[.2,.2,.2])
+ lin_vel=[.2,.2,.2,.2])
@test_throws ErrorException Granular.addGrainCylindrical!(sim, [.1, .1],
10., 1.,
- lin_acc=[.2,.2,.2])
+ lin_acc=[.2,.2,.2,.2])
@test_throws ErrorException Granular.addGrainCylindrical!(sim, [.1, .1],
0., 1.)
@test_throws ErrorException Granular.addGrainCylindrical!(sim, [.1, .1],
@@ -82,18 +80,18 @@ end
Compat.@info "Testing external body force routines"
sim = Granular.createSimulation(id="test")
Granular.addGrainCylindrical!(sim, [ 0., 0.], 10., 1., verbose=false)
-Granular.setBodyForce!(sim.grains[1], [1., 2.])
-@test sim.grains[1].external_body_force ≈ [1., 2.]
-Granular.addBodyForce!(sim.grains[1], [1., 2.])
-@test sim.grains[1].external_body_force ≈ [2., 4.]
+Granular.setBodyForce!(sim.grains[1], [1., 2., 0.])
+@test sim.grains[1].external_body_force ≈ [1., 2., 0.]
+Granular.addBodyForce!(sim.grains[1], [1., 2., 0.])
+@test sim.grains[1].external_body_force ≈ [2., 4., 0.]
Compat.@info "Testing zeroKinematics!()"
-sim.grains[1].force .= ones(2)
-sim.grains[1].lin_acc .= ones(2)
-sim.grains[1].lin_vel .= ones(2)
-sim.grains[1].torque = 1.
-sim.grains[1].ang_acc = 1.
-sim.grains[1].ang_vel = 1.
+sim.grains[1].force .= ones(3)
+sim.grains[1].lin_acc .= ones(3)
+sim.grains[1].lin_vel .= ones(3)
+sim.grains[1].torque .= ones(3)
+sim.grains[1].ang_acc .= ones(3)
+sim.grains[1].ang_vel .= ones(3)
Granular.zeroKinematics!(sim)
@test Granular.totalGrainKineticTranslationalEnergy(sim) ≈ 0.
@test Granular.totalGrainKineticRotationalEnergy(sim) ≈ 0.
diff --git a/test/grid-boundaries.jl b/test/grid-boundaries.jl
@@ -1,8 +1,6 @@
#!/usr/bin/env julia
import Compat
-Compat.@info "#### $(basename(@__FILE__)) ####"
-
verbose=false
Compat.@info "## Inactive/Periodic BCs"
@@ -185,7 +183,7 @@ sim.ocean = Granular.createRegularOceanGrid([5, 5, 2], [1., 1., 1.])
Granular.setGridBoundaryConditions!(sim.ocean, "inactive", verbose=false)
Granular.addGrainCylindrical!(sim, [0.1, 0.5], 0.11, 0.1, verbose=false)
Granular.moveGrainsAcrossPeriodicBoundaries!(sim)
-@test [0.1, 0.5] ≈ sim.grains[1].lin_pos
+@test [0.1, 0.5, 0.] ≈ sim.grains[1].lin_pos
# do not readjust inside grid, periodic boundaries
sim = Granular.createSimulation()
@@ -193,7 +191,7 @@ sim.ocean = Granular.createRegularOceanGrid([5, 5, 2], [1., 1., 1.])
Granular.setGridBoundaryConditions!(sim.ocean, "periodic", verbose=false)
Granular.addGrainCylindrical!(sim, [0.1, 0.5], 0.11, 0.1, verbose=false)
Granular.moveGrainsAcrossPeriodicBoundaries!(sim)
-@test [0.1, 0.5] ≈ sim.grains[1].lin_pos
+@test [0.1, 0.5, 0.] ≈ sim.grains[1].lin_pos
# do not readjust outside grid, inactive boundaries
sim = Granular.createSimulation()
@@ -201,7 +199,7 @@ sim.ocean = Granular.createRegularOceanGrid([5, 5, 2], [1., 1., 1.])
Granular.setGridBoundaryConditions!(sim.ocean, "inactive", verbose=false)
Granular.addGrainCylindrical!(sim, [-0.1, 0.5], 0.11, 0.1, verbose=false)
Granular.moveGrainsAcrossPeriodicBoundaries!(sim)
-@test [-0.1, 0.5] ≈ sim.grains[1].lin_pos
+@test [-0.1, 0.5, 0.] ≈ sim.grains[1].lin_pos
# readjust outside grid, periodic boundaries, -x
sim = Granular.createSimulation()
@@ -209,7 +207,7 @@ sim.ocean = Granular.createRegularOceanGrid([5, 5, 2], [1., 1., 1.])
Granular.setGridBoundaryConditions!(sim.ocean, "periodic", verbose=false)
Granular.addGrainCylindrical!(sim, [-0.1, 0.5], 0.11, 0.1, verbose=false)
Granular.moveGrainsAcrossPeriodicBoundaries!(sim)
-@test [0.9, 0.5] ≈ sim.grains[1].lin_pos
+@test [0.9, 0.5, 0.] ≈ sim.grains[1].lin_pos
# readjust outside grid, periodic boundaries, +x
sim = Granular.createSimulation()
@@ -217,7 +215,7 @@ sim.ocean = Granular.createRegularOceanGrid([5, 5, 2], [1., 1., 1.])
Granular.setGridBoundaryConditions!(sim.ocean, "periodic", verbose=false)
Granular.addGrainCylindrical!(sim, [1.1, 0.5], 0.11, 0.1, verbose=false)
Granular.moveGrainsAcrossPeriodicBoundaries!(sim)
-@test [0.1, 0.5] ≈ sim.grains[1].lin_pos
+@test [0.1, 0.5, 0.] ≈ sim.grains[1].lin_pos
# readjust outside grid, periodic boundaries, -y
sim = Granular.createSimulation()
@@ -225,7 +223,7 @@ sim.ocean = Granular.createRegularOceanGrid([5, 5, 2], [1., 1., 1.])
Granular.setGridBoundaryConditions!(sim.ocean, "periodic", verbose=false)
Granular.addGrainCylindrical!(sim, [0.3, -0.1], 0.11, 0.1, verbose=false)
Granular.moveGrainsAcrossPeriodicBoundaries!(sim)
-@test [0.3, 0.9] ≈ sim.grains[1].lin_pos
+@test [0.3, 0.9, 0.] ≈ sim.grains[1].lin_pos
# readjust outside grid, periodic boundaries, +y
sim = Granular.createSimulation()
@@ -233,7 +231,7 @@ sim.ocean = Granular.createRegularOceanGrid([5, 5, 2], [1., 1., 1.])
Granular.setGridBoundaryConditions!(sim.ocean, "periodic", verbose=false)
Granular.addGrainCylindrical!(sim, [0.3, 1.1], 0.11, 0.1, verbose=false)
Granular.moveGrainsAcrossPeriodicBoundaries!(sim)
-@test [0.3, 0.1] ≈ sim.grains[1].lin_pos
+@test [0.3, 0.1, 0.] ≈ sim.grains[1].lin_pos
# throw error if atmosphere and ocean BCs differ
sim = Granular.createSimulation()
@@ -253,4 +251,4 @@ sim.ocean = Granular.createRegularOceanGrid([5, 5, 2], [1., 1., 1.])
Granular.setGridBoundaryConditions!(sim.ocean, "inactive", verbose=false)
Granular.addGrainCylindrical!(sim, [0.1, 0.5], 0.11, 0.1, verbose=false)
Granular.moveGrainsAcrossPeriodicBoundaries!(sim)
-@test [0.1, 0.5] ≈ sim.grains[1].lin_pos
+@test [0.1, 0.5, 0.] ≈ sim.grains[1].lin_pos
diff --git a/test/grid.jl b/test/grid.jl
@@ -5,8 +5,6 @@ import Granular
# Check the grid interpolation and sorting functions
verbose = false
-Compat.@info "#### $(basename(@__FILE__)) ####"
-
if Granular.hasNetCDF
ocean = Granular.readOceanNetCDF("Baltic/00010101.ocean_month.nc",
"Baltic/ocean_hgrid.nc")
diff --git a/test/jld.jl b/test/jld.jl
@@ -1,8 +1,6 @@
#!/usr/bin/env julia
import Compat
-Compat.@info "#### $(basename(@__FILE__)) ####"
-
Compat.@info "Determining if JLD is installed"
if Granular.hasJLD
Compat.@info "JLD found, proceeding with JLD-specific tests"
diff --git a/test/netcdf.jl b/test/netcdf.jl
@@ -2,8 +2,6 @@
# Check if NetCDF files are read correctly from the disk.
-Compat.@info "#### $(basename(@__FILE__)) ####"
-
@test_throws ErrorException Granular.readOceanStateNetCDF("nonexistentfile")
@test_throws ErrorException Granular.readOceanGridNetCDF("nonexistentfile")
diff --git a/test/ocean.jl b/test/ocean.jl
@@ -3,8 +3,6 @@
# Check if ocean-specific functions and grid operations are functioning
# correctly
-Compat.@info "#### $(basename(@__FILE__)) ####"
-
Compat.@info "Testing regular grid generation"
sim = Granular.createSimulation()
sim.ocean = Granular.createRegularOceanGrid([6, 6, 6], [1., 1., 1.])
@@ -58,15 +56,15 @@ E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
Compat.@info "Testing vortex interaction (static ocean)"
sim = deepcopy(sim_init)
-sim.grains[1].ang_vel = 0.1
+sim.grains[1].ang_vel[3] = 0.1
E_kin_lin_init = Granular.totalGrainKineticTranslationalEnergy(sim)
E_kin_rot_init = Granular.totalGrainKineticRotationalEnergy(sim)
Granular.run!(sim, verbose=false)
E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
@test E_kin_rot_init > E_kin_rot_final # energy lost to ocean
-@test sim.grains[1].ang_vel > 0. # check angular velocity orientation
-@test sim.grains[1].ang_pos > 0. # check angular position orientation
+@test sim.grains[1].ang_vel[3] > 0. # check angular velocity orientation
+@test sim.grains[1].ang_pos[3] > 0. # check angular position orientation
@test E_kin_lin_init ≈ E_kin_lin_final # no linear velocity gained
Compat.@info "Testing vortex interaction (static ice floe)"
@@ -81,8 +79,8 @@ E_kin_rot_init = Granular.totalGrainKineticRotationalEnergy(sim)
Granular.run!(sim, verbose=false)
E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
-@test sim.grains[1].ang_vel > 0. # check angular velocity orientation
-@test sim.grains[1].ang_pos > 0. # check angular position orientation
+@test sim.grains[1].ang_vel[3] > 0. # check angular velocity orientation
+@test sim.grains[1].ang_pos[3] > 0. # check angular position orientation
@test E_kin_rot_init < E_kin_rot_final # rotation after due to ocean vortex
@test E_kin_lin_init ≈ E_kin_lin_final # no linear velocity gained
@@ -97,8 +95,8 @@ E_kin_rot_init = Granular.totalGrainKineticRotationalEnergy(sim)
Granular.run!(sim, verbose=false)
E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
-@test sim.grains[1].ang_vel < 0. # check angular velocity orientation
-@test sim.grains[1].ang_pos < 0. # check angular position orientation
+@test sim.grains[1].ang_vel[3] < 0. # check angular velocity orientation
+@test sim.grains[1].ang_pos[3] < 0. # check angular position orientation
@test E_kin_rot_init < E_kin_rot_final # rotation after due to ocean vortex
@test E_kin_lin_init ≈ E_kin_lin_final # no linear velocity gained
@@ -113,8 +111,8 @@ E_kin_rot_init = Granular.totalGrainKineticRotationalEnergy(sim)
Granular.run!(sim, verbose=false)
E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
-@test sim.grains[1].ang_vel < 0. # check angular velocity orientation
-@test sim.grains[1].ang_pos < 0. # check angular position orientation
+@test sim.grains[1].ang_vel[3] < 0. # check angular velocity orientation
+@test sim.grains[1].ang_pos[3] < 0. # check angular position orientation
@test E_kin_rot_init < E_kin_rot_final # rotation after due to ocean vortex
@test E_kin_lin_init ≈ E_kin_lin_final # no linear velocity gained
@@ -129,7 +127,7 @@ E_kin_rot_init = Granular.totalGrainKineticRotationalEnergy(sim)
Granular.run!(sim, verbose=false)
E_kin_lin_final = Granular.totalGrainKineticTranslationalEnergy(sim)
E_kin_rot_final = Granular.totalGrainKineticRotationalEnergy(sim)
-@test sim.grains[1].ang_vel > 0. # check angular velocity orientation
-@test sim.grains[1].ang_pos > 0. # check angular position orientation
+@test sim.grains[1].ang_vel[3] > 0. # check angular velocity orientation
+@test sim.grains[1].ang_pos[3] > 0. # check angular position orientation
@test E_kin_rot_init < E_kin_rot_final # rotation after due to ocean vortex
@test E_kin_lin_init ≈ E_kin_lin_final # no linear velocity gained
diff --git a/test/packing.jl b/test/packing.jl
@@ -6,8 +6,6 @@ verbose = false
plot = false
plot_packings=false
-Compat.@info "#### $(basename(@__FILE__)) ####"
-
Compat.@info "Testing regular packing generation (power law GSD)"
sim = Granular.createSimulation()
Granular.regularPacking!(sim, [2, 2], 1., 1., size_distribution="powerlaw")
diff --git a/test/profiling.jl b/test/profiling.jl
@@ -9,8 +9,6 @@ import Plots
import Granular
import CurveFit
-Compat.@info "#### $(basename(@__FILE__)) ####"
-
verbose=false
Compat.@info "Testing performance with many interacting grains"
diff --git a/test/runtests.jl b/test/runtests.jl
@@ -3,22 +3,28 @@ using Compat.Test
using Compat.LinearAlgebra
import Granular
-include("collision-2floes-normal.jl")
-include("collision-5floes-normal.jl")
-include("collision-2floes-oblique.jl")
-include("grid.jl")
-include("contact-search-and-geometry.jl")
-include("grid-boundaries.jl")
-include("ocean.jl")
-include("atmosphere.jl")
-include("wall.jl")
-include("grain.jl")
-include("packing.jl")
-include("util.jl")
-include("temporal.jl")
-include("cohesion.jl")
+function run_test(filename::String)
+ Compat.printstyled("Info: #### $filename ####\n", color=:green)
+ include(filename)
+end
+
+run_test("compressive_failure.jl")
+run_test("cohesion.jl")
+run_test("grain.jl")
+run_test("vtk.jl")
+run_test("collision-2floes-normal.jl")
+run_test("collision-5floes-normal.jl")
+run_test("collision-2floes-oblique.jl")
+run_test("grid.jl")
+run_test("contact-search-and-geometry.jl")
+run_test("grid-boundaries.jl")
+run_test("ocean.jl")
+run_test("atmosphere.jl")
+run_test("wall.jl")
+run_test("packing.jl")
+run_test("util.jl")
+run_test("temporal.jl")
if Granular.hasNetCDF
- include("netcdf.jl")
+ run_test("netcdf.jl")
end
-include("vtk.jl")
-include("jld.jl")
+run_test("jld.jl")
diff --git a/test/util.jl b/test/util.jl
@@ -1,8 +1,8 @@
#!/usr/bin/env julia
+import Granular
import Compat
using Compat.Random
-
-Compat.@info "#### $(basename(@__FILE__)) ####"
+using Compat.Test
Compat.@info "Testing power-law RNG"
@@ -26,3 +26,12 @@ for i=1:10^5
@test 0. <= minimum(Granular.randpower(5, 1., 0., 1.))
@test 1. >= minimum(Granular.randpower(5, 1., 0., 1.))
end
+
+@test [1,2,0] == Granular.vecTo3d([1,2])
+@test [1,2,3] == Granular.vecTo3d([1,2], fill=3)
+@test [1,3,3] == Granular.vecTo3d([1], fill=3)
+@test [1,3,3] == Granular.vecTo3d(1, fill=3)
+@test [1.,2.,3.] == Granular.vecTo3d([1.,2.], fill=3.)
+@test [1.,3.,3.] == Granular.vecTo3d(1., fill=3.)
+@test [1.,0.,0.] == Granular.vecTo3d(1.)
+@test [1.,0.,0.] == Granular.vecTo3d([1.])
diff --git a/test/vtk.jl b/test/vtk.jl
@@ -3,8 +3,6 @@ import Compat
# Check the contact search and geometry of a two-particle interaction
-Compat.@info "#### $(basename(@__FILE__)) ####"
-
Compat.@info "Writing simple simulation to VTK file"
sim = Granular.createSimulation(id="test")
Granular.addGrainCylindrical!(sim, [ 0., 0.], 10., 1., verbose=false)
@@ -29,12 +27,12 @@ end
grainpath = "test/test.grains.1.vtu"
grainchecksum =
-"eff130f14975dd5da06186c5b61c0c1a9d679f2188d03019c8144898ddc902b6 " *
+"b1748cb82e8270d951d9c1acaea0d151f0a5c48a1255e75de350bf9bcbc15fe9 " *
grainpath * "\n"
graininteractionpath = "test/test.grain-interaction.1.vtp"
graininteractionchecksum =
-"c61f314a997c4405eaa98e6a4e3f81368ab2e1c60d06d9a0d998b69d7c8fb1bf " *
+"b8e49252a0ac87c2fce05e68ffab46589853429dc9f75d89818e4a37b953b137 " *
graininteractionpath * "\n"
oceanpath = "test/test.ocean.1.vts"
diff --git a/test/wall.jl b/test/wall.jl
@@ -2,14 +2,12 @@
# Check the basic dynamic wall functionality
-Compat.@info "#### $(basename(@__FILE__)) ####"
-
Compat.@info "# Test wall initialization"
Compat.@info "Testing argument value checks"
sim = Granular.createSimulation()
Granular.addGrainCylindrical!(sim, [ 0., 0.], 10., 2., verbose=false)
@test_throws ErrorException Granular.addWallLinearFrictionless!(sim,
- [.1, .1, .1],
+ [.1, .1, .1, .1],
1.)
@test_throws ErrorException Granular.addWallLinearFrictionless!(sim,
[1., 1.],
@@ -18,7 +16,7 @@ Granular.addGrainCylindrical!(sim, [ 0., 0.], 10., 2., verbose=false)
[.1, .1, .1],
1.)
@test_throws ErrorException Granular.addWallLinearFrictionless!(sim,
- [0., 1.], 1.,
+ [0., 1., 1.], 1.,
bc="asdf")
sim = Granular.createSimulation()
@test_throws ErrorException Granular.addWallLinearFrictionless!(sim, [1., 0.],
@@ -56,7 +54,7 @@ sim.walls[1].force = 1.0
@test Granular.getWallNormalStress(sim, wall_index=1, stress_type="effective") ≈ 1.0/(20.0*2.0)
@test_throws ErrorException Granular.getWallNormalStress(sim, wall_index=1, stress_type="nonexistent")
-sim.walls[1].normal = [1.0, 1.0]
+sim.walls[1].normal = [1.0, 1.0, 1.0]
@test_throws ErrorException Granular.getWallSurfaceArea(sim, 1)
@test_throws ErrorException Granular.getWallSurfaceArea(sim, [1.,1.], 0.5)
