tminor changes. Walls added to deformation - granular-basin - tectonic deformation experiments with Granular.jl
 (HTM) git clone git://src.adamsgaard.dk/granular-basin
 (DIR) Log
 (DIR) Files
 (DIR) Refs
 (DIR) README
       ---
 (DIR) commit 0813325bb8a04723d8c26c333d18b44d185be2d5
 (DIR) parent dcec4954da477abbe21e251e1294523b68a6db85
 (HTM) Author: esbenpalmstrom <esbenpalmstroem@gmail.com>
       Date:   Thu, 25 Nov 2021 17:33:43 +0100
       
       minor changes. Walls added to deformation
       
       Diffstat:
         M compact_basin.jl                    |       8 ++++++--
         M deform_basin.jl                     |      58 ++++++++++++++++++++++++-------
         M init_basin.jl                       |       8 +++-----
         M layer_basin.jl                      |      35 +++++++++++++++++++++++++-------
       
       4 files changed, 83 insertions(+), 26 deletions(-)
       ---
 (DIR) diff --git a/compact_basin.jl b/compact_basin.jl
       t@@ -8,12 +8,11 @@ t_start = Dates.now() # Save the start time, print the end time later.
        # lav en lille test? se om dit appendede carpet stadig er forbundet til hoved-
        # simulationsobjektet
        
       -id = "simulation500"    # id of simulation to load
       +id = "simulation1000"    # id of simulation to load
        N = 20e3                # amount of stress to be applied
        t_comp = 3.0            # compaction max duration [s]
        
        sim = Granular.readSimulation("$(id)/init.jld2")
       -carpet = Granular.readSimulation("$(id)/carpet.jld2")
        SimSettings = SimSettings = JLD2.load("$(id)/SimSettings.jld2")
        
        #mkpath("$(id)/compaction-N$(N)Pa")
       t@@ -98,3 +97,8 @@ cd("..")
        JLD2.save("simulation$(ngrains)/SimSettings.jld2", SimSettings)
        
        Granular.writeSimulation(sim,filename = "$(id)/comp.jld2")
       +
       +# print time elapsed
       +t_now = Dates.now()
       +dur = Dates.canonicalize(t_now-t_start)
       +print("Time elapsed: ",dur)
 (DIR) diff --git a/deform_basin.jl b/deform_basin.jl
       t@@ -27,6 +27,13 @@ for grain in sim.grains
            grain.fixed = false
        end
        
       +y_bot = Inf
       +for grain in sim.grains
       +    if y_bot > grain.lin_pos[2] - grain.contact_radius
       +        global y_bot = grain.lin_pos[2] - grain.contact_radius
       +    end
       +end
       +
        # Add Indenter
        temp_indent = Granular.createSimulation("id=temp_indent")
        
       t@@ -64,7 +71,12 @@ Granular.fitGridToGrains!(sim,
        
        sim.time_iteration = 0
        sim.time = 0.0
       -sim.file_time_since_output_file = 0.
       +sim.file_time_since_output_file = 0.y_bot = Inf
       +for grain in sim.grains
       +    if y_bot > grain.lin_pos[2] - grain.contact_radius
       +        global y_bot = grain.lin_pos[2] - grain.contact_radius
       +    end
       +end
        Granular.setTotalTime!(sim,2.0)
        Granular.setTimeStep!(sim)
        Granular.setOutputFileInterval!(sim, .01)
       t@@ -72,24 +84,46 @@ Granular.resetTime!(sim)
        
        cd("$id")
        sim.id = "deformed"
       -sim.walls = Granular.WallLinear[] # remove existing walls
       +sim.walls = Granular.WallLinearFrictionless[] # remove existing walls
        
        
        #find the edge grains of the carpet
       +left_edge = -Inf
       +right_edge = Inf
       +for i = 1:size(sim.grains,1)
       +    if left_edge < sim.grains[i].lin_pos[1] + sim.grains[i].contact_radius
       +        global left_edge = sim.grains[i].lin_pos[1] + sim.grains[i].contact_radius
       +        left_edge_index = deepcopy(i)
       +    end
       +    if right_edge >sim.grains[i].lin_pos[1] - sim.grains[i].contact_radius
       +        global right_edge = sim.grains[i].lin_pos[1] - sim.grains[i].contact_radius
       +        right_edge_index = deepcopy(i)
       +    end
       +end
        
        
       -#add walls
       -#facing east
       -Granular.addWallLinearFrictionless(sim, [-1.,0.]
       -                                    )
       +#add walls to the east and west
       +Granular.addWallLinearFrictionless!(sim,[1.,0.],
       +                                    left_edge,
       +                                    bc = "fixed")
        
       -while sim.time < sim.time_total
       -    for grain in sim.grains
       +Granular.addWallLinearFrictionless!(sim,[1.,0.],
       +                                    right_edge,
       +                                    bc = "fixed")
        
       -        if grain.lin_vel[2] < 0 && grain.color == 1
       -            grain.lin_vel[2] = 0
       -        end
       -    end
       +#add wall beneath the carpet
       +
       +Granular.addWallLinearFrictionless!(sim, [0.,1.],
       +                                    y_bot,
       +                                    bc = "fixed")
       +
       +while sim.time < sim.time_total
       +#    for grain in sim.grains
       +#
       +#        if grain.lin_vel[2] < 0 && grain.color == 1
       +#            grain.lin_vel[2] = 0
       +#        end
       +#    end
            Granular.run!(sim,single_step = true)
        
        end
 (DIR) diff --git a/init_basin.jl b/init_basin.jl
       t@@ -204,6 +204,7 @@ for i = left_edge+(bot_r/2):bot_r*1.99:left_edge+length
                                        color = 1)
        end
        
       +Granular.findContactsAllToAll!(carpet) # find the grain contacts
        
        append!(sim.grains,carpet.grains) # add the carpet grains to the main simulation object
        # since the assignment will point to the carpet object, changes made to the carpet
       t@@ -214,6 +215,7 @@ for grain in sim.grains
            grain.n_contacts = 0
        end
        
       +
        for grain in sim.grains
                for ic=1:size(grain.contact_age,1)
                        grain.contact_age[ic] = 1e16
       t@@ -221,7 +223,7 @@ for grain in sim.grains
            grain.strength_heal_rate = 1 # new bond stengthening
        end
        
       -Granular.findContactsAllToAll!(carpet) # find the grain contacts
       +#Granular.findContactsAllToAll!(carpet) # find the grain contacts
        
        Granular.fitGridToGrains!(sim,sim.ocean,verbose=false)  # fit the ocean to the added grains
        
       t@@ -244,10 +246,6 @@ cd("..")
        Granular.writeSimulation(sim,
                                filename = "simulation$(ngrains)/init.jld2")
        
       -
       -Granular.writeSimulation(carpet,
       -                        filename = "simulation$(ngrains)/carpet.jld2")
       -
        JLD2.save("simulation$(ngrains)/SimSettings.jld2", SimSettings)
        
        # print time elapsed
 (DIR) diff --git a/layer_basin.jl b/layer_basin.jl
       t@@ -3,7 +3,7 @@ import JLD2
        import PyPlot
        import Dates
        
       -id = "simulation500"    # id of simulation to load, just write the folder
       +id = "simulation1000"    # id of simulation to load, just write the folder
                                # name here
        
        # Layer interface positions
       t@@ -16,8 +16,9 @@ interfaces = [0,0.4,0.6,1]
        # mechanical properties for each layer
        youngs_modulus = [2e7,2e7,2e7]              # elastic modulus
        poissons_ratio = [0.185,0.185,0.185]        # shear stiffness ratio
       -tensile_strength = [0.3,0.05,0.3]            # strength of bonds between grains
       -contact_dynamic_friction = [0.4,0.4,0.4]    # friction between grains
       +tensile_strength = [0.3,0.05,0.3]           # strength of bonds between grains
       +shear_strength = [0.3,0.05,0.3]             # shear stregth of bonds
       +contact_dynamic_friction = [0.4,0.05,0.4]   # friction between grains
        rotating = [true,true,true]                 # can grains rotate or not
        color = [0,0,0]
        
       t@@ -26,13 +27,14 @@ carpet_poissons_ratio = 0.185
        carpet_tensile_strength = Inf
        carpet_contact_dynamic_friction = 0.4
        carpet_rotating = true
       +carpet_shear_strength = Inf
        
        sim = Granular.readSimulation("$(id)/comp.jld2")
       -carpet = Granular.readSimulation("$(id)/carpet.jld2")
        SimSettings = SimSettings = JLD2.load("$(id)/SimSettings.jld2")
        
       +sim.walls = Granular.WallLinearFrictionless[] # remove existing walls
       +
        Granular.zeroKinematics!(sim)       # end any movement
       -Granular.zeroKinematics!(carpet)    # end any movement
        
        y_top = -Inf
        for grain in sim.grains
       t@@ -72,6 +74,7 @@ for grain in sim.grains
                    grain.youngs_modulus = youngs_modulus[i-1]
                    grain.poissons_ratio = poissons_ratio[i-1]
                    grain.tensile_strength = tensile_strength[i-1]
       +            grain.shear_strength = shear_strength[i-1]
                    grain.contact_dynamic_friction = contact_dynamic_friction[i-1]
                    grain.rotating = rotating[i-1]
                    grain.color = color[i-1]
       t@@ -79,6 +82,7 @@ for grain in sim.grains
                    grain.youngs_modulus = carpet_youngs_modulus
                    grain.poissons_ratio = carpet_poissons_ratio
                    grain.tensile_strength = carpet_tensile_strength
       +            grain.shear_strength = carpet_shear_strength
                    grain.contact_dynamic_friction = carpet_contact_dynamic_friction
                    grain.rotating = carpet_rotating
                end
       t@@ -103,18 +107,35 @@ for grain in sim.grains
            end
        end
        
       -Granular.findContactsAllToAll!(sim) # find the grain contacts
       +#Granular.findContactsAllToAll!(sim) # find the grain contacts
       +#Granular.run!(sim,single_step=true)
        
        #reduce the contact radius again
        for i = 1:size(sim.grains,1)
            sim.grains[i].contact_radius -= increase_array[i]
        end
        
       +"""
       +for grain in sim.grains
       +    grain.contacts[:] .= 0
       +    grain.n_contacts = 0
       +end
       +
       +#vil det ikke være nødvendigt at køre et enkelt timestep her?
       +
       +for grain in sim.grains
       +        for ic=1:size(grain.contact_age,1)
       +                grain.contact_age[ic] = 1e16
       +        end
       +    grain.strength_heal_rate = 1 # new bond stengthening
       +end
       +"""
       +
        cd("$id")
        sim.id = "layered"
        
        Granular.resetTime!(sim)
       -Granular.setTotalTime!(sim,0.2)
       +Granular.setTotalTime!(sim,1.0)
        Granular.run!(sim)
        
        cd("..")