tresolve contact relative movement - 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 1cbb23226d339f3e3ef3d13fef3ab0141148fcbe
 (DIR) parent 75872d905c2b8343e79dfd89b2122ad4d51b4f63
 (HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
       Date:   Tue, 15 Mar 2016 16:23:06 -0700
       
       resolve contact relative movement
       
       Diffstat:
         M main.c                              |       6 ++++--
         M slider.c                            |      74 ++++++++++++++++++++++++++++---
         M slider.h                            |      43 ++++++++++++++++++++++++++-----
         M vector_math.c                       |      12 ++++++++++++
         M vector_math.h                       |       2 ++
       
       5 files changed, 124 insertions(+), 13 deletions(-)
       ---
 (DIR) diff --git a/main.c b/main.c
       t@@ -13,10 +13,12 @@ int main(int argc, char** argv)
            const int N = nx*ny*nz;
            slider* sliders = create_regular_slider_grid(nx, ny, nz, 1.0, 1.0, 1.0);
        
       -    // set slider masses
       +    // set slider masses and moments of inertia
            int i;
       -    for (i=0; i<N; i++)
       +    for (i=0; i<N; i++) {
                sliders[i].mass = 1.;
       +        sliders[i].moment_of_inertia = 1.;
       +    }
        
            find_and_bond_to_neighbors_n2(sliders, N, 1.5);
        
 (DIR) diff --git a/slider.c b/slider.c
       t@@ -14,30 +14,94 @@ void print_slider_position(slider s)
        void update_kinematics(slider s, Float dt, long int iteration)
        {
            s.acc = divide_float3_scalar(s.force, s.mass);
       +    s.angacc = divide_float3_scalar(s.torque, s.moment_of_inertia);
        
       -    Float3 acc0;
       -    if (iteration == 0)
       +    Float3 acc0, angacc0;
       +    if (iteration == 0) {
                acc0 = make_float3(0., 0., 0.);
       -    else
       +        angacc0 = make_float3(0., 0., 0.);
       +    } else {
                acc0 = s.acc;
       +        angacc0 = s.angacc;
       +    }
        
            const Float3 dacc_dt = make_float3(
                    (s.acc.x - acc0.x)/dt,
                    (s.acc.y - acc0.y)/dt,
                    (s.acc.z - acc0.z)/dt);
        
       +    const Float3 dangacc_dt = make_float3(
       +            (s.angacc.x - angacc0.x)/dt,
       +            (s.angacc.y - angacc0.y)/dt,
       +            (s.angacc.z - angacc0.z)/dt);
       +
            const Float3 dpos_dt = make_float3(
                    s.vel.x*dt + 0.5*s.acc.x*dt*dt + 1./6.*dacc_dt.x*dt*dt*dt,
                    s.vel.y*dt + 0.5*s.acc.y*dt*dt + 1./6.*dacc_dt.y*dt*dt*dt,
                    s.vel.z*dt + 0.5*s.acc.z*dt*dt + 1./6.*dacc_dt.z*dt*dt*dt);
        
       +    const Float3 dangpos_dt = make_float3(
       +            s.angvel.x*dt + 0.5*s.angacc.x*dt*dt + 1./6.*dangacc_dt.x*dt*dt*dt,
       +            s.angvel.y*dt + 0.5*s.angacc.y*dt*dt + 1./6.*dangacc_dt.y*dt*dt*dt,
       +            s.angvel.z*dt + 0.5*s.angacc.z*dt*dt + 1./6.*dangacc_dt.z*dt*dt*dt);
        
            const Float3 dvel_dt = make_float3(
                    s.acc.x*dt + 0.5*dacc_dt.x*dt*dt,
                    s.acc.y*dt + 0.5*dacc_dt.y*dt*dt,
                    s.acc.z*dt + 0.5*dacc_dt.z*dt*dt);
        
       -    s.pos = add_float3(s.pos, dpos_dt);
       -    s.vel = add_float3(s.vel, dvel_dt);
       +    const Float3 dangvel_dt = make_float3(
       +            s.angacc.x*dt + 0.5*dangacc_dt.x*dt*dt,
       +            s.angacc.y*dt + 0.5*dangacc_dt.y*dt*dt,
       +            s.angacc.z*dt + 0.5*dangacc_dt.z*dt*dt);
       +
       +    s.pos    = add_float3(s.pos, dpos_dt);
       +    s.angpos = add_float3(s.angpos, dangpos_dt);
       +    s.vel    = add_float3(s.vel, dvel_dt);
       +    s.angvel = add_float3(s.angvel, dangvel_dt);
       +}
       +
       +// Find the relative displacement and velocity between two sliders
       +void slider_displacement(slider s1, const slider s2, const Float dt)
       +{
       +
       +    const Float3 dist = subtract_float3(s1.pos, s2.pos);
       +    const Float dist_norm = norm_float3(dist);
       +
       +
       +    // relative contact interface velocity w/o rolling
       +    const Float3 vel_linear = subtract_float3(s1.vel, s2.vel);
       +
       +    // relative contact interface velocity with rolling
       +    const Float3 vel = add_float3(vel_linear,
       +            add_float3(
       +                multiply_scalar_float3(dist_norm/2.,
       +                    cross_float3(dist, s1.angvel)),
       +                multiply_scalar_float3(dist_norm/2.,
       +                    cross_float3(dist, s2.angvel))));
       +
       +
        
        }
       +
       +
       +// Resolve bond mechanics between a slider and one of its neighbors
       +
       +
       +// Resolve interaction between a slider and all of its neighbors
       +void slider_neighbor_interaction(
       +        slider s,
       +        const slider* sliders,
       +        const int N)
       +{
       +    int i;
       +    for (i=0; i<MAX_NEIGHBORS; i++) {
       +        if (s.neighbors[i] != -1) {
       +
       +
       +        }
       +    }
       +}
       +
       +
       +
 (DIR) diff --git a/slider.h b/slider.h
       t@@ -6,13 +6,44 @@
        #define MAX_NEIGHBORS 32
        
        typedef struct {
       -    Float3 pos;  // spatial position [m]
       -    Float3 vel;  // spatial velocity [m/s]
       -    Float3 acc;  // spatial acceleration [m/(s*s)]
       -    Float3 force;  // sum of forces [N]
       -    Float mass;  // mass [kg]
       +
       +    // linear position, velocity, acceleration of this slider
       +    Float3 pos;  // [m]
       +    Float3 vel;  // [m/s]
       +    Float3 acc;  // [m/(s*s)]
       +
       +    // angular position, velocity, acceleration of this slider
       +    Float3 angpos;  // [m]
       +    Float3 angvel;  // [m/s]
       +    Float3 angacc;  // [m/(s*s)]
       +
       +    // sum of forces acting on this slider [N]
       +    Float3 force;
       +
       +    // sum of torques acting on this slider [N*m]
       +    Float3 torque;
       +
       +    // slider mass [kg]
       +    Float mass;
       +
       +    // moment of inertia [kg m*m]
       +    Float moment_of_inertia;
       +
       +    // inter-slider contact model parameters
            Float spring_stiffness;  // elastic compressibility [N/m]
       -    int neighbors[MAX_NEIGHBORS]; // indexes of sliders this one is bonded to
       +
       +    // slider-surface contact model parameters
       +    Float coulomb_friction;  // Coulomb frictional value [-]
       +
       +    // The uniquely identifying indexes of the slider neighbors which are bonded 
       +    // to this slider.  A value of -1 denotes an empty field.  Preallocated for 
       +    // speed.
       +    int neighbors[MAX_NEIGHBORS];
       +
       +    // relative spatial movement between this slider and its neighbors
       +    Float3 neighbor_relative_displacement[MAX_NEIGHBORS];
       +    Float3 neighbor_relative_contact_velocity[MAX_NEIGHBORS];
       +
        } slider;
        
        
 (DIR) diff --git a/vector_math.c b/vector_math.c
       t@@ -43,6 +43,13 @@ inline Float3 divide_float3(Float3 v1, Float3 v2)
            return make_float3(v1.x/v2.x, v1.y/v2.y, v1.z/v2.z);
        }
        
       +inline Float3 cross_float3(Float3 v1, Float3 v2)
       +{
       +    return make_float3(
       +            v1.y*v2.z - v1.z*v2.y,
       +            v1.z*v2.x - v1.x*v2.z,
       +            v1.x*v2.y - v1.y*v2.x);
       +}
        
        // vector-scalar operations
        inline Float3 add_float3_scalar(Float3 v, Float s)
       t@@ -70,6 +77,11 @@ inline Float3 multiply_float3_scalar(Float3 v, Float s)
            return make_float3(v.x*s, v.y*s, v.z*s);
        }
        
       +inline Float3 multiply_scalar_float3(Float s, Float3 v)
       +{
       +    return multiply_float3_scalar(v, s);
       +}
       +
        inline Float3 divide_float3_scalar(Float3 v, Float s)
        {
            return make_float3(v.x/s, v.y/s, v.z/s);
 (DIR) diff --git a/vector_math.h b/vector_math.h
       t@@ -12,9 +12,11 @@ Float3 add_float3(Float3 v1, Float3 v2);
        Float3 subtract_float3(Float3 v1, Float3 v2);
        Float3 multiply_float3(Float3 v1, Float3 v2);
        Float3 divide_float3(Float3 v1, Float3 v2);
       +Float3 cross_float3(Float3 v1, Float3 v2);
        
        // vector-scalar operations
        Float3 multiply_float3_scalar(Float3 v, Float s);
       +Float3 multiply_scalar_float3(Float s, Float3 v);
        Float3 divide_float3_scalar(Float3 v, Float s);
        Float3 divide_scalar_float3(Float s, Float3 v);
        Float3 add_float3_scalar(Float3 v, Float s);