Granular.jl

commit ec8211d8369f687eb8c6ac01638559c5809c203e
parent afc9a66cbc465d81c7ec692a18b77dcb310b28a0
Author: Anders Damsgaard <andersd@riseup.net>
Date:   Fri,  2 Jun 2017 14:57:52 -0400

add more fields to interaction output, do not correct effective radius for overlap

Diffstat:
M
src/interaction.jl
 | 
2
+-
M
src/io.jl
 | 
44
++++++++++++++++++++++++++++++++++++++++++++

2 files changed, 45 insertions(+), 1 deletion(-)
diff --git a/src/interaction.jl b/src/interaction.jl
@@ -80,7 +80,7 @@ function interactIceFloes!(simulation::Simulation, i::Int, j::Int, ic::Int)
     δ_t = dot(t, δ_t0 - (n*dot(n, δ_t0))) + vel_t*simulation.time_step
 
     # Effective radius
-    R_ij = harmonicMean(r_i, r_j) - abs(δ_n)/2.
+    R_ij = harmonicMean(r_i, r_j)
 
     # Contact area
     A_ij = R_ij*min(simulation.ice_floes[i].thickness, 
diff --git a/src/io.jl b/src/io.jl
@@ -146,6 +146,10 @@ function writeIceFloeInteractionVTK(simulation::Simulation,
     i2 = Int64[]
     inter_particle_vector = Array{float, 1}[]
     force = float[]
+    effective_radius = float[]
+    contact_area = float[]
+    contact_stiffness = float[]
+    tensile_stress = float[]
     shear_displacement = Array{float, 1}[]
     contact_age = float[]
     for i=1:length(simulation.ice_floes)
@@ -189,6 +193,10 @@ function writeIceFloeInteractionVTK(simulation::Simulation,
                 push!(inter_particle_vector, p)
 
                 push!(force, k_n*δ_n)
+                push!(effective_radius, R_ij)
+                push!(contact_area, A_ij)
+                push!(contact_stiffness, k_n)
+                push!(tensile_stress, k_n*δ_n/A_ij)
 
                 push!(shear_displacement, simulation.ice_floes[i].
                       contact_parallel_displacement[ic])
@@ -246,6 +254,42 @@ function writeIceFloeInteractionVTK(simulation::Simulation,
         write(f, "        </DataArray>\n")
 
         write(f, "        <DataArray type=\"Float32\" " *
+              "Name=\"Effective radius [m]\" " *
+              "NumberOfComponents=\"1\" format=\"ascii\">\n")
+        for i=1:length(i1)
+            write(f, "$(effective_radius[i]) ")
+        end
+        write(f, "\n")
+        write(f, "        </DataArray>\n")
+
+        write(f, "        <DataArray type=\"Float32\" " *
+              "Name=\"Contact area [m^2]\" " *
+              "NumberOfComponents=\"1\" format=\"ascii\">\n")
+        for i=1:length(i1)
+            write(f, "$(contact_area[i]) ")
+        end
+        write(f, "\n")
+        write(f, "        </DataArray>\n")
+
+        write(f, "        <DataArray type=\"Float32\" " *
+              "Name=\"Contact stiffness [N/m]\" " *
+              "NumberOfComponents=\"1\" format=\"ascii\">\n")
+        for i=1:length(i1)
+            write(f, "$(contact_stiffness[i]) ")
+        end
+        write(f, "\n")
+        write(f, "        </DataArray>\n")
+
+        write(f, "        <DataArray type=\"Float32\" " *
+              "Name=\"Tensile stress [Pa]\" " *
+              "NumberOfComponents=\"1\" format=\"ascii\">\n")
+        for i=1:length(i1)
+            write(f, "$(tensile_stress[i]) ")
+        end
+        write(f, "\n")
+        write(f, "        </DataArray>\n")
+
+        write(f, "        <DataArray type=\"Float32\" " *
               "Name=\"Contact age [s]\" NumberOfComponents=\"1\" 
         format=\"ascii\">\n")
         for i=1:length(i1)