tfind and save projected future slider position - 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 9361f609b8753a5ee25000d7ff5269592b8f82ec
 (DIR) parent ed285da822b8be868f38c8cc1419a2d5fa783592
 (HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
       Date:   Tue, 22 Mar 2016 14:31:48 -0700
       
       find and save projected future slider position
       
       Diffstat:
         M Makefile                            |       2 +-
         M slidergrid/main.c                   |      11 +++++++----
         M slidergrid/slider.c                 |      61 +++++++++++++++++++++++++------
         M slidergrid/slider.h                 |       3 +++
       
       4 files changed, 61 insertions(+), 16 deletions(-)
       ---
 (DIR) diff --git a/Makefile b/Makefile
       t@@ -13,7 +13,7 @@ BIN=test
        default: run-test
        
        run-test: test
       -        ./$<
       +        ./$< --verbose
                @#python postprocessing.py --plot-sliders $<-output
                @#rsync -rav test-output /var/www/html/
        
 (DIR) diff --git a/slidergrid/main.c b/slidergrid/main.c
       t@@ -126,10 +126,13 @@ int main(int argc, char** argv)
                    sim.time <= sim.time_end;
                    sim.time += sim.dt) {
        
       -        // loop over all sliders
       +        for (i=0; i<sim.N; i++)
       +            project_slider_position(&sim.sliders[i], sim.dt, sim.iteration);
       +
       +
       +        // resolve slider-to-slider interaction
                for (i=0; i<sim.N; i++) {
        
       -            // resolve slider-to-slider interaction
                    slider_neighbor_interaction(
                            &sim.sliders[i],
                            sim.sliders,
       t@@ -151,10 +154,10 @@ int main(int argc, char** argv)
                            printf("%d ", sim.sliders[i].neighbors[j]);
                    puts("\n");
        #endif
       +        }
        
       -            // update slider kinematics
       +        for (i=0; i<sim.N; i++) {
                    update_kinematics(&sim.sliders[i], sim.dt, sim.iteration);
       -        }
        
                if (time_since_file >= sim.file_interval) {
                    if (write_simulation_output(&sim, output_folder)) {
 (DIR) diff --git a/slidergrid/slider.c b/slidergrid/slider.c
       t@@ -19,8 +19,17 @@ slider create_slider()
        
        void initialize_slider_values(slider* s)
        {
       +    s->pos_future = zeroes_float3();
       +    s->vel = zeroes_float3();
       +    s->acc = zeroes_float3();
       +    s->angpos_future = zeroes_float3();
       +    s->angpos = zeroes_float3();
       +    s->angvel = zeroes_float3();
       +    s->angacc = zeroes_float3();
       +
            s->force = zeroes_float3();
            s->torque = zeroes_float3();
       +
            s->mass = 0.0;
            s->moment_of_inertia = 0.0;
            s->bond_parallel_stiffness = 0.0;
       t@@ -37,6 +46,23 @@ void initialize_slider_values(slider* s)
        
        }
        
       +
       +// find change in inter-slider distance from old position and projected new 
       +// position, based on current kinematics and explicit temporal integration using 
       +// TY2 expansion
       +void project_slider_position(slider* s, Float dt, long int iteration)
       +{
       +    s->pos_future = make_float3(
       +                s->pos.x + s->vel.x*dt + 0.5*s->acc.x*dt*dt,
       +                s->pos.y + s->vel.y*dt + 0.5*s->acc.y*dt*dt,
       +                s->pos.z + s->vel.z*dt + 0.5*s->acc.z*dt*dt);
       +
       +    s->angpos_future = make_float3(
       +                s->angpos.x + s->angvel.x*dt + 0.5*s->angacc.x*dt*dt,
       +                s->angpos.y + s->angvel.y*dt + 0.5*s->angacc.y*dt*dt,
       +                s->angpos.z + s->angvel.z*dt + 0.5*s->angacc.z*dt*dt);
       +}
       +
        /* Explicit temporal integration scheme based on three-term Taylor expansion.
         * Truncation error O(dt^4) for positions, O(dt^3) for velocities.  Acceleration 
         * change is approximated by backwards differences. */
       t@@ -92,7 +118,7 @@ void update_kinematics(slider* s, Float dt, long int iteration)
        
        // Find the relative displacement and velocity between two sliders
        void slider_displacement(slider* s1, const slider s2,
       -        const int idx_neighbor, const int iteration)
       +        const int idx_neighbor, const int iteration, const Float dt)
        {
        
            // vector pointing from neighbor (s2) position to this slider position (s1)
       t@@ -124,14 +150,17 @@ void slider_displacement(slider* s1, const slider s2,
                //subtract_float3(vel, multiply_float3_scalar(n, dot_float3(n, vel)));
        
            // read previous inter-slider distance vector
       -    Float3 dist0;
       -    if (iteration == 0)
       -        dist0 = zeroes_float3();
       -    else
       -        dist0 = s1->neighbor_distance[idx_neighbor];
       +    const Float3 dist0 = s1->neighbor_distance[idx_neighbor];
       +
       +    const Float3 dist_future = subtract_float3(s1_pos_future, s2_pos_future);
        
       -    // increment in inter-slider distance
       -    const Float3 ddist = subtract_float3(dist, dist0);
       +    // increment in inter-slider distance, divide by two to get displacement 
       +    // over 1 time step
       +    Float3 ddist = divide_float3_scalar(
       +            subtract_float3(dist_future, dist0), 2.0);
       +
       +    //if (iteration == 0)
       +        //ddist = zeroes_float3();
        
            // save current inter-slider distance
            s1->neighbor_distance[idx_neighbor] = dist;
       t@@ -141,7 +170,8 @@ void slider_displacement(slider* s1, const slider s2,
                s1->neighbor_relative_distance_displacement[idx_neighbor] = ddist;
            else
                s1->neighbor_relative_distance_displacement[idx_neighbor] =
       -            add_float3(s1->neighbor_relative_distance_displacement[idx_neighbor],
       +            add_float3(
       +                    s1->neighbor_relative_distance_displacement[idx_neighbor],
                            ddist);
        
            // total relative displacement in inter-slider distance
       t@@ -203,7 +233,8 @@ void slider_neighbor_interaction(
                slider* s,
                const slider* sliders,
                const int N,
       -        const int iteration)
       +        const int iteration,
       +        const Float dt)
        {
            int idx_neighbor;
            for (idx_neighbor=0; idx_neighbor<MAX_NEIGHBORS; idx_neighbor++) {
       t@@ -216,10 +247,18 @@ void slider_neighbor_interaction(
        
                    slider_displacement(
                            s, sliders[s->neighbors[idx_neighbor]],
       -                    idx_neighbor, iteration);
       +                    idx_neighbor, iteration, dt);
       +            printf("- %d: rel_disp = %f %f %f\n",
       +                    idx_neighbor, 
       +                    s->neighbor_relative_distance_displacement[idx_neighbor].x,
       +                    s->neighbor_relative_distance_displacement[idx_neighbor].y,
       +                    s->neighbor_relative_distance_displacement[idx_neighbor].z
       +                  );
        
                    slider_interaction(
                            s, sliders[s->neighbors[idx_neighbor]], idx_neighbor);
       +            printf("s.force = %f %f %f\n",
       +                    s->force.x, s->force.y, s->force.z);
                }
            }
        }
 (DIR) diff --git a/slidergrid/slider.h b/slidergrid/slider.h
       t@@ -8,11 +8,13 @@
        typedef struct {
        
            // linear position, velocity, acceleration of this slider
       +    Float3 pos_future;  // [m]
            Float3 pos;  // [m]
            Float3 vel;  // [m/s]
            Float3 acc;  // [m/(s*s)]
        
            // angular position, velocity, acceleration of this slider
       +    Float3 angpos_future;  // [rad]
            Float3 angpos;  // [rad]
            Float3 angvel;  // [rad/s]
            Float3 angacc;  // [rad/(s*s)]
       t@@ -50,6 +52,7 @@ void print_slider_position(slider s);
        void initialize_slider_values(slider* s);
        //void check_slider_values(slider s);
        
       +void project_slider_position(slider* s, Float dt, long int iteration);
        void update_kinematics(slider* s, Float dt, long int iteration);
        
        void slider_neighbor_interaction(