tadd arrays for twisting and bending of the bond - slidergrid - grid of elastic sliders on a frictional surface
 (HTM) git clone git://src.adamsgaard.dk/slidergrid
 (DIR) Log
 (DIR) Files
 (DIR) Refs
 (DIR) README
 (DIR) LICENSE
       ---
 (DIR) commit d852181baa63f618bea35aa6d8d6053763688920
 (DIR) parent 86fc15639894a4f2ebaaba73678afc6245e2d8e7
 (HTM) Author: Anders Damsgaard Christensen <adc@geo.au.dk>
       Date:   Mon, 18 Apr 2016 14:55:29 -0700
       
       add arrays for twisting and bending of the bond
       
       Diffstat:
         M slidergrid/slider.c                 |      79 +++++++++++++++++++++++++++++--
         M slidergrid/slider.h                 |       6 +++++-
         M tests/elasticity/Makefile           |      10 +++++++++-
       
       3 files changed, 90 insertions(+), 5 deletions(-)
       ---
 (DIR) diff --git a/slidergrid/slider.c b/slidergrid/slider.c
       t@@ -45,6 +45,10 @@ void initialize_slider_values(slider* s)
                s->neighbor_relative_distance_velocity[i] = zeroes_float3();
                s->neighbor_relative_tangential_displacement[i] = zeroes_float3();
                s->neighbor_relative_tangential_velocity[i] = zeroes_float3();
       +        s->neighbor_relative_twist[i] = 0.0;
       +        s->neighbor_relative_twist_velocity[i] = 0.0;
       +        s->neighbor_relative_bend[i] = zeroes_float3();
       +        s->neighbor_relative_bend_velocity[i] = zeroes_float3();
            }
        }
        
       t@@ -365,17 +369,17 @@ void bond_shear_kelvin_voigt(slider* s1, const slider s2,
                (s1->bond_shear_kv_viscosity + s2.bond_shear_kv_viscosity
                 + 1.0e-40);
        
       -    // bond-parallel Kelvin-Voigt elasticity
       +    // bond-normal Kelvin-Voigt elasticity
            const Float3 f_t_elastic =
                multiply_scalar_float3(bond_shear_kv_stiffness,
                        s1->neighbor_relative_tangential_displacement[idx_neighbor]);
        
       -    // bond-parallel Kelvin-Voigt viscosity
       +    // bond-normal Kelvin-Voigt viscosity
            const Float3 f_t_viscous =
                multiply_scalar_float3(bond_shear_kv_viscosity,
                        s1->neighbor_relative_tangential_velocity[idx_neighbor]);
        
       -    // bond-parallel Kelvin-Voigt force, counteracts tension and compression
       +    // bond-normal Kelvin-Voigt force, counteracts shear
            const Float3 f_t = multiply_scalar_float3( -1.0,
                    add_float3(f_t_elastic, f_t_viscous));
        
       t@@ -421,12 +425,81 @@ void bond_shear_kelvin_voigt(slider* s1, const slider s2,
        #endif
        }
        
       +void bond_twist_kelvin_voigt(slider* s1, const slider s2,
       +        const int idx_neighbor)
       +{
       +    // determine the bond-shear KV stiffness from the harmonic mean.
       +    const Float bond_shear_kv_stiffness =
       +        2.*s1->bond_shear_kv_stiffness*s2.bond_shear_kv_stiffness/
       +        (s1->bond_shear_kv_stiffness + s2.bond_shear_kv_stiffness
       +         + 1.0e-40);
       +
       +    // determine the bond-shear KV viscosity from the harmonic mean.
       +    const Float bond_shear_kv_viscosity =
       +        2.*s1->bond_shear_kv_viscosity*s2.bond_shear_kv_viscosity/
       +        (s1->bond_shear_kv_viscosity + s2.bond_shear_kv_viscosity
       +         + 1.0e-40);
       +
       +    // bond-parallel Kelvin-Voigt elasticity
       +    const Float3 m_n_elastic =
       +        multiply_scalar_float3(bond_shear_kv_stiffness,
       +                s1->neighbor_relative_twist[idx_neighbor]);
       +
       +    // bond-parallel Kelvin-Voigt viscosity
       +    const Float3 m_n_viscous =
       +        multiply_scalar_float3(bond_shear_kv_viscosity,
       +                s1->neighbor_relative_twist_velocity[idx_neighbor]);
       +
       +    // bond-parallel Kelvin-Voigt moment, counteracts twisting
       +    const Float3 m_n = multiply_scalar_float3( -1.0,
       +            add_float3(m_n_elastic, m_n_viscous));
       +
       +    // determine torque on slider from shear on this bond
       +    /*const Float3 torque = multiply_scalar_float3( -1.0,
       +                cross_float3(
       +                    multiply_float3_scalar(
       +                        s1->neighbor_distance[idx_neighbor], 0.5),
       +                    f_t));*/
       +
       +    // add bond-shear Kelvin-Voigt force to sum of torques on slider
       +    s1->torque = add_float3(s1->torque, torque);
       +
       +#ifdef DEBUG_SLIDER_FORCE_COMPONENTS
       +    printf("  slider_interaction KV-twist:\n"
       +            "\tf_t    = %f %f %f\n"
       +            "\ttorque = %f %f %f\n"
       +            "\tf_t_elastic = %f %f %f\n"
       +            "\tf_t_viscous = %f %f %f\n"
       +            "\ttan_disp = %f %f %f\n"
       +            "\ttan_vel  = %f %f %f\n",
       +            f_t.x,
       +            f_t.y,
       +            f_t.z,
       +            torque.x,
       +            torque.y,
       +            torque.z,
       +            f_t_elastic.x,
       +            f_t_elastic.y,
       +            f_t_elastic.z,
       +            f_t_viscous.x,
       +            f_t_viscous.y,
       +            f_t_viscous.z,
       +            s1->neighbor_relative_tangential_displacement[idx_neighbor].x,
       +            s1->neighbor_relative_tangential_displacement[idx_neighbor].y,
       +            s1->neighbor_relative_tangential_displacement[idx_neighbor].z,
       +            s1->neighbor_relative_tangential_velocity[idx_neighbor].x,
       +            s1->neighbor_relative_tangential_velocity[idx_neighbor].y,
       +            s1->neighbor_relative_tangential_velocity[idx_neighbor].z);
       +#endif
       +}
       +
        // Resolve bond mechanics between a slider and one of its neighbors based on the 
        // incremental linear-elastic model by Potyondy and Cundall 2004
        void slider_interaction(slider* s1, const slider s2, const int idx_neighbor)
        {
            bond_parallel_kelvin_voigt(s1, s2, idx_neighbor);
            bond_shear_kelvin_voigt(s1, s2, idx_neighbor);
       +    bond_twist_kelvin_voigt(s1, s2, idx_neighbor);
        }
        
        
 (DIR) diff --git a/slidergrid/slider.h b/slidergrid/slider.h
       t@@ -52,9 +52,13 @@ typedef struct {
            // relative spatial movement between this slider and its neighbors
            Float3 neighbor_distance[MAX_NEIGHBORS];
            Float neighbor_relative_distance_displacement[MAX_NEIGHBORS];
       -    Float3 neighbor_relative_distance_velocity[MAX_NEIGHBORS];
       +    Float3 neighbor_relative_distance_velocity[MAX_NEIGHBORS];  // Float?
            Float3 neighbor_relative_tangential_displacement[MAX_NEIGHBORS];
            Float3 neighbor_relative_tangential_velocity[MAX_NEIGHBORS];
       +    Float neighbor_relative_twist[MAX_NEIGHBORS];
       +    Float neighbor_relative_twist_velocity[MAX_NEIGHBORS];
       +    Float3 neighbor_relative_bend[MAX_NEIGHBORS];
       +    Float3 neighbor_relative_bend_velocity[MAX_NEIGHBORS];
        } slider;
        
        
 (DIR) diff --git a/tests/elasticity/Makefile b/tests/elasticity/Makefile
       t@@ -9,7 +9,7 @@ ESSENTIALOBJS=$(SRCFOLDER)/main.o \
                                  $(SRCFOLDER)/grid.o \
                                  $(SRCFOLDER)/vector_math.o \
                                  $(SRCFOLDER)/simulation.o
       -BINS=normal shear
       +BINS=normal shear twist
        
        default: normal-output/normal.E_kin.pdf shear-output/shear.E_kin.pdf
        
       t@@ -23,12 +23,20 @@ shear-output/shear.E_kin.pdf: shear
                python $(ROOT)postprocessing.py --plot-kinetic-energy $<-output
                @#python $(ROOT)postprocessing.py --plot-sliders $<-output
        
       +twist-output/twist.E_kin.pdf: twist
       +        ./$< --verbose
       +        python $(ROOT)postprocessing.py --plot-kinetic-energy $<-output
       +        @#python $(ROOT)postprocessing.py --plot-sliders $<-output
       +
        normal: normal.o $(ESSENTIALOBJS)
                $(CC) $(LDLIBS) $^ -o $@
        
        shear: shear.o $(ESSENTIALOBJS)
                $(CC) $(LDLIBS) $^ -o $@
        
       +twist: twist.o $(ESSENTIALOBJS)
       +        $(CC) $(LDLIBS) $^ -o $@
       +
        clean:
                @$(RM) $(BINS)
                @$(RM) -r *-output