diff --git a/v0.2.1/conf.cpp b/v0.2.1/conf.cpp
index a48ee4e..89f6565 100755
--- a/v0.2.1/conf.cpp
+++ b/v0.2.1/conf.cpp
@@ -371,13 +371,14 @@ void parse_smpl_motions(vector<Motion> &motions, const Json::Value &json, bool&
             trans(fid, i) = tmp;
         }
     }
+    int keep_frames = 10;
     // thetas now have F x 72 in axis-angle representations
     // should parse thetas into one motion
     // motion is Spline<transformation>
     // only one part for SMPL
     motions.resize(1);
     // points is a structure with x and v, where x is the transformation
-    motions[0].points.resize(num_frames + initialization_frames);
+    motions[0].points.resize(num_frames + initialization_frames + keep_frames);
     // consider the interpolation process
     for(int fid = 0;fid < initialization_frames;fid++){
         double scale = double(fid) / double(initialization_frames);
@@ -415,9 +416,32 @@ void parse_smpl_motions(vector<Motion> &motions, const Json::Value &json, bool&
         }
         motions[0].points[fid].t = (fid / fps) - 0;
     }
+    for(int fid = initialization_frames; fid < initialization_frames + keep_frames; fid++){
+        motions[0].points[fid].x.translation = Vec3(trans(0, 0), trans(0, 1), trans(0, 2));
+        for(int pid = 0;pid < 24; pid++){
+            float thetax = thetas(0,3 * pid), thetay = thetas(0,3 * pid + 1), thetaz = thetas(0,3 * pid + 2);
+            // calculate the normalize angle and axis
+            float angle = sqrt(thetax * thetax + thetay * thetay + thetaz * thetaz) + 1e-8;
+            Vec3 axis = Vec3(thetax / angle, thetay / angle, thetaz / angle);
+            Quaternion q = Quaternion::from_axisangle(axis, angle);
+            // rotations contain 24 parts
+            motions[0].points[fid].x.rotations.push_back(q);
+            // setup the translation
+            // nonsense
+            motions[0].points[fid].x.rotation = q;
+            motions[0].points[fid].v.rotations.emplace_back(Quaternion());
 
-    for (int fid = initialization_frames; fid < num_frames + initialization_frames; fid++) {
-        int fid2 = fid - initialization_frames;
+        }
+        motions[0].points[fid].v = Transformation();
+        motions[0].points[fid].x.dynamic_betas.resize(10);
+        for(int i=0;i<10;i++){
+            motions[0].points[fid].x.dynamic_betas[i] = smpl_betas[i];
+        }
+        motions[0].points[fid].t = (fid / fps) - 0;     // here 0 indicates the start time
+    }
+
+    for (int fid = initialization_frames + keep_frames; fid < num_frames + initialization_frames + keep_frames; fid++) {
+        int fid2 = fid - initialization_frames - keep_frames;
         motions[0].points[fid].x.translation = Vec3(trans(fid2, 0), trans(fid2, 1), trans(fid2, 2));
         for(int pid = 0;pid < 24; pid++){
             float thetax = thetas(fid2,3 * pid), thetay = thetas(fid2,3 * pid + 1), thetaz = thetas(fid2,3 * pid + 2);
@@ -442,6 +466,14 @@ void parse_smpl_motions(vector<Motion> &motions, const Json::Value &json, bool&
         motions[0].points[fid].t = (fid / fps) - 0;     // here 0 indicates the start time
 
     }
+    // fill in the velocity
+    for(int fid=0;fid<num_frames+initialization_frames+keep_frames;fid++){
+        if(fid == 0 || fid == num_frames+initialization_frames+keep_frames-1){
+            motions[0].points[fid].v = motions[0].points[fid].x * 0;
+        }else{
+            motions[0].points[fid].v = (motions[0].points[fid + 1].x - motions[0].points[fid-1].x) / (motions[0].points[fid+1].t - motions[0].points[fid-1].t);
+        }
+    }
     Eigen::MatrixXf theta_mat;
     theta_mat.resize(24, 3);
     for(int i=0;i<24;i++){
diff --git a/v0.2.1/displayreplay.cpp b/v0.2.1/displayreplay.cpp
index 79a32b2..d317c2f 100755
--- a/v0.2.1/displayreplay.cpp
+++ b/v0.2.1/displayreplay.cpp
@@ -70,7 +70,12 @@ static void reload () {
     }
     for (int o = 0; o < sim.obstacles.size(); o++) {
         if(!sim.is_smpl)    sim.obstacles[o].get_mesh(sim.time);
-        else    sim.obstacles[o].get_smpl_mesh(sim.time);
+        else{
+//            double decay_time = 0.1, blend = sim.step_time / decay_time;
+//            blend = blend / (1 + blend);
+            sim.obstacles[o].get_smpl_mesh(sim.time);
+//            sim.obstacles[o].smpl_blend_with_previous(sim.time, sim.step_time, blend);
+        }
     }
 }
 
diff --git a/v0.2.1/obstacle.cpp b/v0.2.1/obstacle.cpp
index 5a7c694..94d8523 100755
--- a/v0.2.1/obstacle.cpp
+++ b/v0.2.1/obstacle.cpp
@@ -102,7 +102,7 @@ Mesh& Obstacle::get_smpl_mesh(double time) {
 
         smpl->setAllShapes(dyna_betas);
         smpl->updateModel();
-        // std::cout << "smpl model updated" << std::endl;
+//         std::cout << "smpl model updated" << std::endl;
 //        smpl->saveToOBJ("rest_smpl.obj");
 //        exit(0);
 
@@ -142,3 +142,20 @@ void Obstacle::blend_with_previous (double t, double dt, double blend) {
     }
     compute_ws_data(mesh);
 }
+
+void Obstacle::smpl_blend_with_previous (double t, double dt, double blend) {
+    const Motion *spline = smpl_motion;
+//    Transformation trans = (spline)
+//                           ? get_trans(*spline, t)
+//                             * inverse(get_trans(*spline, t-dt))
+//                           : identity_();
+    Mesh &mesh = curr_state_mesh;
+    for (int n = 0; n < mesh.nodes.size(); n++) {
+        Node *node = mesh.nodes[n];
+//        Vec3 x0 = trans.apply(node->x0);
+        Vec3 x0 = node -> x0;
+        node->x = x0 + blend*(node->x - x0);
+    }
+    compute_ws_data(mesh);
+}
+
diff --git a/v0.2.1/obstacle.hpp b/v0.2.1/obstacle.hpp
index 901c980..9895911 100755
--- a/v0.2.1/obstacle.hpp
+++ b/v0.2.1/obstacle.hpp
@@ -53,6 +53,7 @@ struct Obstacle {
 
     // lerp with previous mesh at time t - dt
     void blend_with_previous (double t, double dt, double blend);
+    void smpl_blend_with_previous (double t, double dt, double blend);
 
     const Motion *transform_spline;
 
diff --git a/v0.2.1/runphysics.cpp b/v0.2.1/runphysics.cpp
index 955b2d4..b245b3c 100755
--- a/v0.2.1/runphysics.cpp
+++ b/v0.2.1/runphysics.cpp
@@ -114,7 +114,7 @@ static void save_timings () {
 void save (const Simulation &sim, int frame) {
     save(sim.cloth_meshes, frame);
     // also to save the obstacles
-    // save(sim.obstacle_meshes, frame, true);
+      save(sim.obstacle_meshes, frame, true);
     save_obstacle_transforms(sim.obstacles, frame, sim.time);
     // std::cout << "end saving obs transform" <<std::endl; 
 }
diff --git a/v0.2.1/simulation.cpp b/v0.2.1/simulation.cpp
index 60aaf7f..8870d57 100755
--- a/v0.2.1/simulation.cpp
+++ b/v0.2.1/simulation.cpp
@@ -141,11 +141,12 @@ void advance_step(Simulation &sim) {
   collision_step(sim);
   // std::cout << "collision step finished" << std::endl; 
   if (sim.step % sim.frame_steps == 0) {
-//     remeshing_step(sim);
+    remeshing_step(sim);
     sim.frame++;
   }
   // only remesh at the beginning
 //  if (sim.step == 1) remeshing_step(sim);
+  // remeshing_step(sim); 
   delete_constraints(cons);
 }
 
@@ -362,6 +363,7 @@ void update_obstacles(Simulation &sim, bool update_positions) {
         sim.obstacles[o].blend_with_previous(sim.time, sim.step_time, blend);
     }else{
         sim.obstacles[o].get_smpl_mesh(sim.time);
+//        sim.obstacles[o].smpl_blend_with_previous(sim.time, sim.step_time, blend);
     }
     if (!update_positions and !sim.is_smpl) {
       // put positions back where they were
@@ -376,7 +378,9 @@ void update_obstacles(Simulation &sim, bool update_positions) {
         for (int n = 0; n < mesh.nodes.size(); n++) {
             Node *node = mesh.nodes[n];
             // set to 0.
-            node->v = (node -> x - node -> x) / sim.step_time;
+            // node->v = (node -> x - node -> x) / sim.step_time;
+            node->v = (node->x - node->x0) / sim.step_time;
+            node->x = node->x0; 
         }
     }
   }
diff --git a/v0.2.1/spline.cpp b/v0.2.1/spline.cpp
index c91bd0e..4ceb675 100755
--- a/v0.2.1/spline.cpp
+++ b/v0.2.1/spline.cpp
@@ -57,15 +57,28 @@ T Spline<T>::pos (double t) const {
          const Point &p0 = points[i-1], &p1 = points[i];
 //         return p0.x;
          double s = (t - p0.t)/(p1.t - p0.t), s2 = s*s, s3 = s2*s;
-
-//         T tmp2 = p0.x*(2*s3 - 3*s2 + 1) + p1.x*(-2*s3 + 3*s2)
-//                  + (p0.v*(s3 - 2*s2 + s) + p1.v*(s3 - s2))*(p1.t - p0.t);
-
-//         return tmp2;
+//         std::cout << "p0, p1t t " << p0.t << " " <<p1.t << " " << t << std::endl;
+//         std::cout << "s, s2, s3 " << s << " " << s2 << " " << s3 << std::endl;
+//
+//         T tmp21 = p0.x*(2*s3 - 3*s2 + 1) + p1.x*(-2*s3 + 3*s2);
+//         std::cout << "tmp21 obtained" << std::endl;
+//         T tmp22 = p0.v*(s3 - 2*s2 + s);
+//         std::cout << "tmp22 obtained" << std::endl;
+//         T tmp23 = p1.v*(s3 - s2);
+//         std::cout << "tmp23 obtained" << std::endl;
+//         T tmp24 = (p0.v*(s3 - 2*s2 + s) + p1.v*(s3 - s2))*(p1.t - p0.t);
+//         T tmp25 = p0.x*(2*s3 - 3*s2 + 1) + p1.x*(-2*s3 + 3*s2);
+//         std::cout << "tmp24 obtained" << std::endl;
+//         T tmp26 = tmp24 + tmp25;
+//        std::cout << "tmp24 obtained" << std::endl;
+         T tmp2 = p0.x*(2*s3 - 3*s2 + 1) + p1.x*(-2*s3 + 3*s2)
+                  + (p0.v*(s3 - 2*s2 + s) + p1.v*(s3 - s2))*(p1.t - p0.t);
+//         std::cout << "tmp2 obtained" << std::endl;
+         return tmp2;
 //         return p0.x*(2*s3 - 3*s2 + 1) + p1.x*(-2*s3 + 3*s2);
 
             // don't consider the velocity term
-            return p0.x * (1 - s) + p1.x * s;
+//            return p0.x * (1 - s) + p1.x * s;
     }
 }
 
diff --git a/v0.2.1/transformation.cpp b/v0.2.1/transformation.cpp
index c3c2c41..e51d823 100755
--- a/v0.2.1/transformation.cpp
+++ b/v0.2.1/transformation.cpp
@@ -280,7 +280,7 @@ Transformation get_smpl_transformation(const Motion &motion, double t){
     for(int i=0;i<T.rotations.size();i++){
         T.rotations[i] = normalize(T.rotations[i]);
     }
-    // std::cout << "T.rotations" << T.rotations.size() << std::endl;
+//     std::cout << "L283, T.rotations obtained" << T.rotations.size() << std::endl;
     return T;
 }