Update to bullet master (2.90)

27 files changed, 3659 insertions, 753 deletions

diff --git a/thirdparty/bullet/BulletSoftBody/btConjugateResidual.h b/thirdparty/bullet/BulletSoftBody/btConjugateResidual.h
new file mode 100644
index 0000000000..7b211c4172
--- /dev/null
+++ b/thirdparty/bullet/BulletSoftBody/btConjugateResidual.h
@@ -0,0 +1,188 @@
+/*
+ Written by Xuchen Han <xuchenhan2015@u.northwestern.edu>
+ 
+ Bullet Continuous Collision Detection and Physics Library
+ Copyright (c) 2019 Google Inc. http://bulletphysics.org
+ This software is provided 'as-is', without any express or implied warranty.
+ In no event will the authors be held liable for any damages arising from the use of this software.
+ Permission is granted to anyone to use this software for any purpose,
+ including commercial applications, and to alter it and redistribute it freely,
+ subject to the following restrictions:
+ 1. The origin of this software must not be misrepresented; you must not claim that you wrote the original software. If you use this software in a product, an acknowledgment in the product documentation would be appreciated but is not required.
+ 2. Altered source versions must be plainly marked as such, and must not be misrepresented as being the original software.
+ 3. This notice may not be removed or altered from any source distribution.
+ */
+
+#ifndef BT_CONJUGATE_RESIDUAL_H
+#define BT_CONJUGATE_RESIDUAL_H
+#include <iostream>
+#include <cmath>
+#include <limits>
+#include <LinearMath/btAlignedObjectArray.h>
+#include <LinearMath/btVector3.h>
+#include <LinearMath/btScalar.h>
+#include "LinearMath/btQuickprof.h"
+template <class MatrixX>
+class btConjugateResidual
+{
+    typedef btAlignedObjectArray<btVector3> TVStack;
+    TVStack r,p,z,temp_p, temp_r, best_x;
+    // temp_r = A*r
+    // temp_p = A*p
+    // z = M^(-1) * temp_p = M^(-1) * A * p
+    int max_iterations;
+    btScalar tolerance_squared, best_r;
+public:
+    btConjugateResidual(const int max_it_in)
+    : max_iterations(max_it_in)
+    {
+        tolerance_squared = 1e-2;
+    }
+    
+    virtual ~btConjugateResidual(){}
+    
+    // return the number of iterations taken
+    int solve(MatrixX& A, TVStack& x, const TVStack& b, bool verbose = false)
+    {
+        BT_PROFILE("CRSolve");
+        btAssert(x.size() == b.size());
+        reinitialize(b);
+        // r = b - A * x --with assigned dof zeroed out
+        A.multiply(x, temp_r); // borrow temp_r here to store A*x
+        r = sub(b, temp_r);
+        // z = M^(-1) * r
+        A.precondition(r, z);  // borrow z to store preconditioned r
+        r = z;
+        btScalar residual_norm = norm(r);
+        if (residual_norm <= tolerance_squared) {
+            if (verbose)
+            {
+                std::cout << "Iteration = 0" << std::endl;
+                std::cout << "Two norm of the residual = " << residual_norm << std::endl;
+            }
+            return 0;
+        }
+        p = r;
+        btScalar r_dot_Ar, r_dot_Ar_new;
+        // temp_p = A*p
+        A.multiply(p, temp_p);
+        // temp_r = A*r
+        temp_r = temp_p;
+        r_dot_Ar = dot(r, temp_r);
+        for (int k = 1; k <= max_iterations; k++) {
+            // z = M^(-1) * Ap
+            A.precondition(temp_p, z);
+            // alpha = r^T * A * r / (Ap)^T * M^-1 * Ap)
+            btScalar alpha = r_dot_Ar / dot(temp_p, z);
+            //  x += alpha * p;
+            multAndAddTo(alpha, p, x);
+            //  r -= alpha * z;
+            multAndAddTo(-alpha, z, r);
+            btScalar norm_r = norm(r);
+            if (norm_r < best_r)
+            {
+                best_x = x;
+                best_r = norm_r;
+                if (norm_r < tolerance_squared) {
+                    if (verbose)
+                    {
+                        std::cout << "ConjugateResidual iterations " << k << std::endl;
+                    }
+                    return k;
+                }
+                else
+                {
+                    if (verbose)
+                    {
+                        std::cout << "ConjugateResidual iterations " << k << " has residual "<< norm_r << std::endl;
+                    }
+                }
+            }
+            // temp_r = A * r;
+            A.multiply(r, temp_r);
+            r_dot_Ar_new = dot(r, temp_r);
+            btScalar beta = r_dot_Ar_new/r_dot_Ar;
+            r_dot_Ar = r_dot_Ar_new;
+            // p = beta*p + r;
+            p = multAndAdd(beta, p, r);
+            // temp_p = beta*temp_p + temp_r;
+            temp_p = multAndAdd(beta, temp_p, temp_r);
+        }
+        if (verbose)
+        {
+            std::cout << "ConjugateResidual max iterations reached " << max_iterations << std::endl;
+        }
+        x = best_x;
+        return max_iterations;
+    }
+    
+    void reinitialize(const TVStack& b)
+    {
+        r.resize(b.size());
+        p.resize(b.size());
+        z.resize(b.size());
+        temp_p.resize(b.size());
+        temp_r.resize(b.size());
+        best_x.resize(b.size());
+        best_r = SIMD_INFINITY;
+    }
+    
+    TVStack sub(const TVStack& a, const TVStack& b)
+    {
+        // c = a-b
+        btAssert(a.size() == b.size());
+        TVStack c;
+        c.resize(a.size());
+        for (int i = 0; i < a.size(); ++i)
+        {
+            c[i] = a[i] - b[i];
+        }
+        return c;
+    }
+    
+    btScalar squaredNorm(const TVStack& a)
+    {
+        return dot(a,a);
+    }
+    
+    btScalar norm(const TVStack& a)
+    {
+        btScalar ret = 0;
+        for (int i = 0; i < a.size(); ++i)
+        {
+            for (int d = 0; d < 3; ++d)
+            {
+                ret = btMax(ret, btFabs(a[i][d]));
+            }
+        }
+        return ret;
+    }
+    
+    btScalar dot(const TVStack& a, const TVStack& b)
+    {
+        btScalar ans(0);
+        for (int i = 0; i < a.size(); ++i)
+            ans += a[i].dot(b[i]);
+        return ans;
+    }
+    
+    void multAndAddTo(btScalar s, const TVStack& a, TVStack& result)
+    {
+        //        result += s*a
+        btAssert(a.size() == result.size());
+        for (int i = 0; i < a.size(); ++i)
+            result[i] += s * a[i];
+    }
+    
+    TVStack multAndAdd(btScalar s, const TVStack& a, const TVStack& b)
+    {
+        // result = a*s + b
+        TVStack result;
+        result.resize(a.size());
+        for (int i = 0; i < a.size(); ++i)
+            result[i] = s * a[i] + b[i];
+        return result;
+    }
+};
+#endif /* btConjugateResidual_h */
+
diff --git a/thirdparty/bullet/BulletSoftBody/btDeformableBackwardEulerObjective.cpp b/thirdparty/bullet/BulletSoftBody/btDeformableBackwardEulerObjective.cpp
index 1b247641aa..5381ee6265 100644
--- a/thirdparty/bullet/BulletSoftBody/btDeformableBackwardEulerObjective.cpp
+++ b/thirdparty/bullet/BulletSoftBody/btDeformableBackwardEulerObjective.cpp
@@ -23,12 +23,15 @@ btDeformableBackwardEulerObjective::btDeformableBackwardEulerObjective(btAligned
 , m_backupVelocity(backup_v)
 , m_implicit(false)
 {
-    m_preconditioner = new MassPreconditioner(m_softBodies);
+    m_massPreconditioner = new MassPreconditioner(m_softBodies);
+    m_KKTPreconditioner = new KKTPreconditioner(m_softBodies, m_projection, m_lf, m_dt, m_implicit);
+    m_preconditioner = m_KKTPreconditioner;
 }
 
 btDeformableBackwardEulerObjective::~btDeformableBackwardEulerObjective()
 {
-    delete m_preconditioner;
+    delete m_KKTPreconditioner;
+    delete m_massPreconditioner;
 }
 
 void btDeformableBackwardEulerObjective::reinitialize(bool nodeUpdated, btScalar dt)
@@ -47,7 +50,7 @@ void btDeformableBackwardEulerObjective::reinitialize(bool nodeUpdated, btScalar
         m_lf[i]->reinitialize(nodeUpdated);
     }
     m_projection.reinitialize(nodeUpdated);
-    m_preconditioner->reinitialize(nodeUpdated);
+//    m_preconditioner->reinitialize(nodeUpdated);
 }
 
 void btDeformableBackwardEulerObjective::setDt(btScalar dt)
@@ -80,6 +83,33 @@ void btDeformableBackwardEulerObjective::multiply(const TVStack& x, TVStack& b)
              m_lf[i]->addScaledElasticForceDifferential(-m_dt*m_dt, x, b);
         }
     }
+    int offset = m_nodes.size();
+    for (int i = offset; i < b.size(); ++i)
+    {
+        b[i].setZero();
+    }
+    // add in the lagrange multiplier terms
+    
+    for (int c = 0; c < m_projection.m_lagrangeMultipliers.size(); ++c)
+    {
+        // C^T * lambda
+        const LagrangeMultiplier& lm = m_projection.m_lagrangeMultipliers[c];
+        for (int i = 0; i < lm.m_num_nodes; ++i)
+        {
+            for (int j = 0; j < lm.m_num_constraints; ++j)
+            {
+                b[lm.m_indices[i]] += x[offset+c][j] * lm.m_weights[i] * lm.m_dirs[j];
+            }
+        }
+        // C * x
+        for (int d = 0; d < lm.m_num_constraints; ++d)
+        {
+            for (int i = 0; i < lm.m_num_nodes; ++i)
+            {
+                b[offset+c][d] += lm.m_weights[i] * x[lm.m_indices[i]].dot(lm.m_dirs[d]);
+            }
+        }
+    }
 }
 
 void btDeformableBackwardEulerObjective::updateVelocity(const TVStack& dv)
@@ -134,7 +164,7 @@ void btDeformableBackwardEulerObjective::computeResidual(btScalar dt, TVStack &r
             m_lf[i]->addScaledDampingForce(dt, residual);
         }
     }
-    m_projection.project(residual);
+//    m_projection.project(residual);
 }
 
 btScalar btDeformableBackwardEulerObjective::computeNorm(const TVStack& residual) const
@@ -186,9 +216,9 @@ void btDeformableBackwardEulerObjective::initialGuess(TVStack& dv, const TVStack
 }
 
 //set constraints as projections
-void btDeformableBackwardEulerObjective::setConstraints()
+void btDeformableBackwardEulerObjective::setConstraints(const btContactSolverInfo& infoGlobal)
 {
-    m_projection.setConstraints();
+    m_projection.setConstraints(infoGlobal);
 }
 
 void btDeformableBackwardEulerObjective::applyDynamicFriction(TVStack& r)
diff --git a/thirdparty/bullet/BulletSoftBody/btDeformableBackwardEulerObjective.h b/thirdparty/bullet/BulletSoftBody/btDeformableBackwardEulerObjective.h
index 05ab42ff0a..86579e71ac 100644
--- a/thirdparty/bullet/BulletSoftBody/btDeformableBackwardEulerObjective.h
+++ b/thirdparty/bullet/BulletSoftBody/btDeformableBackwardEulerObjective.h
@@ -15,11 +15,12 @@
 
 #ifndef BT_BACKWARD_EULER_OBJECTIVE_H
 #define BT_BACKWARD_EULER_OBJECTIVE_H
-#include "btConjugateGradient.h"
+//#include "btConjugateGradient.h"
 #include "btDeformableLagrangianForce.h"
 #include "btDeformableMassSpringForce.h"
 #include "btDeformableGravityForce.h"
 #include "btDeformableCorotatedForce.h"
+#include "btDeformableMousePickingForce.h"
 #include "btDeformableLinearElasticityForce.h"
 #include "btDeformableNeoHookeanForce.h"
 #include "btDeformableContactProjection.h"
@@ -39,6 +40,8 @@ public:
     const TVStack& m_backupVelocity;
     btAlignedObjectArray<btSoftBody::Node* > m_nodes;
     bool m_implicit;
+    MassPreconditioner* m_massPreconditioner;
+    KKTPreconditioner* m_KKTPreconditioner;
 
     btDeformableBackwardEulerObjective(btAlignedObjectArray<btSoftBody *>& softBodies, const TVStack& backup_v);
     
@@ -79,7 +82,7 @@ public:
     void updateVelocity(const TVStack& dv);
     
     //set constraints as projections
-    void setConstraints();
+    void setConstraints(const btContactSolverInfo& infoGlobal);
     
     // update the projections and project the residual
     void project(TVStack& r)
@@ -129,6 +132,42 @@ public:
 
     // Calculate the total potential energy in the system
     btScalar totalEnergy(btScalar dt);
+    
+    void addLagrangeMultiplier(const TVStack& vec, TVStack& extended_vec)
+    {
+        extended_vec.resize(vec.size() + m_projection.m_lagrangeMultipliers.size());
+        for (int i = 0; i < vec.size(); ++i)
+        {
+            extended_vec[i] = vec[i];
+        }
+        int offset = vec.size();
+        for (int i = 0; i < m_projection.m_lagrangeMultipliers.size(); ++i)
+        {
+            extended_vec[offset + i].setZero();
+        }
+    }
+    
+    void addLagrangeMultiplierRHS(const TVStack& residual, const TVStack& m_dv, TVStack& extended_residual)
+    {
+        extended_residual.resize(residual.size() + m_projection.m_lagrangeMultipliers.size());
+        for (int i = 0; i < residual.size(); ++i)
+        {
+            extended_residual[i] = residual[i];
+        }
+        int offset = residual.size();
+        for (int i = 0; i < m_projection.m_lagrangeMultipliers.size(); ++i)
+        {
+            const LagrangeMultiplier& lm = m_projection.m_lagrangeMultipliers[i];
+            extended_residual[offset + i].setZero();
+            for (int d = 0; d < lm.m_num_constraints; ++d)
+            {
+                for (int n = 0; n < lm.m_num_nodes; ++n)
+                {
+                    extended_residual[offset + i][d] += lm.m_weights[n] * m_dv[lm.m_indices[n]].dot(lm.m_dirs[d]);
+                }
+            }
+        }
+    }
 };
 
 #endif /* btBackwardEulerObjective_h */
diff --git a/thirdparty/bullet/BulletSoftBody/btDeformableBodySolver.cpp b/thirdparty/bullet/BulletSoftBody/btDeformableBodySolver.cpp
index 7724a8ec69..132699c54f 100644
--- a/thirdparty/bullet/BulletSoftBody/btDeformableBodySolver.cpp
+++ b/thirdparty/bullet/BulletSoftBody/btDeformableBodySolver.cpp
@@ -18,13 +18,15 @@
 #include "btDeformableBodySolver.h"
 #include "btSoftBodyInternals.h"
 #include "LinearMath/btQuickprof.h"
-static const int kMaxConjugateGradientIterations  = 50;
+static const int kMaxConjugateGradientIterations = 50;
 btDeformableBodySolver::btDeformableBodySolver()
 : m_numNodes(0)
 , m_cg(kMaxConjugateGradientIterations)
+, m_cr(kMaxConjugateGradientIterations)
 , m_maxNewtonIterations(5)
 , m_newtonTolerance(1e-4)
 , m_lineSearch(false)
+, m_useProjection(false)
 {
     m_objective = new btDeformableBackwardEulerObjective(m_softBodies, m_backupVelocity);
 }
@@ -41,7 +43,22 @@ void btDeformableBodySolver::solveDeformableConstraints(btScalar solverdt)
     {
         m_objective->computeResidual(solverdt, m_residual);
         m_objective->applyDynamicFriction(m_residual);
-        computeStep(m_dv, m_residual);
+        if (m_useProjection)
+        {
+            computeStep(m_dv, m_residual);
+        }
+        else
+        {
+            TVStack rhs, x;
+            m_objective->addLagrangeMultiplierRHS(m_residual, m_dv, rhs);
+            m_objective->addLagrangeMultiplier(m_dv, x);
+            m_objective->m_preconditioner->reinitialize(true);
+            computeStep(x, rhs);
+            for (int i = 0; i<m_dv.size(); ++i)
+            {
+                    m_dv[i] = x[i];
+            }
+        }
         updateVelocity();
     }
     else
@@ -63,7 +80,7 @@ void btDeformableBodySolver::solveDeformableConstraints(btScalar solverdt)
                     ++counter;
                 }
             }
-            
+
             m_objective->computeResidual(solverdt, m_residual);
             if (m_objective->computeNorm(m_residual) < m_newtonTolerance && i > 0)
             {
@@ -200,7 +217,10 @@ void btDeformableBodySolver::updateDv(btScalar scale)
 
 void btDeformableBodySolver::computeStep(TVStack& ddv, const TVStack& residual)
 {
-    m_cg.solve(*m_objective, ddv, residual);
+    if (m_useProjection)
+        m_cg.solve(*m_objective, ddv, residual, false);
+    else
+        m_cr.solve(*m_objective, ddv, residual, false);
 }
 
 void btDeformableBodySolver::reinitialize(const btAlignedObjectArray<btSoftBody *>& softBodies, btScalar dt)
@@ -226,27 +246,22 @@ void btDeformableBodySolver::reinitialize(const btAlignedObjectArray<btSoftBody
     
     m_dt = dt;
     m_objective->reinitialize(nodeUpdated, dt);
+    updateSoftBodies();
 }
 
-void btDeformableBodySolver::setConstraints()
+void btDeformableBodySolver::setConstraints(const btContactSolverInfo& infoGlobal)
 {
     BT_PROFILE("setConstraint");
-    m_objective->setConstraints();
+    m_objective->setConstraints(infoGlobal);
 }
 
-btScalar btDeformableBodySolver::solveContactConstraints(btCollisionObject** deformableBodies,int numDeformableBodies)
+btScalar btDeformableBodySolver::solveContactConstraints(btCollisionObject** deformableBodies,int numDeformableBodies, const btContactSolverInfo& infoGlobal)
 {
     BT_PROFILE("solveContactConstraints");
-    btScalar maxSquaredResidual = m_objective->m_projection.update(deformableBodies,numDeformableBodies);
+    btScalar maxSquaredResidual = m_objective->m_projection.update(deformableBodies,numDeformableBodies, infoGlobal);
     return maxSquaredResidual;
 }
 
-btScalar btDeformableBodySolver::solveSplitImpulse(const btContactSolverInfo& infoGlobal)
-{
-    BT_PROFILE("solveSplitImpulse");
-    return m_objective->m_projection.solveSplitImpulse(infoGlobal);
-}
-
 void btDeformableBodySolver::splitImpulseSetup(const btContactSolverInfo& infoGlobal)
 {
      m_objective->m_projection.splitImpulseSetup(infoGlobal);
@@ -333,8 +348,10 @@ void btDeformableBodySolver::setupDeformableSolve(bool implicit)
                 m_backupVelocity[counter] = psb->m_nodes[j].m_vn;
             }
             else
+            {
                 m_dv[counter] =  psb->m_nodes[j].m_v - m_backupVelocity[counter];
-            psb->m_nodes[j].m_v = m_backupVelocity[counter] + psb->m_nodes[j].m_vsplit;
+            }
+            psb->m_nodes[j].m_v = m_backupVelocity[counter];
             ++counter;
         }
     }
@@ -385,6 +402,7 @@ void btDeformableBodySolver::predictMotion(btScalar solverdt)
 
 void btDeformableBodySolver::predictDeformableMotion(btSoftBody* psb, btScalar dt)
 {
+    BT_PROFILE("btDeformableBodySolver::predictDeformableMotion");
     int i, ni;
     
     /* Update                */
@@ -423,40 +441,22 @@ void btDeformableBodySolver::predictDeformableMotion(btSoftBody* psb, btScalar d
             n.m_v *= max_v;
         }
         n.m_q = n.m_x + n.m_v * dt;
+        n.m_penetration = 0;
     }
 
     /* Nodes                */
-    ATTRIBUTE_ALIGNED16(btDbvtVolume)
-    vol;
-    for (i = 0, ni = psb->m_nodes.size(); i < ni; ++i)
-    {
-        btSoftBody::Node& n = psb->m_nodes[i];
-        btVector3 points[2] = {n.m_x, n.m_q};
-        vol = btDbvtVolume::FromPoints(points, 2);
-        vol.Expand(btVector3(psb->m_sst.radmrg, psb->m_sst.radmrg, psb->m_sst.radmrg));
-        psb->m_ndbvt.update(n.m_leaf, vol);
-    }
-
+    psb->updateNodeTree(true, true);
     if (!psb->m_fdbvt.empty())
     {
-        for (int i = 0; i < psb->m_faces.size(); ++i)
-        {
-            btSoftBody::Face& f = psb->m_faces[i];
-            btVector3 points[6] = {f.m_n[0]->m_x, f.m_n[0]->m_q,
-                                   f.m_n[1]->m_x, f.m_n[1]->m_q,
-                                   f.m_n[2]->m_x, f.m_n[2]->m_q};
-            vol = btDbvtVolume::FromPoints(points, 6);
-            vol.Expand(btVector3(psb->m_sst.radmrg, psb->m_sst.radmrg, psb->m_sst.radmrg));
-            psb->m_fdbvt.update(f.m_leaf, vol);
-        }
+        psb->updateFaceTree(true, true);
     }
-    /* Clear contacts        */
+    /* Clear contacts */
     psb->m_nodeRigidContacts.resize(0);
     psb->m_faceRigidContacts.resize(0);
     psb->m_faceNodeContacts.resize(0);
     /* Optimize dbvt's        */
-    psb->m_ndbvt.optimizeIncremental(1);
-    psb->m_fdbvt.optimizeIncremental(1);
+//    psb->m_ndbvt.optimizeIncremental(1);
+//    psb->m_fdbvt.optimizeIncremental(1);
 }
 
 
diff --git a/thirdparty/bullet/BulletSoftBody/btDeformableBodySolver.h b/thirdparty/bullet/BulletSoftBody/btDeformableBodySolver.h
index f78a8f696b..d4e5f4c603 100644
--- a/thirdparty/bullet/BulletSoftBody/btDeformableBodySolver.h
+++ b/thirdparty/bullet/BulletSoftBody/btDeformableBodySolver.h
@@ -22,7 +22,8 @@
 #include "btDeformableMultiBodyDynamicsWorld.h"
 #include "BulletDynamics/Featherstone/btMultiBodyLinkCollider.h"
 #include "BulletDynamics/Featherstone/btMultiBodyConstraint.h"
-
+#include "btConjugateResidual.h"
+#include "btConjugateGradient.h"
 struct btCollisionObjectWrapper;
 class btDeformableBackwardEulerObjective;
 class btDeformableMultiBodyDynamicsWorld;
@@ -40,14 +41,15 @@ protected:
     TVStack m_backupVelocity;       // backed up v, equals v_n for implicit, equals v_{n+1}^* for explicit
     btScalar m_dt;                  // dt
     btConjugateGradient<btDeformableBackwardEulerObjective> m_cg;  // CG solver
+    btConjugateResidual<btDeformableBackwardEulerObjective> m_cr;  // CR solver
     bool m_implicit;                // use implicit scheme if true, explicit scheme if false
     int m_maxNewtonIterations;      // max number of newton iterations
     btScalar m_newtonTolerance;     // stop newton iterations if f(x) < m_newtonTolerance
     bool m_lineSearch;              // If true, use newton's method with line search under implicit scheme
-    
 public:
     // handles data related to objective function
     btDeformableBackwardEulerObjective* m_objective;
+    bool m_useProjection;
     
     btDeformableBodySolver();
     
@@ -61,15 +63,11 @@ public:
     // update soft body normals
     virtual void updateSoftBodies();
     
+    virtual btScalar solveContactConstraints(btCollisionObject** deformableBodies,int numDeformableBodies, const btContactSolverInfo& infoGlobal);
+    
     // solve the momentum equation
     virtual void solveDeformableConstraints(btScalar solverdt);
     
-    // solve the contact between deformable and rigid as well as among deformables
-    btScalar solveContactConstraints(btCollisionObject** deformableBodies,int numDeformableBodies);
-    
-    // solve the position error  between deformable and rigid as well as among deformables;
-    btScalar solveSplitImpulse(const btContactSolverInfo& infoGlobal);
-    
     // set up the position error in split impulse
     void splitImpulseSetup(const btContactSolverInfo& infoGlobal);
 
@@ -77,7 +75,7 @@ public:
     void reinitialize(const btAlignedObjectArray<btSoftBody *>& softBodies, btScalar dt);
     
     // set up contact constraints
-    void setConstraints();
+    void setConstraints(const btContactSolverInfo& infoGlobal);
     
     // add in elastic forces and gravity to obtain v_{n+1}^* and calls predictDeformableMotion
     virtual void predictMotion(btScalar solverdt);
diff --git a/thirdparty/bullet/BulletSoftBody/btDeformableContactConstraint.cpp b/thirdparty/bullet/BulletSoftBody/btDeformableContactConstraint.cpp
index e8219dc50e..2864446de6 100644
--- a/thirdparty/bullet/BulletSoftBody/btDeformableContactConstraint.cpp
+++ b/thirdparty/bullet/BulletSoftBody/btDeformableContactConstraint.cpp
@@ -15,9 +15,9 @@
 
 #include "btDeformableContactConstraint.h"
 /* ================   Deformable Node Anchor   =================== */
-btDeformableNodeAnchorConstraint::btDeformableNodeAnchorConstraint(const btSoftBody::DeformableNodeRigidAnchor& a)
+btDeformableNodeAnchorConstraint::btDeformableNodeAnchorConstraint(const btSoftBody::DeformableNodeRigidAnchor& a, const btContactSolverInfo& infoGlobal)
 : m_anchor(&a)
-, btDeformableContactConstraint(a.m_cti.m_normal)
+, btDeformableContactConstraint(a.m_cti.m_normal, infoGlobal)
 {
 }
 
@@ -79,14 +79,14 @@ btVector3 btDeformableNodeAnchorConstraint::getVa() const
     return va;
 }
 
-btScalar btDeformableNodeAnchorConstraint::solveConstraint()
+btScalar btDeformableNodeAnchorConstraint::solveConstraint(const btContactSolverInfo& infoGlobal)
 {
     const btSoftBody::sCti& cti = m_anchor->m_cti;
     btVector3 va = getVa();
     btVector3 vb = getVb();
     btVector3 vr = (vb - va);
     // + (m_anchor->m_node->m_x - cti.m_colObj->getWorldTransform() * m_anchor->m_local) * 10.0
-    const btScalar dn = btDot(vr, cti.m_normal);
+    const btScalar dn = btDot(vr, vr);
     // dn is the normal component of velocity diffrerence. Approximates the residual. // todo xuchenhan@: this prob needs to be scaled by dt
     btScalar residualSquare = dn*dn;
     btVector3 impulse = m_anchor->m_c0 * vr;
@@ -134,14 +134,15 @@ void btDeformableNodeAnchorConstraint::applyImpulse(const btVector3& impulse)
 }
 
 /* ================   Deformable vs. Rigid   =================== */
-btDeformableRigidContactConstraint::btDeformableRigidContactConstraint(const btSoftBody::DeformableRigidContact& c)
+btDeformableRigidContactConstraint::btDeformableRigidContactConstraint(const btSoftBody::DeformableRigidContact& c, const btContactSolverInfo& infoGlobal)
 : m_contact(&c)
-, btDeformableContactConstraint(c.m_cti.m_normal)
+, btDeformableContactConstraint(c.m_cti.m_normal, infoGlobal)
 {
     m_total_normal_dv.setZero();
     m_total_tangent_dv.setZero();
-    // penetration is non-positive. The magnitude of penetration is the depth of penetration.
-    m_penetration = btMin(btScalar(0), c.m_cti.m_offset);
+    // The magnitude of penetration is the depth of penetration.
+    m_penetration = c.m_cti.m_offset;
+//	m_penetration = btMin(btScalar(0),c.m_cti.m_offset);
 }
 
 btDeformableRigidContactConstraint::btDeformableRigidContactConstraint(const btDeformableRigidContactConstraint& other)
@@ -206,16 +207,16 @@ btVector3 btDeformableRigidContactConstraint::getVa() const
     return va;
 }
 
-btScalar btDeformableRigidContactConstraint::solveConstraint()
+btScalar btDeformableRigidContactConstraint::solveConstraint(const btContactSolverInfo& infoGlobal)
 {
     const btSoftBody::sCti& cti = m_contact->m_cti;
     btVector3 va = getVa();
     btVector3 vb = getVb();
     btVector3 vr = vb - va;
-    const btScalar dn = btDot(vr, cti.m_normal);
+    btScalar dn = btDot(vr, cti.m_normal) + m_penetration * infoGlobal.m_deformable_erp / infoGlobal.m_timeStep;
     // dn is the normal component of velocity diffrerence. Approximates the residual. // todo xuchenhan@: this prob needs to be scaled by dt
     btScalar residualSquare = dn*dn;
-    btVector3 impulse = m_contact->m_c0 * vr;
+    btVector3 impulse = m_contact->m_c0 * (vr + m_penetration * infoGlobal.m_deformable_erp / infoGlobal.m_timeStep * cti.m_normal) ;
     const btVector3 impulse_normal = m_contact->m_c0 * (cti.m_normal * dn);
     btVector3 impulse_tangent = impulse - impulse_normal;
     btVector3 old_total_tangent_dv = m_total_tangent_dv;
@@ -256,6 +257,8 @@ btScalar btDeformableRigidContactConstraint::solveConstraint()
     impulse = impulse_normal + impulse_tangent;
     // apply impulse to deformable nodes involved and change their velocities
     applyImpulse(impulse);
+	if (residualSquare < 1e-7)
+		return residualSquare;
     // apply impulse to the rigid/multibodies involved and change their velocities
     if (cti.m_colObj->getInternalType() == btCollisionObject::CO_RIGID_BODY)
     {
@@ -285,43 +288,17 @@ btScalar btDeformableRigidContactConstraint::solveConstraint()
             }
         }
     }
+//    va = getVa();
+//    vb = getVb();
+//    vr = vb - va;
+//    btScalar dn1 = btDot(vr, cti.m_normal) / 150;
+//    m_penetration += dn1;
     return residualSquare;
 }
-
-btScalar btDeformableRigidContactConstraint::solveSplitImpulse(const btContactSolverInfo& infoGlobal)
-{
-    const btSoftBody::sCti& cti = m_contact->m_cti;
-    const btScalar dn = m_penetration;
-    if (dn != 0)
-    {
-        const btVector3 impulse = (m_contact->m_c0 * (cti.m_normal * dn / infoGlobal.m_timeStep));
-        // one iteration of the position impulse corrects all the position error at this timestep
-        m_penetration -= dn;
-        // apply impulse to deformable nodes involved and change their position
-        applySplitImpulse(impulse);
-        // apply impulse to the rigid/multibodies involved and change their position
-        if (cti.m_colObj->getInternalType() == btCollisionObject::CO_RIGID_BODY)
-        {
-            btRigidBody* rigidCol = 0;
-            rigidCol = (btRigidBody*)btRigidBody::upcast(cti.m_colObj);
-            if (rigidCol)
-            {
-                rigidCol->applyPushImpulse(impulse, m_contact->m_c1);
-            }
-        }
-        else if (cti.m_colObj->getInternalType() == btCollisionObject::CO_FEATHERSTONE_LINK)
-        {
-            // todo xuchenhan@
-        }
-        return (m_penetration/infoGlobal.m_timeStep) * (m_penetration/infoGlobal.m_timeStep);
-    }
-    return 0;
-}
-
 /* ================   Node vs. Rigid   =================== */
-btDeformableNodeRigidContactConstraint::btDeformableNodeRigidContactConstraint(const btSoftBody::DeformableNodeRigidContact& contact)
+btDeformableNodeRigidContactConstraint::btDeformableNodeRigidContactConstraint(const btSoftBody::DeformableNodeRigidContact& contact, const btContactSolverInfo& infoGlobal)
     : m_node(contact.m_node)
-    , btDeformableRigidContactConstraint(contact)
+    , btDeformableRigidContactConstraint(contact, infoGlobal)
     {
     }
 
@@ -349,22 +326,17 @@ void btDeformableNodeRigidContactConstraint::applyImpulse(const btVector3& impul
     contact->m_node->m_v -= dv;
 }
 
-void btDeformableNodeRigidContactConstraint::applySplitImpulse(const btVector3& impulse)
-{
-    const btSoftBody::DeformableNodeRigidContact* contact = getContact();
-    btVector3 dv = impulse * contact->m_c2;
-    contact->m_node->m_vsplit -= dv;
-};
-
 /* ================   Face vs. Rigid   =================== */
-btDeformableFaceRigidContactConstraint::btDeformableFaceRigidContactConstraint(const btSoftBody::DeformableFaceRigidContact& contact)
+btDeformableFaceRigidContactConstraint::btDeformableFaceRigidContactConstraint(const btSoftBody::DeformableFaceRigidContact& contact, const btContactSolverInfo& infoGlobal, bool useStrainLimiting)
 : m_face(contact.m_face)
-, btDeformableRigidContactConstraint(contact)
+, m_useStrainLimiting(useStrainLimiting)
+, btDeformableRigidContactConstraint(contact, infoGlobal)
 {
 }
 
 btDeformableFaceRigidContactConstraint::btDeformableFaceRigidContactConstraint(const btDeformableFaceRigidContactConstraint& other)
 : m_face(other.m_face)
+, m_useStrainLimiting(other.m_useStrainLimiting)
 , btDeformableRigidContactConstraint(other)
 {
 }
@@ -411,47 +383,70 @@ void btDeformableFaceRigidContactConstraint::applyImpulse(const btVector3& impul
         v1 -= dv * contact->m_weights[1];
     if (im2 > 0)
         v2 -= dv * contact->m_weights[2];
-    
-    // apply strain limiting to prevent undamped modes
-    btScalar m01 = (btScalar(1)/(im0 + im1));
-    btScalar m02 = (btScalar(1)/(im0 + im2));
-    btScalar m12 = (btScalar(1)/(im1 + im2));
-    
-    btVector3 dv0 = im0 * (m01 * (v1-v0) + m02 * (v2-v0));
-    btVector3 dv1 = im1 * (m01 * (v0-v1) + m12 * (v2-v1));
-    btVector3 dv2 = im2 * (m12 * (v1-v2) + m02 * (v0-v2));
-    
-    v0 += dv0;
-    v1 += dv1;
-    v2 += dv2;
-}
-
-void btDeformableFaceRigidContactConstraint::applySplitImpulse(const btVector3& impulse)
-{
-    const btSoftBody::DeformableFaceRigidContact* contact = getContact();
-    btVector3 dv = impulse * contact->m_c2;
-    btSoftBody::Face* face = contact->m_face;
-
-    btVector3& v0 = face->m_n[0]->m_vsplit;
-    btVector3& v1 = face->m_n[1]->m_vsplit;
-    btVector3& v2 = face->m_n[2]->m_vsplit;
-    const btScalar& im0 = face->m_n[0]->m_im;
-    const btScalar& im1 = face->m_n[1]->m_im;
-    const btScalar& im2 = face->m_n[2]->m_im;
-    if (im0 > 0)
-        v0 -= dv * contact->m_weights[0];
-    if (im1 > 0)
-        v1 -= dv * contact->m_weights[1];
-    if (im2 > 0)
-        v2 -= dv * contact->m_weights[2];
+	if (m_useStrainLimiting)
+	{
+		btScalar relaxation = 1./btScalar(m_infoGlobal->m_numIterations);
+		btScalar m01 = (relaxation/(im0 + im1));
+		btScalar m02 = (relaxation/(im0 + im2));
+		btScalar m12 = (relaxation/(im1 + im2));
+		#ifdef USE_STRAIN_RATE_LIMITING
+		// apply strain limiting to prevent the new velocity to change the current length of the edge by more than 1%.
+		btScalar p = 0.01;
+		btVector3& x0 = face->m_n[0]->m_x;
+		btVector3& x1 = face->m_n[1]->m_x;
+		btVector3& x2 = face->m_n[2]->m_x;
+		const btVector3 x_diff[3] = {x1-x0, x2-x0, x2-x1};
+		const btVector3 v_diff[3] = {v1-v0, v2-v0, v2-v1};
+		btVector3 u[3];
+		btScalar x_diff_dot_u, dn[3];
+		btScalar dt = m_infoGlobal->m_timeStep;
+		for (int i = 0; i < 3; ++i)
+		{
+			btScalar x_diff_norm = x_diff[i].safeNorm();
+			btScalar x_diff_norm_new = (x_diff[i] + v_diff[i] * dt).safeNorm();
+			btScalar strainRate = x_diff_norm_new/x_diff_norm;
+			u[i] = v_diff[i];
+			u[i].safeNormalize();
+			if (x_diff_norm == 0 || (1-p <= strainRate && strainRate <= 1+p))
+			{
+				dn[i] = 0;
+				continue;
+			}
+			x_diff_dot_u = btDot(x_diff[i], u[i]);
+			btScalar s;
+			if (1-p > strainRate)
+			{
+				s = 1/dt * (-x_diff_dot_u - btSqrt(x_diff_dot_u*x_diff_dot_u + (p*p-2*p) * x_diff_norm * x_diff_norm));
+			}
+			else
+			{
+				s = 1/dt * (-x_diff_dot_u + btSqrt(x_diff_dot_u*x_diff_dot_u + (p*p+2*p) * x_diff_norm * x_diff_norm));
+			}
+			//		x_diff_norm_new = (x_diff[i] + s * u[i] * dt).safeNorm();
+			//		strainRate = x_diff_norm_new/x_diff_norm;
+			dn[i] = s - v_diff[i].safeNorm();
+		}
+		btVector3 dv0 = im0 * (m01 * u[0]*(-dn[0]) + m02 * u[1]*-(dn[1]));
+		btVector3 dv1 = im1 * (m01 * u[0]*(dn[0]) + m12 * u[2]*(-dn[2]));
+		btVector3 dv2 = im2 * (m12 * u[2]*(dn[2]) + m02 * u[1]*(dn[1]));
+	#else
+		// apply strain limiting to prevent undamped modes
+		btVector3 dv0 = im0 * (m01 * (v1-v0) + m02 * (v2-v0));
+		btVector3 dv1 = im1 * (m01 * (v0-v1) + m12 * (v2-v1));
+		btVector3 dv2 = im2 * (m12 * (v1-v2) + m02 * (v0-v2));
+	#endif
+		v0 += dv0;
+		v1 += dv1;
+		v2 += dv2;
+	}
 }
 
 /* ================   Face vs. Node   =================== */
-btDeformableFaceNodeContactConstraint::btDeformableFaceNodeContactConstraint(const btSoftBody::DeformableFaceNodeContact& contact)
+btDeformableFaceNodeContactConstraint::btDeformableFaceNodeContactConstraint(const btSoftBody::DeformableFaceNodeContact& contact, const btContactSolverInfo& infoGlobal)
 : m_node(contact.m_node)
 , m_face(contact.m_face)
 , m_contact(&contact)
-, btDeformableContactConstraint(contact.m_normal)
+, btDeformableContactConstraint(contact.m_normal, infoGlobal)
 {
     m_total_normal_dv.setZero();
     m_total_tangent_dv.setZero();
@@ -487,7 +482,7 @@ btVector3 btDeformableFaceNodeContactConstraint::getDv(const btSoftBody::Node* n
     return dv * contact->m_weights[2];
 }
 
-btScalar btDeformableFaceNodeContactConstraint::solveConstraint()
+btScalar btDeformableFaceNodeContactConstraint::solveConstraint(const btContactSolverInfo& infoGlobal)
 {
     btVector3 va = getVa();
     btVector3 vb = getVb();
@@ -577,15 +572,4 @@ void btDeformableFaceNodeContactConstraint::applyImpulse(const btVector3& impuls
     {
         v2 -= dvb * contact->m_weights[2];
     }
-    // todo: Face node constraints needs more work
-//    btScalar m01 = (btScalar(1)/(im0 + im1));
-//    btScalar m02 = (btScalar(1)/(im0 + im2));
-//    btScalar m12 = (btScalar(1)/(im1 + im2));
-//
-//    btVector3 dv0 = im0 * (m01 * (v1-v0) + m02 * (v2-v0));
-//    btVector3 dv1 = im1 * (m01 * (v0-v1) + m12 * (v2-v1));
-//    btVector3 dv2 = im2 * (m12 * (v1-v2) + m02 * (v0-v2));
-//    v0 += dv0;
-//    v1 += dv1;
-//    v2 += dv2;
 }
diff --git a/thirdparty/bullet/BulletSoftBody/btDeformableContactConstraint.h b/thirdparty/bullet/BulletSoftBody/btDeformableContactConstraint.h
index 912119e7c3..9f9d5bf0a3 100644
--- a/thirdparty/bullet/BulletSoftBody/btDeformableContactConstraint.h
+++ b/thirdparty/bullet/BulletSoftBody/btDeformableContactConstraint.h
@@ -24,34 +24,33 @@ public:
     // True if the friction is static
     // False if the friction is dynamic
     bool m_static;
-    
-    // normal of the contact
-    btVector3 m_normal;
-    
-    btDeformableContactConstraint(const btVector3& normal): m_static(false), m_normal(normal)
-    {
-    }
-    
-    btDeformableContactConstraint(bool isStatic, const btVector3& normal): m_static(isStatic), m_normal(normal)
-    {
-    }
-    
-    btDeformableContactConstraint(const btDeformableContactConstraint& other)
-    : m_static(other.m_static)
-    , m_normal(other.m_normal)
-    {
-        
-    }
-    btDeformableContactConstraint(){}
-    
+	const btContactSolverInfo* m_infoGlobal;
+
+	// normal of the contact
+	btVector3 m_normal;
+
+	btDeformableContactConstraint(const btVector3& normal, const btContactSolverInfo& infoGlobal): m_static(false), m_normal(normal), m_infoGlobal(&infoGlobal)
+	{
+	}
+
+	btDeformableContactConstraint(bool isStatic, const btVector3& normal, const btContactSolverInfo& infoGlobal): m_static(isStatic), m_normal(normal), m_infoGlobal(&infoGlobal)
+	{
+	}
+	
+	btDeformableContactConstraint(){}
+
+	btDeformableContactConstraint(const btDeformableContactConstraint& other)
+	: m_static(other.m_static)
+	, m_normal(other.m_normal)
+	, m_infoGlobal(other.m_infoGlobal)
+	{
+	}
+
     virtual ~btDeformableContactConstraint(){}
     
     // solve the constraint with inelastic impulse and return the error, which is the square of normal component of velocity diffrerence
     // the constraint is solved by calculating the impulse between object A and B in the contact and apply the impulse to both objects involved in the contact
-    virtual btScalar solveConstraint() = 0;
-    
-    // solve the position error by applying an inelastic impulse that changes only the position (not velocity)
-    virtual btScalar solveSplitImpulse(const btContactSolverInfo& infoGlobal) = 0;
+    virtual btScalar solveConstraint(const btContactSolverInfo& infoGlobal) = 0;
     
     // get the velocity of the object A in the contact
     virtual btVector3 getVa() const = 0;
@@ -65,9 +64,6 @@ public:
     // apply impulse to the soft body node and/or face involved
     virtual void applyImpulse(const btVector3& impulse) = 0;
     
-    // apply position based impulse to the soft body node and/or face involved
-    virtual void applySplitImpulse(const btVector3& impulse) = 0;
-    
     // scale the penetration depth by erp
     virtual void setPenetrationScale(btScalar scale) = 0;
 };
@@ -77,29 +73,21 @@ public:
 class btDeformableStaticConstraint : public btDeformableContactConstraint
 {
 public:
-    const btSoftBody::Node* m_node;
-    
-    btDeformableStaticConstraint(){}
+    btSoftBody::Node* m_node;
     
-    btDeformableStaticConstraint(const btSoftBody::Node* node): m_node(node), btDeformableContactConstraint(false, btVector3(0,0,0))
+    btDeformableStaticConstraint(btSoftBody::Node* node, const btContactSolverInfo& infoGlobal): m_node(node), btDeformableContactConstraint(false, btVector3(0,0,0), infoGlobal)
     {
     }
-    
+	btDeformableStaticConstraint(){}
     btDeformableStaticConstraint(const btDeformableStaticConstraint& other)
     : m_node(other.m_node)
     , btDeformableContactConstraint(other)
     {
-        
     }
     
     virtual ~btDeformableStaticConstraint(){}
     
-    virtual btScalar solveConstraint()
-    {
-        return 0;
-    }
-    
-    virtual btScalar solveSplitImpulse(const btContactSolverInfo& infoGlobal)
+    virtual btScalar solveConstraint(const btContactSolverInfo& infoGlobal)
     {
         return 0;
     }
@@ -120,7 +108,6 @@ public:
     }
     
     virtual void applyImpulse(const btVector3& impulse){}
-    virtual void applySplitImpulse(const btVector3& impulse){}
     virtual void setPenetrationScale(btScalar scale){}
 };
 
@@ -130,19 +117,15 @@ class btDeformableNodeAnchorConstraint : public btDeformableContactConstraint
 {
 public:
     const btSoftBody::DeformableNodeRigidAnchor* m_anchor;
-    
-    btDeformableNodeAnchorConstraint(){}
-    btDeformableNodeAnchorConstraint(const btSoftBody::DeformableNodeRigidAnchor& c);
+	
+    btDeformableNodeAnchorConstraint(const btSoftBody::DeformableNodeRigidAnchor& c, const btContactSolverInfo& infoGlobal);
     btDeformableNodeAnchorConstraint(const btDeformableNodeAnchorConstraint& other);
+	btDeformableNodeAnchorConstraint(){}
     virtual ~btDeformableNodeAnchorConstraint()
     {
     }
-    virtual btScalar solveConstraint();
-    virtual btScalar solveSplitImpulse(const btContactSolverInfo& infoGlobal)
-    {
-        // todo xuchenhan@
-        return 0;
-    }
+    virtual btScalar solveConstraint(const btContactSolverInfo& infoGlobal);
+
     // object A is the rigid/multi body, and object B is the deformable node/face
     virtual btVector3 getVa() const;
     // get the velocity of the deformable node in contact
@@ -152,10 +135,7 @@ public:
         return btVector3(0,0,0);
     }
     virtual void applyImpulse(const btVector3& impulse);
-    virtual void applySplitImpulse(const btVector3& impulse)
-    {
-        // todo xuchenhan@
-    };
+
     virtual void setPenetrationScale(btScalar scale){}
 };
 
@@ -169,10 +149,10 @@ public:
     btVector3 m_total_tangent_dv;
     btScalar m_penetration;
     const btSoftBody::DeformableRigidContact* m_contact;
-    
-    btDeformableRigidContactConstraint(){}
-    btDeformableRigidContactConstraint(const btSoftBody::DeformableRigidContact& c);
+	
+    btDeformableRigidContactConstraint(const btSoftBody::DeformableRigidContact& c, const btContactSolverInfo& infoGlobal);
     btDeformableRigidContactConstraint(const btDeformableRigidContactConstraint& other);
+	btDeformableRigidContactConstraint(){}
     virtual ~btDeformableRigidContactConstraint()
     {
     }
@@ -180,9 +160,7 @@ public:
     // object A is the rigid/multi body, and object B is the deformable node/face
     virtual btVector3 getVa() const;
     
-    virtual btScalar solveConstraint();
-    
-    virtual btScalar solveSplitImpulse(const btContactSolverInfo& infoGlobal);
+    virtual btScalar solveConstraint(const btContactSolverInfo& infoGlobal);
     
     virtual void setPenetrationScale(btScalar scale)
     {
@@ -196,12 +174,11 @@ class btDeformableNodeRigidContactConstraint : public btDeformableRigidContactCo
 {
 public:
     // the deformable node in contact
-    const btSoftBody::Node* m_node;
-    
-    btDeformableNodeRigidContactConstraint(){}
-    btDeformableNodeRigidContactConstraint(const btSoftBody::DeformableNodeRigidContact& contact);
+    btSoftBody::Node* m_node;
+	
+    btDeformableNodeRigidContactConstraint(const btSoftBody::DeformableNodeRigidContact& contact, const btContactSolverInfo& infoGlobal);
     btDeformableNodeRigidContactConstraint(const btDeformableNodeRigidContactConstraint& other);
-    
+	btDeformableNodeRigidContactConstraint(){}
     virtual ~btDeformableNodeRigidContactConstraint()
     {
     }
@@ -219,7 +196,6 @@ public:
     }
     
     virtual void applyImpulse(const btVector3& impulse);
-    virtual void applySplitImpulse(const btVector3& impulse);
 };
 
 //
@@ -228,10 +204,10 @@ class btDeformableFaceRigidContactConstraint : public btDeformableRigidContactCo
 {
 public:
     const btSoftBody::Face* m_face;
-    btDeformableFaceRigidContactConstraint(){}
-    btDeformableFaceRigidContactConstraint(const btSoftBody::DeformableFaceRigidContact& contact);
+    bool m_useStrainLimiting;
+    btDeformableFaceRigidContactConstraint(const btSoftBody::DeformableFaceRigidContact& contact, const btContactSolverInfo& infoGlobal, bool useStrainLimiting);
     btDeformableFaceRigidContactConstraint(const btDeformableFaceRigidContactConstraint& other);
-    
+    btDeformableFaceRigidContactConstraint(): m_useStrainLimiting(false) {}
     virtual ~btDeformableFaceRigidContactConstraint()
     {
     }
@@ -249,7 +225,6 @@ public:
     }
     
     virtual void applyImpulse(const btVector3& impulse);
-    virtual void applySplitImpulse(const btVector3& impulse);
 };
 
 //
@@ -263,19 +238,11 @@ public:
     btVector3 m_total_normal_dv;
     btVector3 m_total_tangent_dv;
     
-    btDeformableFaceNodeContactConstraint(){}
-    
-    btDeformableFaceNodeContactConstraint(const btSoftBody::DeformableFaceNodeContact& contact);
-    
+    btDeformableFaceNodeContactConstraint(const btSoftBody::DeformableFaceNodeContact& contact, const btContactSolverInfo& infoGlobal);
+	btDeformableFaceNodeContactConstraint(){}
     virtual ~btDeformableFaceNodeContactConstraint(){}
     
-    virtual btScalar solveConstraint();
-    
-    virtual btScalar solveSplitImpulse(const btContactSolverInfo& infoGlobal)
-    {
-        // todo: xuchenhan@
-        return 0;
-    }
+    virtual btScalar solveConstraint(const btContactSolverInfo& infoGlobal);
     
     // get the velocity of the object A in the contact
     virtual btVector3 getVa() const;
@@ -293,10 +260,7 @@ public:
     }
     
     virtual void applyImpulse(const btVector3& impulse);
-    virtual void applySplitImpulse(const btVector3& impulse)
-    {
-        // todo xuchenhan@
-    }
+
     virtual void setPenetrationScale(btScalar scale){}
 };
 #endif /* BT_DEFORMABLE_CONTACT_CONSTRAINT_H */
diff --git a/thirdparty/bullet/BulletSoftBody/btDeformableContactProjection.cpp b/thirdparty/bullet/BulletSoftBody/btDeformableContactProjection.cpp
index 5a4f3241b4..22ca8bf582 100644
--- a/thirdparty/bullet/BulletSoftBody/btDeformableContactProjection.cpp
+++ b/thirdparty/bullet/BulletSoftBody/btDeformableContactProjection.cpp
@@ -17,7 +17,7 @@
 #include "btDeformableMultiBodyDynamicsWorld.h"
 #include <algorithm>
 #include <cmath>
-btScalar btDeformableContactProjection::update(btCollisionObject** deformableBodies,int numDeformableBodies)
+btScalar btDeformableContactProjection::update(btCollisionObject** deformableBodies,int numDeformableBodies, const btContactSolverInfo& infoGlobal)
 {
 	btScalar residualSquare = 0;
 	for (int i = 0; i < numDeformableBodies; ++i)
@@ -32,25 +32,25 @@ btScalar btDeformableContactProjection::update(btCollisionObject** deformableBod
 			for (int k = 0; k < m_nodeRigidConstraints[j].size(); ++k)
 			{
 				btDeformableNodeRigidContactConstraint& constraint = m_nodeRigidConstraints[j][k];
-				btScalar localResidualSquare = constraint.solveConstraint();
+				btScalar localResidualSquare = constraint.solveConstraint(infoGlobal);
 				residualSquare = btMax(residualSquare, localResidualSquare);
 			}
 			for (int k = 0; k < m_nodeAnchorConstraints[j].size(); ++k)
 			{
 				btDeformableNodeAnchorConstraint& constraint = m_nodeAnchorConstraints[j][k];
-				btScalar localResidualSquare = constraint.solveConstraint();
+				btScalar localResidualSquare = constraint.solveConstraint(infoGlobal);
 				residualSquare = btMax(residualSquare, localResidualSquare);
 			}
 			for (int k = 0; k < m_faceRigidConstraints[j].size(); ++k)
 			{
 				btDeformableFaceRigidContactConstraint& constraint = m_faceRigidConstraints[j][k];
-				btScalar localResidualSquare = constraint.solveConstraint();
+				btScalar localResidualSquare = constraint.solveConstraint(infoGlobal);
 				residualSquare = btMax(residualSquare, localResidualSquare);
 			}
 			for (int k = 0; k < m_deformableConstraints[j].size(); ++k)
 			{
 				btDeformableFaceNodeContactConstraint& constraint = m_deformableConstraints[j][k];
-				btScalar localResidualSquare = constraint.solveConstraint();
+				btScalar localResidualSquare = constraint.solveConstraint(infoGlobal);
 				residualSquare = btMax(residualSquare, localResidualSquare);
 			}
 		}
@@ -77,39 +77,8 @@ void btDeformableContactProjection::splitImpulseSetup(const btContactSolverInfo&
 	}
 }
 
-btScalar btDeformableContactProjection::solveSplitImpulse(const btContactSolverInfo& infoGlobal)
-{
-	btScalar residualSquare = 0;
-	for (int i = 0; i < m_softBodies.size(); ++i)
-	{
-		// node constraints
-		for (int j = 0; j < m_nodeRigidConstraints[i].size(); ++j)
-		{
-			btDeformableNodeRigidContactConstraint& constraint = m_nodeRigidConstraints[i][j];
-			btScalar localResidualSquare = constraint.solveSplitImpulse(infoGlobal);
-			residualSquare = btMax(residualSquare, localResidualSquare);
-		}
-		// anchor constraints
-		for (int j = 0; j < m_nodeAnchorConstraints[i].size(); ++j)
-		{
-			btDeformableNodeAnchorConstraint& constraint = m_nodeAnchorConstraints[i][j];
-			btScalar localResidualSquare = constraint.solveSplitImpulse(infoGlobal);
-			residualSquare = btMax(residualSquare, localResidualSquare);
-		}
-		// face constraints
-		for (int j = 0; j < m_faceRigidConstraints[i].size(); ++j)
-		{
-			btDeformableFaceRigidContactConstraint& constraint = m_faceRigidConstraints[i][j];
-			btScalar localResidualSquare = constraint.solveSplitImpulse(infoGlobal);
-			residualSquare = btMax(residualSquare, localResidualSquare);
-		}
-
-	}
-	return residualSquare;
-}
-
-void btDeformableContactProjection::setConstraints()
-{
+void btDeformableContactProjection::setConstraints(const btContactSolverInfo& infoGlobal)
+{  
 	BT_PROFILE("setConstraints");
 	for (int i = 0; i < m_softBodies.size(); ++i)
 	{
@@ -124,7 +93,7 @@ void btDeformableContactProjection::setConstraints()
 		{
 			if (psb->m_nodes[j].m_im == 0)
 			{
-				btDeformableStaticConstraint static_constraint(&psb->m_nodes[j]);
+				btDeformableStaticConstraint static_constraint(&psb->m_nodes[j], infoGlobal);
 				m_staticConstraints[i].push_back(static_constraint);
 			}
 		}
@@ -139,7 +108,7 @@ void btDeformableContactProjection::setConstraints()
 				continue;
 			}
 			anchor.m_c1 = anchor.m_cti.m_colObj->getWorldTransform().getBasis() * anchor.m_local;
-			btDeformableNodeAnchorConstraint constraint(anchor);
+			btDeformableNodeAnchorConstraint constraint(anchor, infoGlobal);
 			m_nodeAnchorConstraints[i].push_back(constraint);
 		}
 		
@@ -152,7 +121,7 @@ void btDeformableContactProjection::setConstraints()
 			{
 				continue;
 			}
-			btDeformableNodeRigidContactConstraint constraint(contact);
+			btDeformableNodeRigidContactConstraint constraint(contact, infoGlobal);
 			btVector3 va = constraint.getVa();
 			btVector3 vb = constraint.getVb();
 			const btVector3 vr = vb - va;
@@ -173,7 +142,7 @@ void btDeformableContactProjection::setConstraints()
 			{
 				continue;
 			}
-			btDeformableFaceRigidContactConstraint constraint(contact);
+			btDeformableFaceRigidContactConstraint constraint(contact, infoGlobal, m_useStrainLimiting);
 			btVector3 va = constraint.getVa();
 			btVector3 vb = constraint.getVb();
 			const btVector3 vr = vb - va;
@@ -184,253 +153,404 @@ void btDeformableContactProjection::setConstraints()
 				m_faceRigidConstraints[i].push_back(constraint);
 			}
 		}
-		
-		// set Deformable Face vs. Deformable Node constraint
-		for (int j = 0; j < psb->m_faceNodeContacts.size(); ++j)
-		{
-			const btSoftBody::DeformableFaceNodeContact& contact = psb->m_faceNodeContacts[j];
-
-			btDeformableFaceNodeContactConstraint constraint(contact);
-			btVector3 va = constraint.getVa();
-			btVector3 vb = constraint.getVb();
-			const btVector3 vr = vb - va;
-			const btScalar dn = btDot(vr, contact.m_normal);
-			if (dn > -SIMD_EPSILON)
-			{
-				m_deformableConstraints[i].push_back(constraint);
-			}
-		}
 	}
 }
 
 void btDeformableContactProjection::project(TVStack& x)
 {
-	const int dim = 3;
-	for (int index = 0; index < m_projectionsDict.size(); ++index)
-	{
-		btAlignedObjectArray<btVector3>& projectionDirs = *m_projectionsDict.getAtIndex(index);
-		size_t i = m_projectionsDict.getKeyAtIndex(index).getUid1();
-		if (projectionDirs.size() >= dim)
-		{
-			// static node
-			x[i].setZero();
-			continue;
-		}
-		else if (projectionDirs.size() == 2)
-		{
-			btVector3 dir0 = projectionDirs[0];
-			btVector3 dir1 = projectionDirs[1];
-			btVector3 free_dir = btCross(dir0, dir1);
-			if (free_dir.safeNorm() < SIMD_EPSILON)
-			{
-				x[i] -= x[i].dot(dir0) * dir0;
-				x[i] -= x[i].dot(dir1) * dir1;
-			}
-			else
-			{
-				free_dir.normalize();
-				x[i] = x[i].dot(free_dir) * free_dir;
-			}
-		}
-		else
-		{
-			btAssert(projectionDirs.size() == 1);
-			btVector3 dir0 = projectionDirs[0];
-			x[i] -= x[i].dot(dir0) * dir0;
-		}
-	}
+#ifndef USE_MGS
+    const int dim = 3;
+    for (int index = 0; index < m_projectionsDict.size(); ++index)
+    {
+        btAlignedObjectArray<btVector3>& projectionDirs = *m_projectionsDict.getAtIndex(index);
+        size_t i = m_projectionsDict.getKeyAtIndex(index).getUid1();
+        if (projectionDirs.size() >= dim)
+        {
+            // static node
+            x[i].setZero();
+            continue;
+        }
+        else if (projectionDirs.size() == 2)
+        {
+            btVector3 dir0 = projectionDirs[0];
+            btVector3 dir1 = projectionDirs[1];
+            btVector3 free_dir = btCross(dir0, dir1);
+            if (free_dir.safeNorm() < SIMD_EPSILON)
+            {
+                x[i] -= x[i].dot(dir0) * dir0;
+                x[i] -= x[i].dot(dir1) * dir1;
+            }
+            else
+            {
+                free_dir.normalize();
+                x[i] = x[i].dot(free_dir) * free_dir;
+            }
+        }
+        else
+        {
+            btAssert(projectionDirs.size() == 1);
+            btVector3 dir0 = projectionDirs[0];
+            x[i] -= x[i].dot(dir0) * dir0;
+        }
+    }
+#else
+    btReducedVector p(x.size());
+    for (int i = 0; i < m_projections.size(); ++i)
+    {
+        p += (m_projections[i].dot(x) * m_projections[i]);
+    }
+    for (int i = 0; i < p.m_indices.size(); ++i)
+    {
+        x[p.m_indices[i]] -= p.m_vecs[i];
+    }
+#endif
 }
 
 void btDeformableContactProjection::setProjection()
 {
-	btAlignedObjectArray<btVector3> units;
-	units.push_back(btVector3(1,0,0));
-	units.push_back(btVector3(0,1,0));
-	units.push_back(btVector3(0,0,1));
-	for (int i = 0; i < m_softBodies.size(); ++i)
-	{
-		btSoftBody* psb = m_softBodies[i];
-		if (!psb->isActive())
-		{
-			continue;
-		}
-		for (int j = 0; j < m_staticConstraints[i].size(); ++j)
-		{
-			int index = m_staticConstraints[i][j].m_node->index;
-			if (m_projectionsDict.find(index) == NULL)
+#ifndef USE_MGS
+    BT_PROFILE("btDeformableContactProjection::setProjection");
+    btAlignedObjectArray<btVector3> units;
+    units.push_back(btVector3(1,0,0));
+    units.push_back(btVector3(0,1,0));
+    units.push_back(btVector3(0,0,1));
+    for (int i = 0; i < m_softBodies.size(); ++i)
+    {
+        btSoftBody* psb = m_softBodies[i];
+        if (!psb->isActive())
+        {
+            continue;
+        }
+        for (int j = 0; j < m_staticConstraints[i].size(); ++j)
+        {
+            int index = m_staticConstraints[i][j].m_node->index;
+            m_staticConstraints[i][j].m_node->m_penetration = SIMD_INFINITY;
+            if (m_projectionsDict.find(index) == NULL)
+            {
+                m_projectionsDict.insert(index, units);
+            }
+            else
+            {
+                btAlignedObjectArray<btVector3>& projections = *m_projectionsDict[index];
+                for (int k = 0; k < 3; ++k)
+                {
+                    projections.push_back(units[k]);
+                }
+            }
+        }
+        for (int j = 0; j < m_nodeAnchorConstraints[i].size(); ++j)
+        {
+            int index = m_nodeAnchorConstraints[i][j].m_anchor->m_node->index;
+            m_nodeAnchorConstraints[i][j].m_anchor->m_node->m_penetration = SIMD_INFINITY;
+            if (m_projectionsDict.find(index) == NULL)
+            {
+                m_projectionsDict.insert(index, units);
+            }
+            else
+            {
+                btAlignedObjectArray<btVector3>& projections = *m_projectionsDict[index];
+                for (int k = 0; k < 3; ++k)
+                {
+                    projections.push_back(units[k]);
+                }
+            }
+        }
+        for (int j = 0; j < m_nodeRigidConstraints[i].size(); ++j)
+        {
+            int index = m_nodeRigidConstraints[i][j].m_node->index;
+            m_nodeRigidConstraints[i][j].m_node->m_penetration = -m_nodeRigidConstraints[i][j].getContact()->m_cti.m_offset;
+            if (m_nodeRigidConstraints[i][j].m_static)
+            {
+                if (m_projectionsDict.find(index) == NULL)
+                {
+                    m_projectionsDict.insert(index, units);
+                }
+                else
+                {
+                    btAlignedObjectArray<btVector3>& projections = *m_projectionsDict[index];
+                    for (int k = 0; k < 3; ++k)
+                    {
+                        projections.push_back(units[k]);
+                    }
+                }
+            }
+            else
+            {
+                if (m_projectionsDict.find(index) == NULL)
+                {
+                    btAlignedObjectArray<btVector3> projections;
+                    projections.push_back(m_nodeRigidConstraints[i][j].m_normal);
+                    m_projectionsDict.insert(index, projections);
+                }
+                else
+                {
+                    btAlignedObjectArray<btVector3>& projections = *m_projectionsDict[index];
+                    projections.push_back(m_nodeRigidConstraints[i][j].m_normal);
+                }
+            }
+        }
+        for (int j = 0; j < m_faceRigidConstraints[i].size(); ++j)
+        {
+            const btSoftBody::Face* face = m_faceRigidConstraints[i][j].m_face;
+            btScalar penetration = -m_faceRigidConstraints[i][j].getContact()->m_cti.m_offset;
+            for (int k = 0; k < 3; ++k)
+            {
+                face->m_n[k]->m_penetration = btMax(face->m_n[k]->m_penetration, penetration);
+            }
+            for (int k = 0; k < 3; ++k)
+            {
+                btSoftBody::Node* node = face->m_n[k];
+                node->m_penetration = true;
+                int index = node->index;
+                if (m_faceRigidConstraints[i][j].m_static)
+                {
+                    if (m_projectionsDict.find(index) == NULL)
+                    {
+                        m_projectionsDict.insert(index, units);
+                    }
+                    else
+                    {
+                        btAlignedObjectArray<btVector3>& projections = *m_projectionsDict[index];
+                        for (int k = 0; k < 3; ++k)
+                        {
+                            projections.push_back(units[k]);
+                        }
+                    }
+                }
+                else
+                {
+                    if (m_projectionsDict.find(index) == NULL)
+                    {
+                        btAlignedObjectArray<btVector3> projections;
+                        projections.push_back(m_faceRigidConstraints[i][j].m_normal);
+                        m_projectionsDict.insert(index, projections);
+                    }
+                    else
+                    {
+                        btAlignedObjectArray<btVector3>& projections = *m_projectionsDict[index];
+                        projections.push_back(m_faceRigidConstraints[i][j].m_normal);
+                    }
+                }
+            }
+        }
+    }
+#else
+    int dof = 0;
+    for (int i = 0; i < m_softBodies.size(); ++i)
+    {
+        dof += m_softBodies[i]->m_nodes.size();
+    }
+    for (int i = 0; i < m_softBodies.size(); ++i)
+    {
+        btSoftBody* psb = m_softBodies[i];
+        if (!psb->isActive())
+        {
+            continue;
+        }
+        for (int j = 0; j < m_staticConstraints[i].size(); ++j)
+        {
+            int index = m_staticConstraints[i][j].m_node->index;
+            m_staticConstraints[i][j].m_node->m_penetration = SIMD_INFINITY;
+            btAlignedObjectArray<int> indices;
+            btAlignedObjectArray<btVector3> vecs1,vecs2,vecs3;
+            indices.push_back(index);
+            vecs1.push_back(btVector3(1,0,0));
+            vecs2.push_back(btVector3(0,1,0));
+            vecs3.push_back(btVector3(0,0,1));
+            m_projections.push_back(btReducedVector(dof, indices, vecs1));
+            m_projections.push_back(btReducedVector(dof, indices, vecs2));
+            m_projections.push_back(btReducedVector(dof, indices, vecs3));
+        }
+        
+        for (int j = 0; j < m_nodeAnchorConstraints[i].size(); ++j)
+        {
+            int index = m_nodeAnchorConstraints[i][j].m_anchor->m_node->index;
+            m_nodeAnchorConstraints[i][j].m_anchor->m_node->m_penetration = SIMD_INFINITY;
+            btAlignedObjectArray<int> indices;
+            btAlignedObjectArray<btVector3> vecs1,vecs2,vecs3;
+            indices.push_back(index);
+            vecs1.push_back(btVector3(1,0,0));
+            vecs2.push_back(btVector3(0,1,0));
+            vecs3.push_back(btVector3(0,0,1));
+            m_projections.push_back(btReducedVector(dof, indices, vecs1));
+            m_projections.push_back(btReducedVector(dof, indices, vecs2));
+            m_projections.push_back(btReducedVector(dof, indices, vecs3));
+        }
+        for (int j = 0; j < m_nodeRigidConstraints[i].size(); ++j)
+        {
+            int index = m_nodeRigidConstraints[i][j].m_node->index;
+            m_nodeRigidConstraints[i][j].m_node->m_penetration = -m_nodeRigidConstraints[i][j].getContact()->m_cti.m_offset;
+            btAlignedObjectArray<int> indices;
+            indices.push_back(index);
+            btAlignedObjectArray<btVector3> vecs1,vecs2,vecs3;
+            if (m_nodeRigidConstraints[i][j].m_static)
+            {
+                vecs1.push_back(btVector3(1,0,0));
+                vecs2.push_back(btVector3(0,1,0));
+                vecs3.push_back(btVector3(0,0,1));
+                m_projections.push_back(btReducedVector(dof, indices, vecs1));
+                m_projections.push_back(btReducedVector(dof, indices, vecs2));
+                m_projections.push_back(btReducedVector(dof, indices, vecs3));
+            }
+            else
+            {
+                vecs1.push_back(m_nodeRigidConstraints[i][j].m_normal);
+                m_projections.push_back(btReducedVector(dof, indices, vecs1));
+            }
+        }
+        for (int j = 0; j < m_faceRigidConstraints[i].size(); ++j)
+        {
+            const btSoftBody::Face* face = m_faceRigidConstraints[i][j].m_face;
+			btVector3 bary = m_faceRigidConstraints[i][j].getContact()->m_bary;
+            btScalar penetration = -m_faceRigidConstraints[i][j].getContact()->m_cti.m_offset;
+            for (int k = 0; k < 3; ++k)
+            {
+                face->m_n[k]->m_penetration = btMax(face->m_n[k]->m_penetration, penetration);
+            }
+			if (m_faceRigidConstraints[i][j].m_static)
 			{
-				m_projectionsDict.insert(index, units);
-			}
-			else
-			{
-				btAlignedObjectArray<btVector3>& projections = *m_projectionsDict[index];
-				for (int k = 0; k < 3; ++k)
+				for (int l = 0; l < 3; ++l)
 				{
-					projections.push_back(units[k]);
-				}
-			}
-		}
-		for (int j = 0; j < m_nodeAnchorConstraints[i].size(); ++j)
-		{
-			int index = m_nodeAnchorConstraints[i][j].m_anchor->m_node->index;
-			if (m_projectionsDict.find(index) == NULL)
-			{
-				m_projectionsDict.insert(index, units);
-			}
-			else
-			{
-				btAlignedObjectArray<btVector3>& projections = *m_projectionsDict[index];
-				for (int k = 0; k < 3; ++k)
-				{
-					projections.push_back(units[k]);
-				}
-			}
-		}
-		for (int j = 0; j < m_nodeRigidConstraints[i].size(); ++j)
-		{
-			int index = m_nodeRigidConstraints[i][j].m_node->index;
-			if (m_nodeRigidConstraints[i][j].m_static)
-			{
-				if (m_projectionsDict.find(index) == NULL)
-				{
-					m_projectionsDict.insert(index, units);
-				}
-				else
-				{
-					btAlignedObjectArray<btVector3>& projections = *m_projectionsDict[index];
+					
+					btReducedVector rv(dof);
 					for (int k = 0; k < 3; ++k)
 					{
-						projections.push_back(units[k]);
+						rv.m_indices.push_back(face->m_n[k]->index);
+						btVector3 v(0,0,0);
+						v[l] = bary[k];
+						rv.m_vecs.push_back(v);
+                        rv.sort();
 					}
+					m_projections.push_back(rv);
 				}
 			}
 			else
 			{
-				if (m_projectionsDict.find(index) == NULL)
-				{
-					btAlignedObjectArray<btVector3> projections;
-					projections.push_back(m_nodeRigidConstraints[i][j].m_normal);
-					m_projectionsDict.insert(index, projections);
-				}
-				else
-				{
-					btAlignedObjectArray<btVector3>& projections = *m_projectionsDict[index];
-					projections.push_back(m_nodeRigidConstraints[i][j].m_normal);
-				}
-			}
-		}
-		for (int j = 0; j < m_faceRigidConstraints[i].size(); ++j)
-		{
-			const btSoftBody::Face* face = m_faceRigidConstraints[i][j].m_face;
-			for (int k = 0; k < 3; ++k)
-			{
-				const btSoftBody::Node* node = face->m_n[k];
-				int index = node->index;
-				if (m_faceRigidConstraints[i][j].m_static)
+				btReducedVector rv(dof);
+				for (int k = 0; k < 3; ++k)
 				{
-					if (m_projectionsDict.find(index) == NULL)
-					{
-						m_projectionsDict.insert(index, units);
-					}
-					else
-					{
-						btAlignedObjectArray<btVector3>& projections = *m_projectionsDict[index];
-						for (int k = 0; k < 3; ++k)
-						{
-							projections.push_back(units[k]);
-						}
-					}
-				}
-				else
-				{
-					if (m_projectionsDict.find(index) == NULL)
-					{
-						btAlignedObjectArray<btVector3> projections;
-						projections.push_back(m_faceRigidConstraints[i][j].m_normal);
-						m_projectionsDict.insert(index, projections);
-					}
-					else
-					{
-						btAlignedObjectArray<btVector3>& projections = *m_projectionsDict[index];
-						projections.push_back(m_faceRigidConstraints[i][j].m_normal);
-					}
+					rv.m_indices.push_back(face->m_n[k]->index);
+					rv.m_vecs.push_back(bary[k] * m_faceRigidConstraints[i][j].m_normal);
+                    rv.sort();
 				}
+				m_projections.push_back(rv);
 			}
 		}
-		for (int j = 0; j < m_deformableConstraints[i].size(); ++j)
-		{
-			const btSoftBody::Face* face = m_deformableConstraints[i][j].m_face;
-			for (int k = 0; k < 3; ++k)
-			{
-				const btSoftBody::Node* node = face->m_n[k];
-				int index = node->index;
-				if (m_deformableConstraints[i][j].m_static)
-				{
-					if (m_projectionsDict.find(index) == NULL)
-					{
-						m_projectionsDict.insert(index, units);
-					}
-					else
-					{
-						btAlignedObjectArray<btVector3>& projections = *m_projectionsDict[index];
-						for (int k = 0; k < 3; ++k)
-						{
-							projections.push_back(units[k]);
-						}
-					}
-				}
-				else
-				{
-					if (m_projectionsDict.find(index) == NULL)
-					{
-						btAlignedObjectArray<btVector3> projections;
-						projections.push_back(m_deformableConstraints[i][j].m_normal);
-						m_projectionsDict.insert(index, projections);
-					}
-					else
-					{
-						btAlignedObjectArray<btVector3>& projections = *m_projectionsDict[index];
-						projections.push_back(m_deformableConstraints[i][j].m_normal);
-					}
-				}
-			}
+    }
+    btModifiedGramSchmidt<btReducedVector> mgs(m_projections);
+    mgs.solve();
+    m_projections = mgs.m_out;
+#endif
+}
+
+void btDeformableContactProjection::checkConstraints(const TVStack& x)
+{
+    for (int i = 0; i < m_lagrangeMultipliers.size(); ++i)
+    {
+        btVector3 d(0,0,0);
+        const LagrangeMultiplier& lm = m_lagrangeMultipliers[i];
+        for (int j = 0; j < lm.m_num_constraints; ++j)
+        {
+            for (int k = 0; k < lm.m_num_nodes; ++k)
+            {
+                d[j] += lm.m_weights[k] * x[lm.m_indices[k]].dot(lm.m_dirs[j]);
+            }
+        }
+        printf("d = %f, %f, %f\n",d[0],d[1],d[2]);
+    }
+}
+
+void btDeformableContactProjection::setLagrangeMultiplier()
+{
+    for (int i = 0; i < m_softBodies.size(); ++i)
+    {
+        btSoftBody* psb = m_softBodies[i];
+        if (!psb->isActive())
+        {
+            continue;
+        }
+        for (int j = 0; j < m_staticConstraints[i].size(); ++j)
+        {
+            int index = m_staticConstraints[i][j].m_node->index;
+            m_staticConstraints[i][j].m_node->m_penetration = SIMD_INFINITY;
+            LagrangeMultiplier lm;
+            lm.m_num_nodes = 1;
+            lm.m_indices[0] = index;
+            lm.m_weights[0] = 1.0;
+            lm.m_num_constraints = 3;
+            lm.m_dirs[0] = btVector3(1,0,0);
+            lm.m_dirs[1] = btVector3(0,1,0);
+            lm.m_dirs[2] = btVector3(0,0,1);
+            m_lagrangeMultipliers.push_back(lm);
+        }
+        for (int j = 0; j < m_nodeAnchorConstraints[i].size(); ++j)
+        {
+            int index = m_nodeAnchorConstraints[i][j].m_anchor->m_node->index;
+            m_nodeAnchorConstraints[i][j].m_anchor->m_node->m_penetration = SIMD_INFINITY;
+            LagrangeMultiplier lm;
+            lm.m_num_nodes = 1;
+            lm.m_indices[0] = index;
+            lm.m_weights[0] = 1.0;
+            lm.m_num_constraints = 3;
+            lm.m_dirs[0] = btVector3(1,0,0);
+            lm.m_dirs[1] = btVector3(0,1,0);
+            lm.m_dirs[2] = btVector3(0,0,1);
+            m_lagrangeMultipliers.push_back(lm);
+        }
+        for (int j = 0; j < m_nodeRigidConstraints[i].size(); ++j)
+        {
+            int index = m_nodeRigidConstraints[i][j].m_node->index;
+            m_nodeRigidConstraints[i][j].m_node->m_penetration = -m_nodeRigidConstraints[i][j].getContact()->m_cti.m_offset;
+            LagrangeMultiplier lm;
+            lm.m_num_nodes = 1;
+            lm.m_indices[0] = index;
+            lm.m_weights[0] = 1.0;
+            if (m_nodeRigidConstraints[i][j].m_static)
+            {
+                lm.m_num_constraints = 3;
+                lm.m_dirs[0] = btVector3(1,0,0);
+                lm.m_dirs[1] = btVector3(0,1,0);
+                lm.m_dirs[2] = btVector3(0,0,1);
+            }
+            else
+            {
+                lm.m_num_constraints = 1;
+                lm.m_dirs[0] = m_nodeRigidConstraints[i][j].m_normal;
+            }
+            m_lagrangeMultipliers.push_back(lm);
+        }
+        for (int j = 0; j < m_faceRigidConstraints[i].size(); ++j)
+        {
+            const btSoftBody::Face* face = m_faceRigidConstraints[i][j].m_face;
 			
-			const btSoftBody::Node* node = m_deformableConstraints[i][j].m_node;
-			int index = node->index;
-			if (m_deformableConstraints[i][j].m_static)
+            btVector3 bary = m_faceRigidConstraints[i][j].getContact()->m_bary;
+            btScalar penetration = -m_faceRigidConstraints[i][j].getContact()->m_cti.m_offset;
+			LagrangeMultiplier lm;
+			lm.m_num_nodes = 3;
+			for (int k = 0; k<3; ++k)
 			{
-				if (m_projectionsDict.find(index) == NULL)
-				{
-					m_projectionsDict.insert(index, units);
-				}
-				else
-				{
-					btAlignedObjectArray<btVector3>& projections = *m_projectionsDict[index];
-					for (int k = 0; k < 3; ++k)
-					{
-						projections.push_back(units[k]);
-					}
-				}
+				face->m_n[k]->m_penetration = btMax(face->m_n[k]->m_penetration, penetration);
+				lm.m_indices[k] = face->m_n[k]->index;
+				lm.m_weights[k] = bary[k];
+			}
+            if (m_faceRigidConstraints[i][j].m_static)
+            {
+				lm.m_num_constraints = 3;
+				lm.m_dirs[0] = btVector3(1,0,0);
+				lm.m_dirs[1] = btVector3(0,1,0);
+				lm.m_dirs[2] = btVector3(0,0,1);
 			}
 			else
 			{
-				if (m_projectionsDict.find(index) == NULL)
-				{
-					btAlignedObjectArray<btVector3> projections;
-					projections.push_back(m_deformableConstraints[i][j].m_normal);
-					m_projectionsDict.insert(index, projections);
-				}
-				else
-				{
-					btAlignedObjectArray<btVector3>& projections = *m_projectionsDict[index];
-					projections.push_back(m_deformableConstraints[i][j].m_normal);
-				}
+				lm.m_num_constraints = 1;
+				lm.m_dirs[0] = m_faceRigidConstraints[i][j].m_normal;
 			}
+            m_lagrangeMultipliers.push_back(lm);
 		}
 	}
 }
 
-
+//
 void btDeformableContactProjection::applyDynamicFriction(TVStack& f)
 {
 	for (int i = 0; i < m_softBodies.size(); ++i)
@@ -502,7 +622,12 @@ void btDeformableContactProjection::reinitialize(bool nodeUpdated)
 		m_faceRigidConstraints[i].clear();
 		m_deformableConstraints[i].clear();
 	}
-	m_projectionsDict.clear();
+#ifndef USE_MGS
+    m_projectionsDict.clear();
+#else
+    m_projections.clear();
+#endif
+    m_lagrangeMultipliers.clear();
 }
 
 
diff --git a/thirdparty/bullet/BulletSoftBody/btDeformableContactProjection.h b/thirdparty/bullet/BulletSoftBody/btDeformableContactProjection.h
index 3c4490765e..8d7e94d4fb 100644
--- a/thirdparty/bullet/BulletSoftBody/btDeformableContactProjection.h
+++ b/thirdparty/bullet/BulletSoftBody/btDeformableContactProjection.h
@@ -21,30 +21,37 @@
 #include "BulletDynamics/Featherstone/btMultiBodyConstraint.h"
 #include "btDeformableContactConstraint.h"
 #include "LinearMath/btHashMap.h"
+#include "LinearMath/btReducedVector.h"
+#include "LinearMath/btModifiedGramSchmidt.h"
 #include <vector>
+
+struct LagrangeMultiplier
+{
+    int m_num_constraints;        // Number of constraints
+    int m_num_nodes;              // Number of nodes in these constraints
+    btScalar m_weights[3];        // weights of the nodes involved, same size as m_num_nodes
+    btVector3 m_dirs[3];          // Constraint directions, same size of m_num_constraints;
+    int m_indices[3];             // indices of the nodes involved, same size as m_num_nodes;
+};
+
+
 class btDeformableContactProjection
 {
 public:
     typedef btAlignedObjectArray<btVector3> TVStack;
     btAlignedObjectArray<btSoftBody *>& m_softBodies;
-    
-//    // map from node index to static constraint
-//    btHashMap<btHashInt, btDeformableStaticConstraint> m_staticConstraints;
-//    // map from node index to node rigid constraint
-//    btHashMap<btHashInt, btAlignedObjectArray<btDeformableNodeRigidContactConstraint> > m_nodeRigidConstraints;
-//    // map from node index to face rigid constraint
-//    btHashMap<btHashInt, btAlignedObjectArray<btDeformableFaceRigidContactConstraint*> > m_faceRigidConstraints;
-//    // map from node index to deformable constraint
-//    btHashMap<btHashInt, btAlignedObjectArray<btDeformableFaceNodeContactConstraint*> > m_deformableConstraints;
-//    // map from node index to node anchor constraint
-//    btHashMap<btHashInt, btDeformableNodeAnchorConstraint> m_nodeAnchorConstraints;
 	
     // all constraints involving face
     btAlignedObjectArray<btDeformableContactConstraint*> m_allFaceConstraints;
-    
+#ifndef USE_MGS
     // map from node index to projection directions
     btHashMap<btHashInt, btAlignedObjectArray<btVector3> > m_projectionsDict;
-	
+#else
+    btAlignedObjectArray<btReducedVector> m_projections;
+#endif
+    
+    btAlignedObjectArray<LagrangeMultiplier> m_lagrangeMultipliers;
+    
 	// map from node index to static constraint
 	btAlignedObjectArray<btAlignedObjectArray<btDeformableStaticConstraint> > m_staticConstraints;
 	// map from node index to node rigid constraint
@@ -56,6 +63,8 @@ public:
 	// map from node index to node anchor constraint
 	btAlignedObjectArray<btAlignedObjectArray<btDeformableNodeAnchorConstraint> > m_nodeAnchorConstraints;
     
+    bool m_useStrainLimiting;
+    
     btDeformableContactProjection(btAlignedObjectArray<btSoftBody *>& softBodies)
     : m_softBodies(softBodies)
     {
@@ -72,13 +81,10 @@ public:
     virtual void applyDynamicFriction(TVStack& f);
     
     // update and solve the constraints
-    virtual btScalar update(btCollisionObject** deformableBodies,int numDeformableBodies);
-    
-    // solve the position error using split impulse
-    virtual btScalar solveSplitImpulse(const btContactSolverInfo& infoGlobal);
+    virtual btScalar update(btCollisionObject** deformableBodies,int numDeformableBodies, const btContactSolverInfo& infoGlobal);
     
     // Add constraints to m_constraints. In addition, the constraints that each vertex own are recorded in m_constraintsDict.
-    virtual void setConstraints();
+    virtual void setConstraints(const btContactSolverInfo& infoGlobal);
     
     // Set up projections for each vertex by adding the projection direction to
     virtual void setProjection();
@@ -86,5 +92,9 @@ public:
     virtual void reinitialize(bool nodeUpdated);
     
     virtual void splitImpulseSetup(const btContactSolverInfo& infoGlobal);
+    
+    virtual void setLagrangeMultiplier();
+    
+    void checkConstraints(const TVStack& x);
 };
 #endif /* btDeformableContactProjection_h */
diff --git a/thirdparty/bullet/BulletSoftBody/btDeformableCorotatedForce.h b/thirdparty/bullet/BulletSoftBody/btDeformableCorotatedForce.h
index c2a26338e7..2d042df729 100644
--- a/thirdparty/bullet/BulletSoftBody/btDeformableCorotatedForce.h
+++ b/thirdparty/bullet/BulletSoftBody/btDeformableCorotatedForce.h
@@ -114,6 +114,8 @@ public:
     {
     }
     
+    virtual void buildDampingForceDifferentialDiagonal(btScalar scale, TVStack& diagA){}
+    
     virtual btDeformableLagrangianForceType getForceType()
     {
         return BT_COROTATED_FORCE;
diff --git a/thirdparty/bullet/BulletSoftBody/btDeformableGravityForce.h b/thirdparty/bullet/BulletSoftBody/btDeformableGravityForce.h
index 33e5a8564a..13ee3eacb6 100644
--- a/thirdparty/bullet/BulletSoftBody/btDeformableGravityForce.h
+++ b/thirdparty/bullet/BulletSoftBody/btDeformableGravityForce.h
@@ -50,6 +50,8 @@ public:
     {
     }
     
+    virtual void buildDampingForceDifferentialDiagonal(btScalar scale, TVStack& diagA){}
+    
     virtual void addScaledGravityForce(btScalar scale, TVStack& force)
     {
         int numNodes = getNumNodes();
diff --git a/thirdparty/bullet/BulletSoftBody/btDeformableLagrangianForce.h b/thirdparty/bullet/BulletSoftBody/btDeformableLagrangianForce.h
index 64e80e23b3..0b6447442d 100644
--- a/thirdparty/bullet/BulletSoftBody/btDeformableLagrangianForce.h
+++ b/thirdparty/bullet/BulletSoftBody/btDeformableLagrangianForce.h
@@ -26,7 +26,8 @@ enum btDeformableLagrangianForceType
     BT_MASSSPRING_FORCE = 2,
     BT_COROTATED_FORCE = 3,
     BT_NEOHOOKEAN_FORCE = 4,
-    BT_LINEAR_ELASTICITY_FORCE = 5
+    BT_LINEAR_ELASTICITY_FORCE = 5,
+    BT_MOUSE_PICKING_FORCE = 6
 };
 
 static inline double randomDouble(double low, double high)
@@ -53,6 +54,9 @@ public:
     // add damping df
     virtual void addScaledDampingForceDifferential(btScalar scale, const TVStack& dv, TVStack& df) = 0;
     
+    // build diagonal of A matrix
+    virtual void buildDampingForceDifferentialDiagonal(btScalar scale, TVStack& diagA) = 0;
+    
     // add elastic df
     virtual void addScaledElasticForceDifferential(btScalar scale, const TVStack& dx, TVStack& df) = 0;
     
@@ -85,6 +89,11 @@ public:
         m_softBodies.push_back(psb);
     }
     
+    virtual void removeSoftBody(btSoftBody* psb)
+    {
+        m_softBodies.remove(psb);
+    }
+    
     virtual void setIndices(const btAlignedObjectArray<btSoftBody::Node*>* nodes)
     {
         m_nodes = nodes;
diff --git a/thirdparty/bullet/BulletSoftBody/btDeformableMassSpringForce.h b/thirdparty/bullet/BulletSoftBody/btDeformableMassSpringForce.h
index 54b4e4481d..b128df92cc 100644
--- a/thirdparty/bullet/BulletSoftBody/btDeformableMassSpringForce.h
+++ b/thirdparty/bullet/BulletSoftBody/btDeformableMassSpringForce.h
@@ -149,6 +149,52 @@ public:
         }
     }
     
+    virtual void buildDampingForceDifferentialDiagonal(btScalar scale, TVStack& diagA)
+    {
+        // implicit damping force differential
+        for (int i = 0; i < m_softBodies.size(); ++i)
+        {
+            btSoftBody* psb = m_softBodies[i];
+            if (!psb->isActive())
+            {
+                continue;
+            }
+            btScalar scaled_k_damp = m_dampingStiffness * scale;
+            for (int j = 0; j < psb->m_links.size(); ++j)
+            {
+                const btSoftBody::Link& link = psb->m_links[j];
+                btSoftBody::Node* node1 = link.m_n[0];
+                btSoftBody::Node* node2 = link.m_n[1];
+                size_t id1 = node1->index;
+                size_t id2 = node2->index;
+                if (m_momentum_conserving)
+                {
+                    if ((node2->m_x - node1->m_x).norm() > SIMD_EPSILON)
+                    {
+                        btVector3 dir = (node2->m_x - node1->m_x).normalized();
+                        for (int d = 0; d < 3; ++d)
+                        {
+                            if (node1->m_im > 0)
+                                diagA[id1][d] -= scaled_k_damp * dir[d] * dir[d];
+                            if (node2->m_im > 0)
+                                diagA[id2][d] -= scaled_k_damp * dir[d] * dir[d];
+                        }
+                    }
+                }
+                else
+                {
+                    for (int d = 0; d < 3; ++d)
+                    {
+                        if (node1->m_im > 0)
+                            diagA[id1][d] -= scaled_k_damp;
+                        if (node2->m_im > 0)
+                            diagA[id2][d] -= scaled_k_damp;
+                    }
+                }
+            }
+        }
+    }
+    
     virtual double totalElasticEnergy(btScalar dt)
     {
         double energy = 0;
diff --git a/thirdparty/bullet/BulletSoftBody/btDeformableMousePickingForce.h b/thirdparty/bullet/BulletSoftBody/btDeformableMousePickingForce.h
new file mode 100644
index 0000000000..07c10935f4
--- /dev/null
+++ b/thirdparty/bullet/BulletSoftBody/btDeformableMousePickingForce.h
@@ -0,0 +1,145 @@
+/*
+ Written by Xuchen Han <xuchenhan2015@u.northwestern.edu>
+ 
+ Bullet Continuous Collision Detection and Physics Library
+ Copyright (c) 2019 Google Inc. http://bulletphysics.org
+ This software is provided 'as-is', without any express or implied warranty.
+ In no event will the authors be held liable for any damages arising from the use of this software.
+ Permission is granted to anyone to use this software for any purpose,
+ including commercial applications, and to alter it and redistribute it freely,
+ subject to the following restrictions:
+ 1. The origin of this software must not be misrepresented; you must not claim that you wrote the original software. If you use this software in a product, an acknowledgment in the product documentation would be appreciated but is not required.
+ 2. Altered source versions must be plainly marked as such, and must not be misrepresented as being the original software.
+ 3. This notice may not be removed or altered from any source distribution.
+ */
+
+#ifndef BT_MOUSE_PICKING_FORCE_H
+#define BT_MOUSE_PICKING_FORCE_H
+
+#include "btDeformableLagrangianForce.h"
+
+class btDeformableMousePickingForce : public btDeformableLagrangianForce
+{
+    // If true, the damping force will be in the direction of the spring
+    // If false, the damping force will be in the direction of the velocity
+    btScalar m_elasticStiffness, m_dampingStiffness;
+    const btSoftBody::Face& m_face;
+    btVector3 m_mouse_pos;
+    btScalar m_maxForce;
+public:
+    typedef btAlignedObjectArray<btVector3> TVStack;
+    btDeformableMousePickingForce(btScalar k, btScalar d, const btSoftBody::Face& face, btVector3 mouse_pos, btScalar maxForce = 0.3) : m_elasticStiffness(k), m_dampingStiffness(d), m_face(face), m_mouse_pos(mouse_pos),  m_maxForce(maxForce)
+    {
+    }
+    
+    virtual void addScaledForces(btScalar scale, TVStack& force)
+    {
+        addScaledDampingForce(scale, force);
+        addScaledElasticForce(scale, force);
+    }
+    
+    virtual void addScaledExplicitForce(btScalar scale, TVStack& force)
+    {
+        addScaledElasticForce(scale, force);
+    }
+    
+    virtual void addScaledDampingForce(btScalar scale, TVStack& force)
+    {
+        for (int i = 0; i < 3; ++i)
+        {
+            btVector3 v_diff = m_face.m_n[i]->m_v;
+            btVector3 scaled_force = scale * m_dampingStiffness * v_diff;
+            if ((m_face.m_n[i]->m_x - m_mouse_pos).norm() > SIMD_EPSILON)
+            {
+                btVector3 dir = (m_face.m_n[i]->m_x - m_mouse_pos).normalized();
+                scaled_force = scale * m_dampingStiffness * v_diff.dot(dir) * dir;
+            }
+            force[m_face.m_n[i]->index] -= scaled_force;
+        }
+    }
+    
+    virtual void addScaledElasticForce(btScalar scale, TVStack& force)
+    {
+        btScalar scaled_stiffness = scale * m_elasticStiffness;
+        for (int i = 0; i < 3; ++i)
+        {
+            btVector3 dir = (m_face.m_n[i]->m_q - m_mouse_pos);
+            btVector3 scaled_force = scaled_stiffness * dir;
+            if (scaled_force.safeNorm() > m_maxForce)
+            {
+                scaled_force.safeNormalize();
+                scaled_force *= m_maxForce;
+            }
+            force[m_face.m_n[i]->index] -= scaled_force;
+        }
+    }
+    
+    virtual void addScaledDampingForceDifferential(btScalar scale, const TVStack& dv, TVStack& df)
+    {
+        btScalar scaled_k_damp = m_dampingStiffness * scale;
+        for (int i = 0; i < 3; ++i)
+        {
+            btVector3 local_scaled_df = scaled_k_damp * dv[m_face.m_n[i]->index];
+            if ((m_face.m_n[i]->m_x - m_mouse_pos).norm() > SIMD_EPSILON)
+            {
+                btVector3 dir = (m_face.m_n[i]->m_x - m_mouse_pos).normalized();
+                local_scaled_df= scaled_k_damp * dv[m_face.m_n[i]->index].dot(dir) * dir;
+            }
+            df[m_face.m_n[i]->index] -= local_scaled_df;
+        }
+    }
+    
+    virtual void buildDampingForceDifferentialDiagonal(btScalar scale, TVStack& diagA){}
+    
+    virtual double totalElasticEnergy(btScalar dt)
+    {
+        double energy = 0;
+        for (int i = 0; i < 3; ++i)
+        {
+            btVector3 dir = (m_face.m_n[i]->m_q - m_mouse_pos);
+            btVector3 scaled_force = m_elasticStiffness * dir;
+            if (scaled_force.safeNorm() > m_maxForce)
+            {
+                scaled_force.safeNormalize();
+                scaled_force *= m_maxForce;
+            }
+            energy += 0.5 * scaled_force.dot(dir);
+        }
+        return energy;
+    }
+    
+    virtual double totalDampingEnergy(btScalar dt)
+    {
+        double energy = 0;
+        for (int i = 0; i < 3; ++i)
+        {
+            btVector3 v_diff = m_face.m_n[i]->m_v;
+            btVector3 scaled_force = m_dampingStiffness * v_diff;
+            if ((m_face.m_n[i]->m_x - m_mouse_pos).norm() > SIMD_EPSILON)
+            {
+                btVector3 dir = (m_face.m_n[i]->m_x - m_mouse_pos).normalized();
+                scaled_force = m_dampingStiffness * v_diff.dot(dir) * dir;
+            }
+            energy -= scaled_force.dot(m_face.m_n[i]->m_v) / dt;
+        }
+        return energy;
+    }
+    
+    virtual void addScaledElasticForceDifferential(btScalar scale, const TVStack& dx, TVStack& df)
+    {
+        //TODO
+    }
+    
+    void setMousePos(const btVector3& p)
+    {
+        m_mouse_pos = p;
+    }
+    
+    virtual btDeformableLagrangianForceType getForceType()
+    {
+        return BT_MOUSE_PICKING_FORCE;
+    }
+    
+};
+
+#endif /* btMassSpring_h */
diff --git a/thirdparty/bullet/BulletSoftBody/btDeformableMultiBodyConstraintSolver.cpp b/thirdparty/bullet/BulletSoftBody/btDeformableMultiBodyConstraintSolver.cpp
index 06f95d69f6..c8cc47923e 100644
--- a/thirdparty/bullet/BulletSoftBody/btDeformableMultiBodyConstraintSolver.cpp
+++ b/thirdparty/bullet/BulletSoftBody/btDeformableMultiBodyConstraintSolver.cpp
@@ -32,7 +32,7 @@ btScalar btDeformableMultiBodyConstraintSolver::solveDeformableGroupIterations(b
             m_leastSquaresResidual = solveSingleIteration(iteration, bodies, numBodies, manifoldPtr, numManifolds, constraints, numConstraints, infoGlobal, debugDrawer);
             // solver body velocity -> rigid body velocity
             solverBodyWriteBack(infoGlobal);
-            btScalar deformableResidual = m_deformableSolver->solveContactConstraints(deformableBodies,numDeformableBodies);
+            btScalar deformableResidual = m_deformableSolver->solveContactConstraints(deformableBodies,numDeformableBodies, infoGlobal);
             // update rigid body velocity in rigid/deformable contact
             m_leastSquaresResidual = btMax(m_leastSquaresResidual, deformableResidual);
             // solver body velocity <- rigid body velocity
@@ -112,7 +112,7 @@ void btDeformableMultiBodyConstraintSolver::solveGroupCacheFriendlySplitImpulseI
     if (infoGlobal.m_splitImpulse)
     {
         {
-            m_deformableSolver->splitImpulseSetup(infoGlobal);
+//            m_deformableSolver->splitImpulseSetup(infoGlobal);
             for (iteration = 0; iteration < infoGlobal.m_numIterations; iteration++)
             {
                 btScalar leastSquaresResidual = 0.f;
@@ -127,8 +127,8 @@ void btDeformableMultiBodyConstraintSolver::solveGroupCacheFriendlySplitImpulseI
                         leastSquaresResidual = btMax(leastSquaresResidual, residual * residual);
                     }
                     // solve the position correction between deformable and rigid/multibody
-                    btScalar residual = m_deformableSolver->solveSplitImpulse(infoGlobal);
-                    leastSquaresResidual = btMax(leastSquaresResidual, residual * residual);
+//                    btScalar residual = m_deformableSolver->solveSplitImpulse(infoGlobal);
+//                    leastSquaresResidual = btMax(leastSquaresResidual, residual * residual);
                 }
                 if (leastSquaresResidual <= infoGlobal.m_leastSquaresResidualThreshold || iteration >= (infoGlobal.m_numIterations - 1))
                 {
diff --git a/thirdparty/bullet/BulletSoftBody/btDeformableMultiBodyDynamicsWorld.cpp b/thirdparty/bullet/BulletSoftBody/btDeformableMultiBodyDynamicsWorld.cpp
index 618e5c0d7b..6b742978ef 100644
--- a/thirdparty/bullet/BulletSoftBody/btDeformableMultiBodyDynamicsWorld.cpp
+++ b/thirdparty/bullet/BulletSoftBody/btDeformableMultiBodyDynamicsWorld.cpp
@@ -22,7 +22,6 @@ Call internalStepSimulation multiple times, to achieve 240Hz (4 steps of 60Hz).
 2. Detect discrete collisions between rigid and deformable bodies at position x_{n+1}^* = x_n + dt * v_{n+1}^*.
 
 3a. Solve all constraints, including LCP. Contact, position correction due to numerical drift, friction, and anchors for deformable.
-    TODO: add option for positional drift correction (using vel_target += erp * pos_error/dt
 
 3b. 5 Newton steps (multiple step). Conjugent Gradient solves linear system. Deformable Damping: Then velocities of deformable bodies v_{n+1} are solved in
         M(v_{n+1} - v_{n+1}^*) = damping_force * dt / mass,
@@ -58,14 +57,20 @@ m_deformableBodySolver(deformableBodySolver), m_solverCallback(0)
 	m_sbi.water_density = 0;
 	m_sbi.water_offset = 0;
 	m_sbi.water_normal = btVector3(0, 0, 0);
-	m_sbi.m_gravity.setValue(0, -10, 0);
+	m_sbi.m_gravity.setValue(0, -9.8, 0);
 	m_internalTime = 0.0;
 	m_implicit = false;
 	m_lineSearch = false;
-	m_selfCollision = true;
+    m_useProjection = true;
+	m_ccdIterations = 5;
 	m_solverDeformableBodyIslandCallback = new DeformableBodyInplaceSolverIslandCallback(constraintSolver, dispatcher);
 }
 
+btDeformableMultiBodyDynamicsWorld::~btDeformableMultiBodyDynamicsWorld()
+{
+    delete m_solverDeformableBodyIslandCallback;
+}
+
 void btDeformableMultiBodyDynamicsWorld::internalSingleStepSimulation(btScalar timeStep)
 {
     BT_PROFILE("internalSingleStepSimulation");
@@ -74,20 +79,16 @@ void btDeformableMultiBodyDynamicsWorld::internalSingleStepSimulation(btScalar t
         (*m_internalPreTickCallback)(this, timeStep);
     }
     reinitialize(timeStep);
+    
     // add gravity to velocity of rigid and multi bodys
     applyRigidBodyGravity(timeStep);
     
     ///apply gravity and explicit force to velocity, predict motion
     predictUnconstraintMotion(timeStep);
     
-    ///perform collision detection
+    ///perform collision detection that involves rigid/multi bodies
     btMultiBodyDynamicsWorld::performDiscreteCollisionDetection();
     
-    if (m_selfCollision)
-    {
-        softBodySelfCollision();
-    }
-    
     btMultiBodyDynamicsWorld::calculateSimulationIslands();
     
     beforeSolverCallbacks(timeStep);
@@ -96,7 +97,13 @@ void btDeformableMultiBodyDynamicsWorld::internalSingleStepSimulation(btScalar t
     solveConstraints(timeStep);
     
     afterSolverCallbacks(timeStep);
-    
+
+	performDeformableCollisionDetection();
+
+    applyRepulsionForce(timeStep);
+
+    performGeometricCollisions(timeStep);
+
     integrateTransforms(timeStep);
     
     ///update vehicle simulation
@@ -107,6 +114,27 @@ void btDeformableMultiBodyDynamicsWorld::internalSingleStepSimulation(btScalar t
     // ///////////////////////////////
 }
 
+void btDeformableMultiBodyDynamicsWorld::performDeformableCollisionDetection()
+{
+	for (int i = 0; i < m_softBodies.size(); ++i)
+	{
+		m_softBodies[i]->m_softSoftCollision = true;
+	}
+	
+	for (int i = 0; i < m_softBodies.size(); ++i)
+	{
+		for (int j = i; j < m_softBodies.size(); ++j)
+		{
+			m_softBodies[i]->defaultCollisionHandler(m_softBodies[j]);
+		}
+	}
+	
+	for (int i = 0; i < m_softBodies.size(); ++i)
+	{
+		m_softBodies[i]->m_softSoftCollision = false;
+	}
+}
+
 void btDeformableMultiBodyDynamicsWorld::updateActivationState(btScalar timeStep)
 {
     for (int i = 0; i < m_softBodies.size(); i++)
@@ -131,10 +159,106 @@ void btDeformableMultiBodyDynamicsWorld::updateActivationState(btScalar timeStep
     btMultiBodyDynamicsWorld::updateActivationState(timeStep);
 }
 
+void btDeformableMultiBodyDynamicsWorld::applyRepulsionForce(btScalar timeStep)
+{
+    BT_PROFILE("btDeformableMultiBodyDynamicsWorld::applyRepulsionForce");
+    for (int i = 0; i < m_softBodies.size(); i++)
+    {
+        btSoftBody* psb = m_softBodies[i];
+        if (psb->isActive())
+        {
+			psb->applyRepulsionForce(timeStep, true);
+        }
+    }
+}
+
+void btDeformableMultiBodyDynamicsWorld::performGeometricCollisions(btScalar timeStep)
+{
+	BT_PROFILE("btDeformableMultiBodyDynamicsWorld::performGeometricCollisions");
+    // refit the BVH tree for CCD
+	for (int i = 0; i < m_softBodies.size(); ++i)
+	{
+		btSoftBody* psb = m_softBodies[i];
+		if (psb->isActive())
+		{
+			m_softBodies[i]->updateFaceTree(true, false);
+			m_softBodies[i]->updateNodeTree(true, false);
+			for (int j = 0; j < m_softBodies[i]->m_faces.size(); ++j)
+			{
+				btSoftBody::Face& f = m_softBodies[i]->m_faces[j];
+				f.m_n0 = (f.m_n[1]->m_x - f.m_n[0]->m_x).cross(f.m_n[2]->m_x - f.m_n[0]->m_x);
+			}
+		}
+	}
+
+	// clear contact points & update DBVT
+	for (int r = 0; r < m_ccdIterations; ++r)
+	{
+		for (int i = 0; i < m_softBodies.size(); ++i)
+		{
+			btSoftBody* psb = m_softBodies[i];
+			if (psb->isActive())
+			{
+				// clear contact points in the previous iteration
+				psb->m_faceNodeContacts.clear();
+
+				// update m_q and normals for CCD calculation
+				for (int j = 0; j < psb->m_nodes.size(); ++j)
+				{
+					psb->m_nodes[j].m_q = psb->m_nodes[j].m_x + timeStep * psb->m_nodes[j].m_v;
+				}
+				for (int j = 0; j < psb->m_faces.size(); ++j)
+				{
+					btSoftBody::Face& f = psb->m_faces[j];
+					f.m_n1 = (f.m_n[1]->m_q - f.m_n[0]->m_q).cross(f.m_n[2]->m_q - f.m_n[0]->m_q);
+					f.m_vn = (f.m_n[1]->m_v - f.m_n[0]->m_v).cross(f.m_n[2]->m_v - f.m_n[0]->m_v) * timeStep * timeStep;
+				}
+			}
+        }
+
+		// apply CCD to register new contact points
+		for (int i = 0; i < m_softBodies.size(); ++i)
+		{
+			for (int j = i; j < m_softBodies.size(); ++j)
+			{
+				btSoftBody* psb1 = m_softBodies[i];
+				btSoftBody* psb2 = m_softBodies[j];
+				if (psb1->isActive() && psb2->isActive())
+				{
+					m_softBodies[i]->geometricCollisionHandler(m_softBodies[j]);
+				}
+			}
+        }
+
+		int penetration_count = 0;
+		for (int i = 0; i < m_softBodies.size(); ++i)
+		{
+			btSoftBody* psb = m_softBodies[i];
+			if (psb->isActive())
+			{
+				penetration_count += psb->m_faceNodeContacts.size();
+			}
+		}
+		if (penetration_count == 0)
+		{
+			break;
+		}
+
+		// apply inelastic impulse
+		for (int i = 0; i < m_softBodies.size(); ++i)
+		{
+			btSoftBody* psb = m_softBodies[i];
+			if (psb->isActive())
+			{
+				psb->applyRepulsionForce(timeStep, false);
+			}
+		}
+	}
+}
 
 void btDeformableMultiBodyDynamicsWorld::softBodySelfCollision()
 {
-    m_deformableBodySolver->updateSoftBodies();
+    BT_PROFILE("btDeformableMultiBodyDynamicsWorld::softBodySelfCollision");
     for (int i = 0; i < m_softBodies.size(); i++)
     {
         btSoftBody* psb = m_softBodies[i];
@@ -192,8 +316,6 @@ void btDeformableMultiBodyDynamicsWorld::integrateTransforms(btScalar timeStep)
                 }
             }
             node.m_x  =  node.m_x + timeStep * node.m_v;
-            node.m_v -= node.m_vsplit;
-            node.m_vsplit.setZero();
             node.m_q = node.m_x;
             node.m_vn = node.m_v;
         }
@@ -255,6 +377,7 @@ void btDeformableMultiBodyDynamicsWorld::integrateTransforms(btScalar timeStep)
 
 void btDeformableMultiBodyDynamicsWorld::solveConstraints(btScalar timeStep)
 {
+    BT_PROFILE("btDeformableMultiBodyDynamicsWorld::solveConstraints");
     // save v_{n+1}^* velocity after explicit forces
     m_deformableBodySolver->backupVelocity();
     
@@ -265,8 +388,11 @@ void btDeformableMultiBodyDynamicsWorld::solveConstraints(btScalar timeStep)
     solveContactConstraints();
     
     // set up the directions in which the velocity does not change in the momentum solve
-    m_deformableBodySolver->m_objective->m_projection.setProjection();
-    
+    if (m_useProjection)
+        m_deformableBodySolver->m_objective->m_projection.setProjection();
+    else
+        m_deformableBodySolver->m_objective->m_projection.setLagrangeMultiplier();
+
     // for explicit scheme, m_backupVelocity = v_{n+1}^*
     // for implicit scheme, m_backupVelocity = v_n
     // Here, set dv = v_{n+1} - v_n for nodes in contact
@@ -280,7 +406,7 @@ void btDeformableMultiBodyDynamicsWorld::solveConstraints(btScalar timeStep)
 void btDeformableMultiBodyDynamicsWorld::setupConstraints()
 {
     // set up constraints between multibody and deformable bodies
-    m_deformableBodySolver->setConstraints();
+    m_deformableBodySolver->setConstraints(m_solverInfo);
     
     // set up constraints among multibodies
     {
@@ -403,6 +529,17 @@ void btDeformableMultiBodyDynamicsWorld::reinitialize(btScalar timeStep)
     dispatchInfo.m_stepCount = 0;
     dispatchInfo.m_debugDraw = btMultiBodyDynamicsWorld::getDebugDrawer();
     btMultiBodyDynamicsWorld::getSolverInfo().m_timeStep = timeStep;
+    if (m_useProjection)
+    {
+        m_deformableBodySolver->m_useProjection = true;
+//        m_deformableBodySolver->m_objective->m_projection.m_useStrainLimiting = true;
+        m_deformableBodySolver->m_objective->m_preconditioner =  m_deformableBodySolver->m_objective->m_massPreconditioner;
+    }
+    else
+    {
+        m_deformableBodySolver->m_objective->m_preconditioner =  m_deformableBodySolver->m_objective->m_KKTPreconditioner;
+    }
+        
 }
 
 
@@ -566,6 +703,24 @@ void btDeformableMultiBodyDynamicsWorld::addForce(btSoftBody* psb, btDeformableL
     }
 }
 
+void btDeformableMultiBodyDynamicsWorld::removeForce(btSoftBody* psb, btDeformableLagrangianForce* force)
+{
+    btAlignedObjectArray<btDeformableLagrangianForce*>& forces = m_deformableBodySolver->m_objective->m_lf;
+    int removed_index = -1;
+    for (int i = 0; i < forces.size(); ++i)
+    {
+        if (forces[i]->getForceType() == force->getForceType())
+        {
+            forces[i]->removeSoftBody(psb);
+            if (forces[i]->m_softBodies.size() == 0)
+                removed_index = i;
+            break;
+        }
+    }
+    if (removed_index >= 0)
+        forces.removeAtIndex(removed_index);
+}
+
 void btDeformableMultiBodyDynamicsWorld::removeSoftBody(btSoftBody* body)
 {
     m_softBodies.remove(body);
diff --git a/thirdparty/bullet/BulletSoftBody/btDeformableMultiBodyDynamicsWorld.h b/thirdparty/bullet/BulletSoftBody/btDeformableMultiBodyDynamicsWorld.h
index 7630385767..76b58a0378 100644
--- a/thirdparty/bullet/BulletSoftBody/btDeformableMultiBodyDynamicsWorld.h
+++ b/thirdparty/bullet/BulletSoftBody/btDeformableMultiBodyDynamicsWorld.h
@@ -46,10 +46,10 @@ class btDeformableMultiBodyDynamicsWorld : public btMultiBodyDynamicsWorld
     bool m_drawClusterTree;
     btSoftBodyWorldInfo m_sbi;
     btScalar m_internalTime;
-    int m_contact_iterations;
+    int m_ccdIterations;
     bool m_implicit;
     bool m_lineSearch;
-    bool m_selfCollision;
+    bool m_useProjection;
     DeformableBodyInplaceSolverIslandCallback* m_solverDeformableBodyIslandCallback;
     
     typedef void (*btSolverCallback)(btScalar time, btDeformableMultiBodyDynamicsWorld* world);
@@ -80,9 +80,7 @@ public:
         m_solverCallback = cb;
     }
     
-    virtual ~btDeformableMultiBodyDynamicsWorld()
-    {
-    }
+    virtual ~btDeformableMultiBodyDynamicsWorld();
     
     virtual btMultiBodyDynamicsWorld* getMultiBodyDynamicsWorld()
     {
@@ -133,6 +131,8 @@ public:
     
     void addForce(btSoftBody* psb, btDeformableLagrangianForce* force);
     
+    void removeForce(btSoftBody* psb, btDeformableLagrangianForce* force);
+    
     void removeSoftBody(btSoftBody* body);
     
     void removeCollisionObject(btCollisionObject* collisionObject);
@@ -142,6 +142,8 @@ public:
     
     void setupConstraints();
     
+    void performDeformableCollisionDetection();
+    
     void solveMultiBodyConstraints();
     
     void solveContactConstraints();
@@ -159,7 +161,151 @@ public:
     {
         m_lineSearch = lineSearch;
     }
+    
+    void applyRepulsionForce(btScalar timeStep);
+    
+    void performGeometricCollisions(btScalar timeStep);
+    
+    struct btDeformableSingleRayCallback : public btBroadphaseRayCallback
+    {
+        btVector3 m_rayFromWorld;
+        btVector3 m_rayToWorld;
+        btTransform m_rayFromTrans;
+        btTransform m_rayToTrans;
+        btVector3 m_hitNormal;
+        
+        const btDeformableMultiBodyDynamicsWorld* m_world;
+        btCollisionWorld::RayResultCallback& m_resultCallback;
+        
+        btDeformableSingleRayCallback(const btVector3& rayFromWorld, const btVector3& rayToWorld, const btDeformableMultiBodyDynamicsWorld* world, btCollisionWorld::RayResultCallback& resultCallback)
+        : m_rayFromWorld(rayFromWorld),
+        m_rayToWorld(rayToWorld),
+        m_world(world),
+        m_resultCallback(resultCallback)
+        {
+            m_rayFromTrans.setIdentity();
+            m_rayFromTrans.setOrigin(m_rayFromWorld);
+            m_rayToTrans.setIdentity();
+            m_rayToTrans.setOrigin(m_rayToWorld);
+            
+            btVector3 rayDir = (rayToWorld - rayFromWorld);
+            
+            rayDir.normalize();
+            ///what about division by zero? --> just set rayDirection[i] to INF/1e30
+            m_rayDirectionInverse[0] = rayDir[0] == btScalar(0.0) ? btScalar(1e30) : btScalar(1.0) / rayDir[0];
+            m_rayDirectionInverse[1] = rayDir[1] == btScalar(0.0) ? btScalar(1e30) : btScalar(1.0) / rayDir[1];
+            m_rayDirectionInverse[2] = rayDir[2] == btScalar(0.0) ? btScalar(1e30) : btScalar(1.0) / rayDir[2];
+            m_signs[0] = m_rayDirectionInverse[0] < 0.0;
+            m_signs[1] = m_rayDirectionInverse[1] < 0.0;
+            m_signs[2] = m_rayDirectionInverse[2] < 0.0;
+            
+            m_lambda_max = rayDir.dot(m_rayToWorld - m_rayFromWorld);
+        }
+        
+        virtual bool process(const btBroadphaseProxy* proxy)
+        {
+            ///terminate further ray tests, once the closestHitFraction reached zero
+            if (m_resultCallback.m_closestHitFraction == btScalar(0.f))
+                return false;
+            
+            btCollisionObject* collisionObject = (btCollisionObject*)proxy->m_clientObject;
+            
+            //only perform raycast if filterMask matches
+            if (m_resultCallback.needsCollision(collisionObject->getBroadphaseHandle()))
+            {
+                //RigidcollisionObject* collisionObject = ctrl->GetRigidcollisionObject();
+                //btVector3 collisionObjectAabbMin,collisionObjectAabbMax;
+#if 0
+#ifdef RECALCULATE_AABB
+                btVector3 collisionObjectAabbMin,collisionObjectAabbMax;
+                collisionObject->getCollisionShape()->getAabb(collisionObject->getWorldTransform(),collisionObjectAabbMin,collisionObjectAabbMax);
+#else
+                //getBroadphase()->getAabb(collisionObject->getBroadphaseHandle(),collisionObjectAabbMin,collisionObjectAabbMax);
+                const btVector3& collisionObjectAabbMin = collisionObject->getBroadphaseHandle()->m_aabbMin;
+                const btVector3& collisionObjectAabbMax = collisionObject->getBroadphaseHandle()->m_aabbMax;
+#endif
+#endif
+                //btScalar hitLambda = m_resultCallback.m_closestHitFraction;
+                //culling already done by broadphase
+                //if (btRayAabb(m_rayFromWorld,m_rayToWorld,collisionObjectAabbMin,collisionObjectAabbMax,hitLambda,m_hitNormal))
+                {
+                    m_world->rayTestSingle(m_rayFromTrans, m_rayToTrans,
+                                           collisionObject,
+                                           collisionObject->getCollisionShape(),
+                                           collisionObject->getWorldTransform(),
+                                           m_resultCallback);
+                }
+            }
+            return true;
+        }
+    };
 
+    
+    
+    void rayTest(const btVector3& rayFromWorld, const btVector3& rayToWorld, RayResultCallback& resultCallback) const
+    {
+        BT_PROFILE("rayTest");
+        /// use the broadphase to accelerate the search for objects, based on their aabb
+        /// and for each object with ray-aabb overlap, perform an exact ray test
+        btDeformableSingleRayCallback rayCB(rayFromWorld, rayToWorld, this, resultCallback);
+        
+#ifndef USE_BRUTEFORCE_RAYBROADPHASE
+        m_broadphasePairCache->rayTest(rayFromWorld, rayToWorld, rayCB);
+#else
+        for (int i = 0; i < this->getNumCollisionObjects(); i++)
+        {
+            rayCB.process(m_collisionObjects[i]->getBroadphaseHandle());
+        }
+#endif  //USE_BRUTEFORCE_RAYBROADPHASE
+    }
+    
+    void rayTestSingle(const btTransform& rayFromTrans, const btTransform& rayToTrans,
+                                                     btCollisionObject* collisionObject,
+                                                     const btCollisionShape* collisionShape,
+                                                     const btTransform& colObjWorldTransform,
+                                                     RayResultCallback& resultCallback) const
+    {
+        if (collisionShape->isSoftBody())
+        {
+            btSoftBody* softBody = btSoftBody::upcast(collisionObject);
+            if (softBody)
+            {
+                btSoftBody::sRayCast softResult;
+                if (softBody->rayFaceTest(rayFromTrans.getOrigin(), rayToTrans.getOrigin(), softResult))
+                {
+                    if (softResult.fraction <= resultCallback.m_closestHitFraction)
+                    {
+                        btCollisionWorld::LocalShapeInfo shapeInfo;
+                        shapeInfo.m_shapePart = 0;
+                        shapeInfo.m_triangleIndex = softResult.index;
+                        // get the normal
+                        btVector3 rayDir = rayToTrans.getOrigin() - rayFromTrans.getOrigin();
+                        btVector3 normal = -rayDir;
+                        normal.normalize();
+                        {
+                            normal = softBody->m_faces[softResult.index].m_normal;
+                            if (normal.dot(rayDir) > 0)
+                            {
+                                // normal always point toward origin of the ray
+                                normal = -normal;
+                            }
+                        }
+                        
+                        btCollisionWorld::LocalRayResult rayResult(collisionObject,
+                                                                   &shapeInfo,
+                                                                   normal,
+                                                                   softResult.fraction);
+                        bool normalInWorldSpace = true;
+                        resultCallback.addSingleResult(rayResult, normalInWorldSpace);
+                    }
+                }
+            }
+        }
+        else
+        {
+            btCollisionWorld::rayTestSingle(rayFromTrans, rayToTrans, collisionObject, collisionShape, colObjWorldTransform, resultCallback);
+        }
+    }
 };
 
 #endif  //BT_DEFORMABLE_MULTIBODY_DYNAMICS_WORLD_H
diff --git a/thirdparty/bullet/BulletSoftBody/btDeformableNeoHookeanForce.h b/thirdparty/bullet/BulletSoftBody/btDeformableNeoHookeanForce.h
index 3d06e304d2..d89bc4aca4 100644
--- a/thirdparty/bullet/BulletSoftBody/btDeformableNeoHookeanForce.h
+++ b/thirdparty/bullet/BulletSoftBody/btDeformableNeoHookeanForce.h
@@ -24,21 +24,65 @@ class btDeformableNeoHookeanForce : public btDeformableLagrangianForce
 {
 public:
     typedef btAlignedObjectArray<btVector3> TVStack;
-    btScalar m_mu, m_lambda;
+	btScalar m_mu, m_lambda; // Lame Parameters
+	btScalar m_E, m_nu;  // Young's modulus and Poisson ratio
     btScalar m_mu_damp, m_lambda_damp;
     btDeformableNeoHookeanForce(): m_mu(1), m_lambda(1)
     {
         btScalar damping = 0.05;
         m_mu_damp = damping * m_mu;
         m_lambda_damp = damping * m_lambda;
+		updateYoungsModulusAndPoissonRatio();
     }
     
     btDeformableNeoHookeanForce(btScalar mu, btScalar lambda, btScalar damping = 0.05): m_mu(mu), m_lambda(lambda)
     {
         m_mu_damp = damping * m_mu;
         m_lambda_damp = damping * m_lambda;
+		updateYoungsModulusAndPoissonRatio();
     }
-    
+
+	void updateYoungsModulusAndPoissonRatio()
+	{
+		// conversion from Lame Parameters to Young's modulus and Poisson ratio
+		// https://en.wikipedia.org/wiki/Lam%C3%A9_parameters
+		m_E  = m_mu * (3*m_lambda + 2*m_mu)/(m_lambda + m_mu);
+		m_nu = m_lambda * 0.5 / (m_mu + m_lambda);
+	}
+
+	void updateLameParameters()
+	{
+		// conversion from Young's modulus and Poisson ratio to Lame Parameters
+		// https://en.wikipedia.org/wiki/Lam%C3%A9_parameters
+		m_mu = m_E * 0.5 / (1 + m_nu);
+		m_lambda = m_E * m_nu / ((1 + m_nu) * (1- 2*m_nu));
+	}
+
+    void setYoungsModulus(btScalar E)
+    {
+		m_E = E;
+		updateLameParameters();
+    }
+
+	void setPoissonRatio(btScalar nu)
+	{
+		m_nu = nu;
+		updateLameParameters();
+	}
+	
+	void setDamping(btScalar damping)
+	{
+		m_mu_damp = damping * m_mu;
+		m_lambda_damp = damping * m_lambda;
+	}
+
+	void setLameParameters(btScalar mu, btScalar lambda)
+	{
+		m_mu = mu;
+		m_lambda = lambda;
+		updateYoungsModulusAndPoissonRatio();
+	}
+
     virtual void addScaledForces(btScalar scale, TVStack& force)
     {
         addScaledDampingForce(scale, force);
@@ -269,6 +313,8 @@ public:
         }
     }
     
+    virtual void buildDampingForceDifferentialDiagonal(btScalar scale, TVStack& diagA){}
+    
     virtual void addScaledElasticForceDifferential(btScalar scale, const TVStack& dx, TVStack& df)
     {
         int numNodes = getNumNodes();
diff --git a/thirdparty/bullet/BulletSoftBody/btPreconditioner.h b/thirdparty/bullet/BulletSoftBody/btPreconditioner.h
index d712420381..c2db448ef8 100644
--- a/thirdparty/bullet/BulletSoftBody/btPreconditioner.h
+++ b/thirdparty/bullet/BulletSoftBody/btPreconditioner.h
@@ -68,12 +68,221 @@ public:
     virtual void operator()(const TVStack& x, TVStack& b)
     {
         btAssert(b.size() == x.size());
-        btAssert(m_inv_mass.size() == x.size());
-        for (int i = 0; i < b.size(); ++i)
+        btAssert(m_inv_mass.size() <= x.size());
+        for (int i = 0; i < m_inv_mass.size(); ++i)
         {
             b[i] = x[i] * m_inv_mass[i];
         }
+        for (int i = m_inv_mass.size(); i < b.size(); ++i)
+        {
+            b[i] = x[i];
+        }
+    }
+};
+
+
+class KKTPreconditioner : public Preconditioner
+{
+    const btAlignedObjectArray<btSoftBody *>& m_softBodies;
+    const btDeformableContactProjection& m_projections;
+    const btAlignedObjectArray<btDeformableLagrangianForce*>& m_lf;
+    TVStack m_inv_A, m_inv_S;
+    const btScalar& m_dt;
+    const bool& m_implicit;
+public:
+    KKTPreconditioner(const btAlignedObjectArray<btSoftBody *>& softBodies, const btDeformableContactProjection& projections, const btAlignedObjectArray<btDeformableLagrangianForce*>& lf, const btScalar& dt, const bool& implicit)
+    : m_softBodies(softBodies)
+    , m_projections(projections)
+    , m_lf(lf)
+    , m_dt(dt)
+    , m_implicit(implicit)
+    {
+    }
+    
+    virtual void reinitialize(bool nodeUpdated)
+    {
+        if (nodeUpdated)
+        {
+            int num_nodes = 0;
+            for (int i = 0; i < m_softBodies.size(); ++i)
+            {
+                btSoftBody* psb = m_softBodies[i];
+                num_nodes += psb->m_nodes.size();
+            }
+            m_inv_A.resize(num_nodes);
+        }
+        buildDiagonalA(m_inv_A);
+        for (int i = 0; i < m_inv_A.size(); ++i)
+        {
+//            printf("A[%d] = %f, %f, %f \n", i, m_inv_A[i][0], m_inv_A[i][1], m_inv_A[i][2]);
+            for (int d = 0; d < 3; ++d)
+            {
+                m_inv_A[i][d] = (m_inv_A[i][d] == 0) ? 0.0 : 1.0/ m_inv_A[i][d];
+            }
+        }
+        m_inv_S.resize(m_projections.m_lagrangeMultipliers.size());
+//        printf("S.size() = %d \n", m_inv_S.size());
+        buildDiagonalS(m_inv_A, m_inv_S);
+        for (int i = 0; i < m_inv_S.size(); ++i)
+        {
+//            printf("S[%d] = %f, %f, %f \n", i, m_inv_S[i][0], m_inv_S[i][1], m_inv_S[i][2]);
+            for (int d = 0; d < 3; ++d)
+            {
+                m_inv_S[i][d] = (m_inv_S[i][d] == 0) ? 0.0 : 1.0/ m_inv_S[i][d];
+            }
+        }
+    }
+    
+    void buildDiagonalA(TVStack& diagA) const
+    {
+        size_t counter = 0;
+        for (int i = 0; i < m_softBodies.size(); ++i)
+        {
+            btSoftBody* psb = m_softBodies[i];
+            for (int j = 0; j < psb->m_nodes.size(); ++j)
+            {
+                const btSoftBody::Node& node = psb->m_nodes[j];
+                diagA[counter] = (node.m_im == 0) ? btVector3(0,0,0) : btVector3(1.0/node.m_im, 1.0 / node.m_im, 1.0 / node.m_im);
+                ++counter;
+            }
+        }
+        if (m_implicit)
+        {
+            printf("implicit not implemented\n");
+            btAssert(false);
+        }
+        for (int i = 0; i < m_lf.size(); ++i)
+        {
+            // add damping matrix
+            m_lf[i]->buildDampingForceDifferentialDiagonal(-m_dt, diagA);
+        }
+    }
+    
+    void buildDiagonalS(const TVStack& inv_A, TVStack& diagS)
+    {
+        for (int c = 0; c < m_projections.m_lagrangeMultipliers.size(); ++c)
+        {
+            // S[k,k] = e_k^T * C A_d^-1 C^T * e_k
+            const LagrangeMultiplier& lm = m_projections.m_lagrangeMultipliers[c];
+            btVector3& t = diagS[c];
+            t.setZero();
+            for (int j = 0; j < lm.m_num_constraints; ++j)
+            {
+                for (int i = 0; i < lm.m_num_nodes; ++i)
+                {
+                    for (int d = 0; d < 3; ++d)
+                    {
+                        t[j] += inv_A[lm.m_indices[i]][d] * lm.m_dirs[j][d] * lm.m_dirs[j][d] * lm.m_weights[i] * lm.m_weights[i];
+                    }
+                }
+            }
+        }
+    }
+#define USE_FULL_PRECONDITIONER
+#ifndef USE_FULL_PRECONDITIONER
+    virtual void operator()(const TVStack& x, TVStack& b)
+    {
+        btAssert(b.size() == x.size());
+        for (int i = 0; i < m_inv_A.size(); ++i)
+        {
+            b[i] = x[i] * m_inv_A[i];
+        }
+        int offset = m_inv_A.size();
+        for (int i = 0; i < m_inv_S.size(); ++i)
+        {
+            b[i+offset] = x[i+offset] * m_inv_S[i];
+        }
+    }
+#else
+    virtual void operator()(const TVStack& x, TVStack& b)
+    {
+        btAssert(b.size() == x.size());
+        int offset = m_inv_A.size();
+
+        for (int i = 0; i < m_inv_A.size(); ++i)
+        {
+            b[i] = x[i] * m_inv_A[i];
+        }
+
+        for (int i = 0; i < m_inv_S.size(); ++i)
+        {
+            b[i+offset].setZero();
+        }
+
+        for (int c = 0; c < m_projections.m_lagrangeMultipliers.size(); ++c)
+        {
+            const LagrangeMultiplier& lm = m_projections.m_lagrangeMultipliers[c];
+            // C * x
+            for (int d = 0; d < lm.m_num_constraints; ++d)
+            {
+                for (int i = 0; i < lm.m_num_nodes; ++i)
+                {
+                    b[offset+c][d] += lm.m_weights[i] * b[lm.m_indices[i]].dot(lm.m_dirs[d]);
+                }
+            }
+        }
+
+        for (int i = 0; i < m_inv_S.size(); ++i)
+        {
+            b[i+offset] = b[i+offset] * m_inv_S[i];
+        }
+
+        for (int i = 0; i < m_inv_A.size(); ++i)
+        {
+            b[i].setZero();
+        }
+
+        for (int c = 0; c < m_projections.m_lagrangeMultipliers.size(); ++c)
+        {
+            // C^T * lambda
+            const LagrangeMultiplier& lm = m_projections.m_lagrangeMultipliers[c];
+            for (int i = 0; i < lm.m_num_nodes; ++i)
+            {
+                for (int j = 0; j < lm.m_num_constraints; ++j)
+                {
+                    b[lm.m_indices[i]] += b[offset+c][j] * lm.m_weights[i] * lm.m_dirs[j];
+                }
+            }
+        }
+
+        for (int i = 0; i < m_inv_A.size(); ++i)
+        {
+            b[i] = (x[i] - b[i]) * m_inv_A[i];
+        }
+
+        TVStack t;
+        t.resize(b.size());
+        for (int i = 0; i < m_inv_S.size(); ++i)
+        {
+            t[i+offset] = x[i+offset] * m_inv_S[i];
+        }
+        for (int i = 0; i < m_inv_A.size(); ++i)
+        {
+            t[i].setZero();
+        }
+        for (int c = 0; c < m_projections.m_lagrangeMultipliers.size(); ++c)
+        {
+            // C^T * lambda
+            const LagrangeMultiplier& lm = m_projections.m_lagrangeMultipliers[c];
+            for (int i = 0; i < lm.m_num_nodes; ++i)
+            {
+                for (int j = 0; j < lm.m_num_constraints; ++j)
+                {
+                    t[lm.m_indices[i]] += t[offset+c][j] * lm.m_weights[i] * lm.m_dirs[j];
+                }
+            }
+        }
+        for (int i = 0; i < m_inv_A.size(); ++i)
+        {
+            b[i] += t[i] * m_inv_A[i];
+        }
+
+        for (int i = 0; i < m_inv_S.size(); ++i)
+        {
+            b[i+offset] -= x[i+offset] * m_inv_S[i];
+        }
     }
+#endif
 };
 
 #endif /* BT_PRECONDITIONER_H */
diff --git a/thirdparty/bullet/BulletSoftBody/btSoftBody.cpp b/thirdparty/bullet/BulletSoftBody/btSoftBody.cpp
index 2a458b1d80..81b846d7f8 100644
--- a/thirdparty/bullet/BulletSoftBody/btSoftBody.cpp
+++ b/thirdparty/bullet/BulletSoftBody/btSoftBody.cpp
@@ -18,6 +18,7 @@ subject to the following restrictions:
 #include "BulletSoftBody/btSoftBodySolvers.h"
 #include "btSoftBodyData.h"
 #include "LinearMath/btSerializer.h"
+#include "LinearMath/btImplicitQRSVD.h"
 #include "LinearMath/btAlignedAllocator.h"
 #include "BulletDynamics/Featherstone/btMultiBodyLinkCollider.h"
 #include "BulletDynamics/Featherstone/btMultiBodyConstraint.h"
@@ -25,6 +26,107 @@ subject to the following restrictions:
 #include "BulletCollision/CollisionShapes/btTriangleShape.h"
 #include <iostream>
 //
+static inline btDbvtNode* buildTreeBottomUp(btAlignedObjectArray<btDbvtNode*>& leafNodes, btAlignedObjectArray<btAlignedObjectArray<int> >& adj)
+{
+	int N = leafNodes.size();
+	if (N == 0)
+	{
+		return NULL;
+	}
+	while (N > 1)
+	{
+		btAlignedObjectArray<bool> marked;
+		btAlignedObjectArray<btDbvtNode*> newLeafNodes;
+		btAlignedObjectArray<std::pair<int,int> > childIds;
+		btAlignedObjectArray<btAlignedObjectArray<int> > newAdj;
+		marked.resize(N);
+		for (int i = 0; i < N; ++i)
+			marked[i] = false;
+		
+		// pair adjacent nodes into new(parent) node
+		for (int i = 0; i < N; ++i)
+		{
+			if (marked[i])
+				continue;
+			bool merged = false;
+			for (int j = 0; j < adj[i].size(); ++j)
+			{
+				int n = adj[i][j];
+				if (!marked[adj[i][j]])
+				{
+					btDbvtNode* node = new (btAlignedAlloc(sizeof(btDbvtNode), 16)) btDbvtNode();
+					node->parent = NULL;
+					node->childs[0] = leafNodes[i];
+					node->childs[1] = leafNodes[n];
+					leafNodes[i]->parent = node;
+					leafNodes[n]->parent = node;
+					newLeafNodes.push_back(node);
+					childIds.push_back(std::make_pair(i,n));
+					merged = true;
+					marked[n] = true;
+					break;
+				}
+			}
+			if (!merged)
+			{
+				newLeafNodes.push_back(leafNodes[i]);
+				childIds.push_back(std::make_pair(i,-1));
+			}
+			marked[i] = true;
+		}
+		// update adjacency matrix
+		newAdj.resize(newLeafNodes.size());
+		for (int i = 0; i < newLeafNodes.size(); ++i)
+		{
+			for (int j = i+1; j < newLeafNodes.size(); ++j)
+			{
+				bool neighbor = false;
+				const btAlignedObjectArray<int>& leftChildNeighbors = adj[childIds[i].first];
+				for (int k = 0; k < leftChildNeighbors.size(); ++k)
+				{
+					if (leftChildNeighbors[k] == childIds[j].first || leftChildNeighbors[k] == childIds[j].second)
+					{
+						neighbor = true;
+						break;
+					}
+				}
+				if (!neighbor && childIds[i].second != -1)
+				{
+					const btAlignedObjectArray<int>& rightChildNeighbors = adj[childIds[i].second];
+					for (int k = 0; k < rightChildNeighbors.size(); ++k)
+					{
+						if (rightChildNeighbors[k] == childIds[j].first || rightChildNeighbors[k] == childIds[j].second)
+						{
+							neighbor = true;
+							break;
+						}
+					}
+				}
+				if (neighbor)
+				{
+					newAdj[i].push_back(j);
+					newAdj[j].push_back(i);
+				}
+			}
+		}
+		leafNodes = newLeafNodes;
+		//this assignment leaks memory, the assignment doesn't do a deep copy, for now a manual copy
+		//adj = newAdj;
+		adj.clear();
+		adj.resize(newAdj.size());
+		for (int i = 0; i < newAdj.size(); i++)
+		{
+			for (int j = 0; j < newAdj[i].size(); j++)
+			{
+				adj[i].push_back(newAdj[i][j]);
+			}
+		}
+		N = leafNodes.size();
+	}
+	return leafNodes[0];
+}
+
+//
 btSoftBody::btSoftBody(btSoftBodyWorldInfo* worldInfo, int node_count, const btVector3* x, const btScalar* m)
 	: m_softBodySolver(0), m_worldInfo(worldInfo)
 {
@@ -41,6 +143,7 @@ btSoftBody::btSoftBody(btSoftBodyWorldInfo* worldInfo, int node_count, const btV
 	/* Nodes			*/
 	const btScalar margin = getCollisionShape()->getMargin();
 	m_nodes.resize(node_count);
+    m_X.resize(node_count);
 	for (int i = 0, ni = node_count; i < ni; ++i)
 	{
 		Node& n = m_nodes[i];
@@ -51,8 +154,11 @@ btSoftBody::btSoftBody(btSoftBodyWorldInfo* worldInfo, int node_count, const btV
 		n.m_im = n.m_im > 0 ? 1 / n.m_im : 0;
 		n.m_leaf = m_ndbvt.insert(btDbvtVolume::FromCR(n.m_x, margin), &n);
 		n.m_material = pm;
+        m_X[i] = n.m_x;
 	}
 	updateBounds();
+	setCollisionQuadrature(3);
+	m_fdbvnt = 0;
 }
 
 btSoftBody::btSoftBody(btSoftBodyWorldInfo* worldInfo)
@@ -111,15 +217,18 @@ void btSoftBody::initDefaults()
 	m_collisionShape = new btSoftBodyCollisionShape(this);
 	m_collisionShape->setMargin(0.25f);
 
-	m_initialWorldTransform.setIdentity();
+	m_worldTransform.setIdentity();
 
 	m_windVelocity = btVector3(0, 0, 0);
 	m_restLengthScale = btScalar(1.0);
-    m_dampingCoefficient = 1;
-    m_sleepingThreshold = 0.1;
-    m_useFaceContact = true;
+	m_dampingCoefficient = 1.0;
+	m_sleepingThreshold = .4;
 	m_useSelfCollision = false;
-    m_collisionFlags = 0;
+	m_collisionFlags = 0;
+	m_softSoftCollision = false;
+	m_maxSpeedSquared = 0;
+	m_repulsionStiffness = 0.5;
+	m_fdbvnt = 0;
 }
 
 //
@@ -134,6 +243,8 @@ btSoftBody::~btSoftBody()
 		btAlignedFree(m_materials[i]);
 	for (i = 0; i < m_joints.size(); ++i)
 		btAlignedFree(m_joints[i]);
+	if (m_fdbvnt)
+		delete m_fdbvnt;
 }
 
 //
@@ -897,6 +1008,71 @@ void btSoftBody::setVolumeDensity(btScalar density)
 }
 
 //
+btVector3 btSoftBody::getLinearVelocity()
+{
+    btVector3 total_momentum = btVector3(0,0,0);
+    for (int i = 0; i < m_nodes.size(); ++i)
+    {
+        btScalar mass = m_nodes[i].m_im == 0 ? 0 : 1.0/m_nodes[i].m_im;
+        total_momentum += mass * m_nodes[i].m_v;
+    }
+    btScalar total_mass = getTotalMass();
+    return total_mass == 0 ? total_momentum : total_momentum / total_mass;
+}
+
+//
+void btSoftBody::setLinearVelocity(const btVector3& linVel)
+{
+    btVector3 old_vel = getLinearVelocity();
+    btVector3 diff = linVel - old_vel;
+    for (int i = 0; i < m_nodes.size(); ++i)
+        m_nodes[i].m_v += diff;
+}
+
+//
+void btSoftBody::setAngularVelocity(const btVector3& angVel)
+{
+    btVector3 old_vel = getLinearVelocity();
+    btVector3 com = getCenterOfMass();
+    for (int i = 0; i < m_nodes.size(); ++i)
+    {
+        m_nodes[i].m_v = angVel.cross(m_nodes[i].m_x - com) + old_vel;
+    }
+}
+
+//
+btTransform btSoftBody::getRigidTransform()
+{
+    btVector3 t = getCenterOfMass();
+    btMatrix3x3 S;
+    S.setZero();
+    // get rotation that minimizes L2 difference: \sum_i || RX_i + t - x_i ||
+    for (int i = 0; i < m_nodes.size(); ++i)
+    {
+        S += OuterProduct(m_X[i], t-m_nodes[i].m_x);
+    }
+    btVector3 sigma;
+    btMatrix3x3 U,V;
+    singularValueDecomposition(S,U,sigma,V);
+    btMatrix3x3 R = V * U.transpose();
+    btTransform trs;
+    trs.setIdentity();
+    trs.setOrigin(t);
+    trs.setBasis(R);
+    return trs;
+}
+
+//
+void btSoftBody::transformTo(const btTransform& trs)
+{
+    // get the current best rigid fit
+    btTransform current_transform = getRigidTransform();
+    // apply transform in material space
+    btTransform new_transform = trs * current_transform.inverse();
+    transform(new_transform);
+}
+
+//
 void btSoftBody::transform(const btTransform& trs)
 {
 	const btScalar margin = getCollisionShape()->getMargin();
@@ -916,7 +1092,6 @@ void btSoftBody::transform(const btTransform& trs)
 	updateNormals();
 	updateBounds();
 	updateConstants();
-	m_initialWorldTransform = trs;
 }
 
 //
@@ -1834,6 +2009,25 @@ bool btSoftBody::rayTest(const btVector3& rayFrom,
 	return (rayTest(rayFrom, rayTo, results.fraction, results.feature, results.index, false) != 0);
 }
 
+bool btSoftBody::rayFaceTest(const btVector3& rayFrom,
+                         const btVector3& rayTo,
+                         sRayCast& results)
+{
+	if (m_faces.size() == 0)
+		return false;
+	else
+	{
+    	if (m_fdbvt.empty())
+        	initializeFaceTree();
+	}
+    
+    results.body = this;
+    results.fraction = 1.f;
+    results.index = -1;
+    
+    return (rayFaceTest(rayFrom, rayTo, results.fraction,  results.index) != 0);
+}
+
 //
 void btSoftBody::setSolver(eSolverPresets::_ preset)
 {
@@ -2339,15 +2533,160 @@ int btSoftBody::rayTest(const btVector3& rayFrom, const btVector3& rayTo,
 	return (cnt);
 }
 
+int btSoftBody::rayFaceTest(const btVector3& rayFrom, const btVector3& rayTo,
+						btScalar& mint, int& index) const
+{
+	int cnt = 0;
+	{ /* Use dbvt	*/
+		RayFromToCaster collider(rayFrom, rayTo, mint);
+		
+		btDbvt::rayTest(m_fdbvt.m_root, rayFrom, rayTo, collider);
+		if (collider.m_face)
+		{
+			mint = collider.m_mint;
+			index = (int)(collider.m_face - &m_faces[0]);
+			cnt = 1;
+		}
+	}
+	return (cnt);
+}
+
+
 //
+static inline btDbvntNode* copyToDbvnt(const btDbvtNode* n)
+{
+	if (n == 0)
+		return 0;
+	btDbvntNode* root = new btDbvntNode(n);
+	if (n->isinternal())
+	{
+		btDbvntNode* c0 = copyToDbvnt(n->childs[0]);
+		root->childs[0] = c0;
+		btDbvntNode* c1 = copyToDbvnt(n->childs[1]);
+		root->childs[1] = c1;
+	}
+	return root;
+}
+
+static inline void calculateNormalCone(btDbvntNode* root)
+{
+	if (!root)
+		return;
+	if (root->isleaf())
+	{
+		const btSoftBody::Face* face = (btSoftBody::Face*)root->data;
+		root->normal = face->m_normal;
+		root->angle = 0;
+	}
+	else
+	{
+		btVector3 n0(0,0,0), n1(0,0,0);
+		btScalar a0 = 0, a1 = 0;
+		if (root->childs[0])
+		{
+			calculateNormalCone(root->childs[0]);
+			n0 = root->childs[0]->normal;
+			a0 = root->childs[0]->angle;
+		}
+		if (root->childs[1])
+		{
+			calculateNormalCone(root->childs[1]);
+			n1 = root->childs[1]->normal;
+			a1 = root->childs[1]->angle;
+		}
+		root->normal = (n0+n1).safeNormalize();
+		root->angle = btMax(a0,a1) + btAngle(n0, n1)*0.5;
+	}
+}
+
 void btSoftBody::initializeFaceTree()
 {
+	BT_PROFILE("btSoftBody::initializeFaceTree");
 	m_fdbvt.clear();
+	// create leaf nodes;
+	btAlignedObjectArray<btDbvtNode*> leafNodes;
+	leafNodes.resize(m_faces.size());
 	for (int i = 0; i < m_faces.size(); ++i)
 	{
 		Face& f = m_faces[i];
-		f.m_leaf = m_fdbvt.insert(VolumeOf(f, 0), &f);
+		ATTRIBUTE_ALIGNED16(btDbvtVolume) vol = VolumeOf(f, 0);
+		btDbvtNode* node = new (btAlignedAlloc(sizeof(btDbvtNode), 16)) btDbvtNode();
+		node->parent = NULL;
+		node->data = &f;
+		node->childs[1] = 0;
+		node->volume = vol;
+		leafNodes[i] = node;
+		f.m_leaf = node;
 	}
+	btAlignedObjectArray<btAlignedObjectArray<int> > adj;
+	adj.resize(m_faces.size());
+	// construct the adjacency list for triangles
+	for (int i = 0; i < adj.size(); ++i)
+	{
+		for (int j = i+1; j < adj.size(); ++j)
+		{
+			int dup = 0;
+			for (int k = 0; k < 3; ++k)
+			{
+				for (int l = 0; l < 3; ++l)
+				{
+					if (m_faces[i].m_n[k] == m_faces[j].m_n[l])
+					{
+						++dup;
+						break;
+					}
+				}
+				if (dup == 2)
+				{
+					adj[i].push_back(j);
+					adj[j].push_back(i);
+				}
+			}
+		}
+	}
+	m_fdbvt.m_root = buildTreeBottomUp(leafNodes, adj);
+	if (m_fdbvnt)
+		delete m_fdbvnt;
+	m_fdbvnt = copyToDbvnt(m_fdbvt.m_root);
+	updateFaceTree(false, false);
+	rebuildNodeTree();
+}
+
+//
+void btSoftBody::rebuildNodeTree()
+{
+	m_ndbvt.clear();
+	btAlignedObjectArray<btDbvtNode*> leafNodes;
+	leafNodes.resize(m_nodes.size());
+	for (int i = 0; i < m_nodes.size(); ++i)
+	{
+		Node& n = m_nodes[i];
+		ATTRIBUTE_ALIGNED16(btDbvtVolume) vol = btDbvtVolume::FromCR(n.m_x, 0);
+		btDbvtNode* node = new (btAlignedAlloc(sizeof(btDbvtNode), 16)) btDbvtNode();
+		node->parent = NULL;
+		node->data = &n;
+		node->childs[1] = 0;
+		node->volume = vol;
+		leafNodes[i] = node;
+		n.m_leaf = node;
+	}
+	btAlignedObjectArray<btAlignedObjectArray<int> > adj;
+	adj.resize(m_nodes.size());
+	btAlignedObjectArray<int> old_id;
+	old_id.resize(m_nodes.size());
+	for (int i = 0; i < m_nodes.size(); ++i)
+		old_id[i] = m_nodes[i].index;
+	for (int i = 0; i < m_nodes.size(); ++i)
+		m_nodes[i].index = i;
+	for (int i = 0; i < m_links.size(); ++i)
+	{
+		Link& l = m_links[i];
+		adj[l.m_n[0]->index].push_back(l.m_n[1]->index);
+		adj[l.m_n[1]->index].push_back(l.m_n[0]->index);
+	}
+	m_ndbvt.m_root = buildTreeBottomUp(leafNodes, adj);
+	for (int i = 0; i < m_nodes.size(); ++i)
+		m_nodes[i].index = old_id[i];
 }
 
 //
@@ -2403,10 +2742,9 @@ bool btSoftBody::checkDeformableContact(const btCollisionObjectWrapper* colObjWr
     const btCollisionObject* tmpCollisionObj = colObjWrap->getCollisionObject();
     // use the position x_{n+1}^* = x_n + dt * v_{n+1}^* where v_{n+1}^* = v_n + dtg for collision detect
     // but resolve contact at x_n
-//    btTransform wtr = (predict) ?
-//    (colObjWrap->m_preTransform != NULL ? tmpCollisionObj->getInterpolationWorldTransform()*(*colObjWrap->m_preTransform) : tmpCollisionObj->getInterpolationWorldTransform())
-//                 : colObjWrap->getWorldTransform();
-    const btTransform& wtr = colObjWrap->getWorldTransform();
+    btTransform wtr = (predict) ?
+    (colObjWrap->m_preTransform != NULL ? tmpCollisionObj->getInterpolationWorldTransform()*(*colObjWrap->m_preTransform) : tmpCollisionObj->getInterpolationWorldTransform())
+                 : colObjWrap->getWorldTransform();
 	btScalar dst =
 		m_worldInfo->m_sparsesdf.Evaluate(
 			wtr.invXform(x),
@@ -2457,7 +2795,6 @@ bool btSoftBody::checkDeformableFaceContact(const btCollisionObjectWrapper* colO
     btTransform wtr = (predict) ?
     (colObjWrap->m_preTransform != NULL ? tmpCollisionObj->getInterpolationWorldTransform()*(*colObjWrap->m_preTransform) : tmpCollisionObj->getInterpolationWorldTransform())
     : colObjWrap->getWorldTransform();
-//    const btTransform& wtr = colObjWrap->getWorldTransform();
     btScalar dst;
     
 //#define USE_QUADRATURE 1
@@ -2476,6 +2813,7 @@ bool btSoftBody::checkDeformableFaceContact(const btCollisionObjectWrapper* colO
                                           nrm,
                                           margin);
         nrm = wtr.getBasis() * nrm;
+        cti.m_colObj = colObjWrap->getCollisionObject();
         // use cached contact point
     }
     else
@@ -2492,10 +2830,11 @@ bool btSoftBody::checkDeformableFaceContact(const btCollisionObjectWrapper* colO
         contact_point = results.witnesses[0];
         getBarycentric(contact_point, f.m_n[0]->m_x, f.m_n[1]->m_x, f.m_n[2]->m_x, bary);
         nrm = results.normal;
+        cti.m_colObj = colObjWrap->getCollisionObject();
         for (int i = 0; i < 3; ++i)
             f.m_pcontact[i] = bary[i];
     }
-
+    return (dst < 0);
 #endif
 
     // use collision quadrature point
@@ -2505,7 +2844,11 @@ bool btSoftBody::checkDeformableFaceContact(const btCollisionObjectWrapper* colO
         btVector3 local_nrm;
         for (int q = 0; q < m_quads.size(); ++q)
         {
-            btVector3 p = BaryEval(f.m_n[0]->m_x, f.m_n[1]->m_x, f.m_n[2]->m_x, m_quads[q]);
+            btVector3 p;
+            if (predict)
+                p = BaryEval(f.m_n[0]->m_q, f.m_n[1]->m_q, f.m_n[2]->m_q, m_quads[q]);
+            else
+                p = BaryEval(f.m_n[0]->m_x, f.m_n[1]->m_x, f.m_n[2]->m_x, m_quads[q]);
             btScalar local_dst = m_worldInfo->m_sparsesdf.Evaluate(
                                                     wtr.invXform(p),
                                                     shp,
@@ -2513,43 +2856,83 @@ bool btSoftBody::checkDeformableFaceContact(const btCollisionObjectWrapper* colO
                                                     margin);
             if (local_dst < dst)
             {
+                if (local_dst < 0 && predict)
+                    return true;
                 dst = local_dst;
                 contact_point = p;
                 bary = m_quads[q];
-                nrm = wtr.getBasis() * local_nrm;
+                nrm = local_nrm;
+            }
+            if (!predict)
+            {
+                cti.m_colObj = colObjWrap->getCollisionObject();
+                cti.m_normal = wtr.getBasis() * nrm;
+                cti.m_offset = dst;
             }
         }
+        return (dst < 0);
     }
 #endif
     
+//    // regular face contact
+//    {
+//        btGjkEpaSolver2::sResults results;
+//        btTransform triangle_transform;
+//        triangle_transform.setIdentity();
+//        triangle_transform.setOrigin(f.m_n[0]->m_x);
+//        btTriangleShape triangle(btVector3(0,0,0), f.m_n[1]->m_x-f.m_n[0]->m_x, f.m_n[2]->m_x-f.m_n[0]->m_x);
+//        btVector3 guess(0,0,0);
+//        if (predict)
+//        {
+//            triangle_transform.setOrigin(f.m_n[0]->m_q);
+//            triangle = btTriangleShape(btVector3(0,0,0), f.m_n[1]->m_q-f.m_n[0]->m_q, f.m_n[2]->m_q-f.m_n[0]->m_q);
+//        }
+//        const btConvexShape* csh = static_cast<const btConvexShape*>(shp);
+////        btGjkEpaSolver2::SignedDistance(&triangle, triangle_transform, csh, wtr, guess, results);
+////        dst = results.distance - margin;
+////        contact_point = results.witnesses[0];
+//        btGjkEpaSolver2::Penetration(&triangle, triangle_transform, csh, wtr, guess, results);
+//        if (results.status == btGjkEpaSolver2::sResults::Separated)
+//            return false;
+//        dst = results.distance - margin;
+//        contact_point = results.witnesses[1];
+//        getBarycentric(contact_point, f.m_n[0]->m_x, f.m_n[1]->m_x, f.m_n[2]->m_x, bary);
+//        nrm = results.normal;
+//        for (int i = 0; i < 3; ++i)
+//            f.m_pcontact[i] = bary[i];
+//    }
+//
+//    if (!predict)
+//    {
+//        cti.m_colObj = colObjWrap->getCollisionObject();
+//        cti.m_normal = nrm;
+//        cti.m_offset = dst;
+//    }
+//
+    
     // regular face contact
     {
         btGjkEpaSolver2::sResults results;
         btTransform triangle_transform;
         triangle_transform.setIdentity();
-        triangle_transform.setOrigin(f.m_n[0]->m_x);
-        btTriangleShape triangle(btVector3(0,0,0), f.m_n[1]->m_x-f.m_n[0]->m_x, f.m_n[2]->m_x-f.m_n[0]->m_x);
+        triangle_transform.setOrigin(f.m_n[0]->m_q);
+        btTriangleShape triangle(btVector3(0,0,0), f.m_n[1]->m_q-f.m_n[0]->m_q, f.m_n[2]->m_q-f.m_n[0]->m_q);
         btVector3 guess(0,0,0);
         const btConvexShape* csh = static_cast<const btConvexShape*>(shp);
         btGjkEpaSolver2::SignedDistance(&triangle, triangle_transform, csh, wtr, guess, results);
-        dst = results.distance - margin;
+        dst = results.distance-csh->getMargin();
+        dst -= margin;
+        if (dst >= 0)
+            return false;
         contact_point = results.witnesses[0];
-        getBarycentric(contact_point, f.m_n[0]->m_x, f.m_n[1]->m_x, f.m_n[2]->m_x, bary);
+        getBarycentric(contact_point, f.m_n[0]->m_q, f.m_n[1]->m_q, f.m_n[2]->m_q, bary);
+        btVector3 curr = BaryEval(f.m_n[0]->m_x, f.m_n[1]->m_x, f.m_n[2]->m_x, bary);
         nrm = results.normal;
-        for (int i = 0; i < 3; ++i)
-            f.m_pcontact[i] = bary[i];
-    }
-    
-    if (!predict)
-    {
         cti.m_colObj = colObjWrap->getCollisionObject();
         cti.m_normal = nrm;
-        cti.m_offset = dst;
+        cti.m_offset = dst + (curr - contact_point).dot(nrm);
     }
-    
-    if (dst < 0)
-        return true;
-    return (false);
+    return (dst < 0);
 }
 
 //
@@ -3075,6 +3458,7 @@ void btSoftBody::setSpringStiffness(btScalar k)
     {
         m_links[i].Feature::m_material->m_kLST = k;
     }
+    m_repulsionStiffness = k;
 }
 
 void btSoftBody::initializeDmInverse()
@@ -3372,18 +3756,39 @@ void btSoftBody::setMaxStress(btScalar maxStress)
 //
 void btSoftBody::interpolateRenderMesh()
 {
-    for (int i = 0; i < m_renderNodes.size(); ++i)
-    {
-        Node& n = m_renderNodes[i];
-        n.m_x.setZero();
-        for (int j = 0; j < 4; ++j)
-        {
-			if (m_renderNodesParents[i].size())
+	if (m_z.size() > 0)
+	{
+		for (int i = 0; i < m_renderNodes.size(); ++i)
+		{
+			const Node* p0 = m_renderNodesParents[i][0];
+			const Node* p1 = m_renderNodesParents[i][1];
+			const Node* p2 = m_renderNodesParents[i][2];
+			btVector3 normal = btCross(p1->m_x - p0->m_x, p2->m_x - p0->m_x);
+			btVector3 unit_normal = normal.normalized();
+			Node& n = m_renderNodes[i];
+			n.m_x.setZero();
+			for (int j = 0; j < 3; ++j)
 			{
 				n.m_x += m_renderNodesParents[i][j]->m_x * m_renderNodesInterpolationWeights[i][j];
 			}
-        }
-    }
+			n.m_x += m_z[i] * unit_normal;
+		}
+	}
+	else
+	{
+		for (int i = 0; i < m_renderNodes.size(); ++i)
+		{
+			Node& n = m_renderNodes[i];
+			n.m_x.setZero();
+			for (int j = 0; j < 4; ++j)
+			{
+				if (m_renderNodesParents[i].size())
+				{
+					n.m_x += m_renderNodesParents[i][j]->m_x * m_renderNodesInterpolationWeights[i][j];
+				}
+			}
+		}
+	}
 }
 
 void btSoftBody::setCollisionQuadrature(int N)
@@ -3649,13 +4054,10 @@ void btSoftBody::defaultCollisionHandler(const btCollisionObjectWrapper* pcoWrap
 		break;
         case fCollision::SDF_RD:
         {
-            
             btRigidBody* prb1 = (btRigidBody*)btRigidBody::upcast(pcoWrap->getCollisionObject());
             if (pcoWrap->getCollisionObject()->isActive() || this->isActive())
             {
                 const btTransform wtr = pcoWrap->getWorldTransform();
-//                const btTransform ctr = pcoWrap->getWorldTransform();
-//                const btScalar timemargin = (wtr.getOrigin() - ctr.getOrigin()).length();
                 const btScalar timemargin = 0;
                 const btScalar basemargin = getCollisionShape()->getMargin();
                 btVector3 mins;
@@ -3667,22 +4069,25 @@ void btSoftBody::defaultCollisionHandler(const btCollisionObjectWrapper* pcoWrap
                                                       maxs);
                 volume = btDbvtVolume::FromMM(mins, maxs);
                 volume.Expand(btVector3(basemargin, basemargin, basemargin));
-                btSoftColliders::CollideSDF_RD docollideNode;
-                docollideNode.psb = this;
-                docollideNode.m_colObj1Wrap = pcoWrap;
-                docollideNode.m_rigidBody = prb1;
-                docollideNode.dynmargin = basemargin + timemargin;
-                docollideNode.stamargin = basemargin;
-                m_ndbvt.collideTV(m_ndbvt.m_root, volume, docollideNode);
-                
-                if (this->m_useFaceContact)
+				if (m_cfg.collisions & fCollision::SDF_RDN)
+				{
+					btSoftColliders::CollideSDF_RD docollideNode;
+					docollideNode.psb = this;
+					docollideNode.m_colObj1Wrap = pcoWrap;
+					docollideNode.m_rigidBody = prb1;
+					docollideNode.dynmargin = basemargin + timemargin;
+					docollideNode.stamargin = basemargin;
+					m_ndbvt.collideTV(m_ndbvt.m_root, volume, docollideNode);
+				}
+
+                if (((pcoWrap->getCollisionObject()->getInternalType() == CO_RIGID_BODY) && (m_cfg.collisions & fCollision::SDF_RDF)) || ((pcoWrap->getCollisionObject()->getInternalType() == CO_FEATHERSTONE_LINK) && (m_cfg.collisions & fCollision::SDF_MDF)))
                 {
                     btSoftColliders::CollideSDF_RDF docollideFace;
                     docollideFace.psb = this;
                     docollideFace.m_colObj1Wrap = pcoWrap;
                     docollideFace.m_rigidBody = prb1;
-                    docollideFace.dynmargin = basemargin + timemargin;
-                    docollideFace.stamargin = basemargin;
+					docollideFace.dynmargin = basemargin + timemargin;
+					docollideFace.stamargin = basemargin;
                     m_fdbvt.collideTV(m_fdbvt.m_root, volume, docollideFace);
                 }
             }
@@ -3691,51 +4096,6 @@ void btSoftBody::defaultCollisionHandler(const btCollisionObjectWrapper* pcoWrap
 	}
 }
 
-static inline btDbvntNode* copyToDbvnt(const btDbvtNode* n)
-{
-    if (n == 0)
-        return 0;
-    btDbvntNode* root = new btDbvntNode(n);
-    if (n->isinternal())
-    {
-        btDbvntNode* c0 = copyToDbvnt(n->childs[0]);
-        root->childs[0] = c0;
-        btDbvntNode* c1 = copyToDbvnt(n->childs[1]);
-        root->childs[1] = c1;
-    }
-    return root;
-}
-
-static inline void calculateNormalCone(btDbvntNode* root)
-{
-    if (!root)
-        return;
-    if (root->isleaf())
-    {
-        const btSoftBody::Face* face = (btSoftBody::Face*)root->data;
-        root->normal = face->m_normal;
-        root->angle = 0;
-    }
-    else
-    {
-        btVector3 n0(0,0,0), n1(0,0,0);
-        btScalar a0 = 0, a1 = 0;
-        if (root->childs[0])
-        {
-            calculateNormalCone(root->childs[0]);
-            n0 = root->childs[0]->normal;
-            a0 = root->childs[0]->angle;
-        }
-        if (root->childs[1])
-        {
-            calculateNormalCone(root->childs[1]);
-            n1 = root->childs[1]->normal;
-            a1 = root->childs[1]->angle;
-        }
-        root->normal = (n0+n1).safeNormalize();
-        root->angle = btMax(a0,a1) + btAngle(n0, n1)*0.5;
-    }
-}
 //
 void btSoftBody::defaultCollisionHandler(btSoftBody* psb)
 {
@@ -3779,6 +4139,8 @@ void btSoftBody::defaultCollisionHandler(btSoftBody* psb)
 		break;
         case fCollision::VF_DD:
         {
+            if (!psb->m_softSoftCollision)
+                return;
             if (psb->isActive() || this->isActive())
             {
                 if (this != psb)
@@ -3797,6 +4159,7 @@ void btSoftBody::defaultCollisionHandler(btSoftBody* psb)
                     docollide.psb[0]->m_ndbvt.collideTT(docollide.psb[0]->m_ndbvt.m_root,
                                                         docollide.psb[1]->m_fdbvt.m_root,
                                                         docollide);
+
                     /* psb1 nodes vs psb0 faces    */
                     if (this->m_tetras.size() > 0)
                         docollide.useFaceNormal = true;
@@ -3812,20 +4175,17 @@ void btSoftBody::defaultCollisionHandler(btSoftBody* psb)
                 {
                     if (psb->useSelfCollision())
                     {
-                      btSoftColliders::CollideFF_DD docollide;
-                      docollide.mrg = getCollisionShape()->getMargin() +
-                      psb->getCollisionShape()->getMargin();
-                      docollide.psb[0] = this;
-                      docollide.psb[1] = psb;
-                      if (this->m_tetras.size() > 0)
-                          docollide.useFaceNormal = true;
-                      else
-                          docollide.useFaceNormal = false;
-                      /* psb0 faces vs psb0 faces    */
-                      btDbvntNode* root = copyToDbvnt(this->m_fdbvt.m_root);
-                      calculateNormalCone(root);
-                      this->m_fdbvt.selfCollideT(root,docollide);
-                      delete root;
+                        btSoftColliders::CollideFF_DD docollide;
+                        docollide.mrg = 2*getCollisionShape()->getMargin();
+                        docollide.psb[0] = this;
+                        docollide.psb[1] = psb;
+                        if (this->m_tetras.size() > 0)
+                            docollide.useFaceNormal = true;
+                        else
+                            docollide.useFaceNormal = false;
+                        /* psb0 faces vs psb0 faces    */
+                        calculateNormalCone(this->m_fdbvnt);
+                        this->m_fdbvt.selfCollideT(m_fdbvnt,docollide);
                     }
                 }
             }
@@ -3837,6 +4197,58 @@ void btSoftBody::defaultCollisionHandler(btSoftBody* psb)
 	}
 }
 
+void btSoftBody::geometricCollisionHandler(btSoftBody* psb)
+{
+	if (psb->isActive() || this->isActive())
+	{
+		if (this != psb)
+		{
+			btSoftColliders::CollideCCD docollide;
+			/* common                    */
+			docollide.mrg = SAFE_EPSILON; // for rounding error instead of actual margin
+			docollide.dt = psb->m_sst.sdt;
+			/* psb0 nodes vs psb1 faces    */
+			if (psb->m_tetras.size() > 0)
+				docollide.useFaceNormal = true;
+			else
+				docollide.useFaceNormal = false;
+			docollide.psb[0] = this;
+			docollide.psb[1] = psb;
+			docollide.psb[0]->m_ndbvt.collideTT(docollide.psb[0]->m_ndbvt.m_root,
+												 docollide.psb[1]->m_fdbvt.m_root,
+												 docollide);
+			/* psb1 nodes vs psb0 faces    */
+			if (this->m_tetras.size() > 0)
+				docollide.useFaceNormal = true;
+			else
+				docollide.useFaceNormal = false;
+			docollide.psb[0] = psb;
+			docollide.psb[1] = this;
+			docollide.psb[0]->m_ndbvt.collideTT(docollide.psb[0]->m_ndbvt.m_root,
+												 docollide.psb[1]->m_fdbvt.m_root,
+												 docollide);
+		}
+		else
+		{
+			if (psb->useSelfCollision())
+			{
+				btSoftColliders::CollideCCD docollide;
+				docollide.mrg = SAFE_EPSILON;
+				docollide.psb[0] = this;
+				docollide.psb[1] = psb;
+                docollide.dt = psb->m_sst.sdt;
+				if (this->m_tetras.size() > 0)
+					docollide.useFaceNormal = true;
+				else
+					docollide.useFaceNormal = false;
+				/* psb0 faces vs psb0 faces    */
+				calculateNormalCone(this->m_fdbvnt);  // should compute this outside of this scope
+				this->m_fdbvt.selfCollideT(m_fdbvnt,docollide);
+			}
+		}
+	}
+}
+
 void btSoftBody::setWindVelocity(const btVector3& velocity)
 {
 	m_windVelocity = velocity;
diff --git a/thirdparty/bullet/BulletSoftBody/btSoftBody.h b/thirdparty/bullet/BulletSoftBody/btSoftBody.h
index 2b048c1118..6a55eccbd2 100644
--- a/thirdparty/bullet/BulletSoftBody/btSoftBody.h
+++ b/thirdparty/bullet/BulletSoftBody/btSoftBody.h
@@ -35,6 +35,8 @@ subject to the following restrictions:
 //#else
 #define btSoftBodyData btSoftBodyFloatData
 #define btSoftBodyDataName "btSoftBodyFloatData"
+static const btScalar  OVERLAP_REDUCTION_FACTOR = 0.1;
+static unsigned long seed = 243703;
 //#endif //BT_USE_DOUBLE_PRECISION
 
 class btBroadphaseInterface;
@@ -161,14 +163,18 @@ public:
 			RVSmask = 0x000f,  ///Rigid versus soft mask
 			SDF_RS = 0x0001,   ///SDF based rigid vs soft
 			CL_RS = 0x0002,    ///Cluster vs convex rigid vs soft
-            SDF_RD = 0x0003,   ///DF based rigid vs deformable
-            SDF_RDF = 0x0004,   ///DF based rigid vs deformable faces
+			SDF_RD = 0x0004,   ///rigid vs deformable
 
-			SVSmask = 0x00F0,  ///Rigid versus soft mask
+			SVSmask = 0x00f0,  ///Rigid versus soft mask
 			VF_SS = 0x0010,    ///Vertex vs face soft vs soft handling
 			CL_SS = 0x0020,    ///Cluster vs cluster soft vs soft handling
 			CL_SELF = 0x0040,  ///Cluster soft body self collision
-            VF_DD = 0x0050,    ///Vertex vs face soft vs soft handling
+			VF_DD = 0x0080,    ///Vertex vs face soft vs soft handling
+
+			RVDFmask = 0x0f00, /// Rigid versus deformable face mask
+			SDF_RDF = 0x0100,  /// GJK based Rigid vs. deformable face
+			SDF_MDF = 0x0200,  /// GJK based Multibody vs. deformable face
+            SDF_RDN = 0x0400,  /// SDF based Rigid vs. deformable node
 			/* presets	*/
 			Default = SDF_RS,
 			END
@@ -257,13 +263,13 @@ public:
 		btVector3 m_x;       // Position
 		btVector3 m_q;       // Previous step position/Test position
 		btVector3 m_v;       // Velocity
-        btVector3 m_vsplit;  // Temporary Velocity in addintion to velocity used in split impulse
         btVector3 m_vn;      // Previous step velocity
 		btVector3 m_f;       // Force accumulator
 		btVector3 m_n;       // Normal
 		btScalar m_im;       // 1/mass
 		btScalar m_area;     // Area
 		btDbvtNode* m_leaf;  // Leaf data
+		btScalar m_penetration;   // depth of penetration
 		int m_battach : 1;   // Attached
         int index;
 	};
@@ -289,6 +295,7 @@ public:
 		btScalar m_ra;       // Rest area
 		btDbvtNode* m_leaf;  // Leaf data
         btVector4 m_pcontact; // barycentric weights of the persistent contact
+        btVector3 m_n0, m_n1, m_vn;
         int m_index;
 	};
 	/* Tetra		*/
@@ -717,6 +724,15 @@ public:
 	/* SolverState	*/
 	struct SolverState
 	{
+		//if you add new variables, always initialize them!
+		SolverState()
+			:sdt(0),
+			isdt(0),
+			velmrg(0),
+			radmrg(0),
+			updmrg(0)
+		{
+		}
 		btScalar sdt;     // dt*timescale
 		btScalar isdt;    // 1/sdt
 		btScalar velmrg;  // velocity margin
@@ -796,22 +812,24 @@ public:
 	bool m_bUpdateRtCst;               // Update runtime constants
 	btDbvt m_ndbvt;                    // Nodes tree
 	btDbvt m_fdbvt;                    // Faces tree
+	btDbvntNode* m_fdbvnt;              // Faces tree with normals
 	btDbvt m_cdbvt;                    // Clusters tree
 	tClusterArray m_clusters;          // Clusters
-    btScalar m_dampingCoefficient;     // Damping Coefficient
-    btScalar m_sleepingThreshold;
-    btScalar m_maxSpeedSquared;
-    bool m_useFaceContact;
-    btAlignedObjectArray<btVector3> m_quads; // quadrature points for collision detection
-    
-    btAlignedObjectArray<btVector4> m_renderNodesInterpolationWeights;
-    btAlignedObjectArray<btAlignedObjectArray<const btSoftBody::Node*> > m_renderNodesParents;
-    bool m_useSelfCollision;
+	btScalar m_dampingCoefficient;     // Damping Coefficient
+	btScalar m_sleepingThreshold;
+	btScalar m_maxSpeedSquared;
+	btAlignedObjectArray<btVector3> m_quads; // quadrature points for collision detection
+	btScalar m_repulsionStiffness;
+    btAlignedObjectArray<btVector3> m_X;   // initial positions
+
+	btAlignedObjectArray<btVector4> m_renderNodesInterpolationWeights;
+	btAlignedObjectArray<btAlignedObjectArray<const btSoftBody::Node*> > m_renderNodesParents;
+	btAlignedObjectArray<btScalar> m_z; // vertical distance used in extrapolation
+	bool m_useSelfCollision;
+	bool m_softSoftCollision;
 
 	btAlignedObjectArray<bool> m_clusterConnectivity;  //cluster connectivity, for self-collision
 
-	btTransform m_initialWorldTransform;
-
 	btVector3 m_windVelocity;
 
 	btScalar m_restLengthScale;
@@ -843,11 +861,6 @@ public:
     {
         m_dampingCoefficient = damping_coeff;
     }
-    
-    void setUseFaceContact(bool useFaceContact)
-    {
-        m_useFaceContact = false;
-    }
 
 	///@todo: avoid internal softbody shape hack and move collision code to collision library
 	virtual void setCollisionShape(btCollisionShape* collisionShape)
@@ -957,6 +970,16 @@ public:
 	void setVolumeMass(btScalar mass);
 	/* Set volume density (using tetrahedrons)								*/
 	void setVolumeDensity(btScalar density);
+	/* Get the linear velocity of the center of mass                        */
+	btVector3 getLinearVelocity();
+	/* Set the linear velocity of the center of mass                        */
+	void setLinearVelocity(const btVector3& linVel);
+	/* Set the angular velocity of the center of mass                       */
+	void setAngularVelocity(const btVector3& angVel);
+    /* Get best fit rigid transform                                         */
+    btTransform getRigidTransform();
+    /* Transform to given pose                                              */
+    void transformTo(const btTransform& trs);
 	/* Transform															*/
 	void transform(const btTransform& trs);
 	/* Translate															*/
@@ -1023,6 +1046,11 @@ public:
 	bool rayTest(const btVector3& rayFrom,
 				 const btVector3& rayTo,
 				 sRayCast& results);
+	bool rayFaceTest(const btVector3& rayFrom,
+					 const btVector3& rayTo,
+					 sRayCast& results);
+	int rayFaceTest(const btVector3& rayFrom, const btVector3& rayTo,
+					btScalar& mint, int& index) const;
 	/* Solver presets														*/
 	void setSolver(eSolverPresets::_ preset);
 	/* predictMotion														*/
@@ -1120,6 +1148,7 @@ public:
 	int rayTest(const btVector3& rayFrom, const btVector3& rayTo,
 				btScalar& mint, eFeature::_& feature, int& index, bool bcountonly) const;
 	void initializeFaceTree();
+	void rebuildNodeTree();
 	btVector3 evaluateCom() const;
 	bool checkDeformableContact(const btCollisionObjectWrapper* colObjWrap, const btVector3& x, btScalar margin, btSoftBody::sCti& cti, bool predict = false) const;
     bool checkDeformableFaceContact(const btCollisionObjectWrapper* colObjWrap, Face& f, btVector3& contact_point, btVector3& bary, btScalar margin, btSoftBody::sCti& cti, bool predict = false) const;
@@ -1152,7 +1181,180 @@ public:
 	static void VSolve_Links(btSoftBody* psb, btScalar kst);
 	static psolver_t getSolver(ePSolver::_ solver);
 	static vsolver_t getSolver(eVSolver::_ solver);
+	void geometricCollisionHandler(btSoftBody* psb);
+#define SAFE_EPSILON SIMD_EPSILON*100.0
+	void updateNode(btDbvtNode* node, bool use_velocity, bool margin)
+	{
+		if (node->isleaf())
+		{
+			btSoftBody::Node* n = (btSoftBody::Node*)(node->data);
+			ATTRIBUTE_ALIGNED16(btDbvtVolume) vol;
+			btScalar pad = margin ? m_sst.radmrg : SAFE_EPSILON; // use user defined margin or margin for floating point precision
+			if (use_velocity)
+			{
+				btVector3 points[2] = {n->m_x, n->m_x + m_sst.sdt * n->m_v};
+				vol = btDbvtVolume::FromPoints(points, 2);
+				vol.Expand(btVector3(pad, pad, pad));
+			}
+			else
+			{
+				vol = btDbvtVolume::FromCR(n->m_x, pad);
+			}
+			node->volume = vol;
+			return;
+		}
+		else
+		{
+			updateNode(node->childs[0], use_velocity, margin);
+			updateNode(node->childs[1], use_velocity, margin);
+			ATTRIBUTE_ALIGNED16(btDbvtVolume) vol;
+			Merge(node->childs[0]->volume, node->childs[1]->volume, vol);
+			node->volume = vol;
+		}
+	}
+	
+    void updateNodeTree(bool use_velocity, bool margin)
+	{
+		if (m_ndbvt.m_root)
+			updateNode(m_ndbvt.m_root, use_velocity, margin);
+	}
+
+	template <class DBVTNODE> // btDbvtNode or btDbvntNode
+	void updateFace(DBVTNODE* node, bool use_velocity, bool margin)
+	{
+		if (node->isleaf())
+		{
+			btSoftBody::Face* f = (btSoftBody::Face*)(node->data);
+			btScalar pad = margin ? m_sst.radmrg : SAFE_EPSILON; // use user defined margin or margin for floating point precision
+			ATTRIBUTE_ALIGNED16(btDbvtVolume) vol;
+			if (use_velocity)
+			{
+				btVector3 points[6] = {f->m_n[0]->m_x, f->m_n[0]->m_x + m_sst.sdt * f->m_n[0]->m_v,
+					f->m_n[1]->m_x, f->m_n[1]->m_x + m_sst.sdt * f->m_n[1]->m_v,
+					f->m_n[2]->m_x, f->m_n[2]->m_x + m_sst.sdt * f->m_n[2]->m_v};
+				vol = btDbvtVolume::FromPoints(points, 6);
+			}
+			else
+			{
+				btVector3 points[3] = {f->m_n[0]->m_x,
+					f->m_n[1]->m_x,
+					f->m_n[2]->m_x};
+				vol = btDbvtVolume::FromPoints(points, 3);
+			}
+			vol.Expand(btVector3(pad, pad, pad));
+			node->volume = vol;
+			return;
+		}
+		else
+		{
+			updateFace(node->childs[0], use_velocity, margin);
+			updateFace(node->childs[1], use_velocity, margin);
+			ATTRIBUTE_ALIGNED16(btDbvtVolume) vol;
+			Merge(node->childs[0]->volume, node->childs[1]->volume, vol);
+			node->volume = vol;
+		}
+	}
+	void updateFaceTree(bool use_velocity, bool margin)
+	{
+		if (m_fdbvt.m_root)
+			updateFace(m_fdbvt.m_root, use_velocity, margin);
+		if (m_fdbvnt)
+			updateFace(m_fdbvnt, use_velocity, margin);
+	}
+
+	template <typename T>
+	static inline T BaryEval(const T& a,
+							 const T& b,
+							 const T& c,
+							 const btVector3& coord)
+	{
+		return (a * coord.x() + b * coord.y() + c * coord.z());
+	}
 
+    void applyRepulsionForce(btScalar timeStep, bool applySpringForce)
+	{
+		btAlignedObjectArray<int> indices;
+		{
+			// randomize the order of repulsive force
+			indices.resize(m_faceNodeContacts.size());
+			for (int i = 0; i < m_faceNodeContacts.size(); ++i)
+				indices[i] = i;
+#define NEXTRAND (seed = (1664525L * seed + 1013904223L) & 0xffffffff)
+			int i, ni;
+
+			for (i = 0, ni = indices.size(); i < ni; ++i)
+			{
+				btSwap(indices[i], indices[NEXTRAND % ni]);
+			}
+		}
+		for (int k = 0; k < m_faceNodeContacts.size(); ++k)
+		{
+			int i = indices[k];
+			btSoftBody::DeformableFaceNodeContact& c = m_faceNodeContacts[i];
+			btSoftBody::Node* node = c.m_node;
+			btSoftBody::Face* face = c.m_face;
+			const btVector3& w = c.m_bary;
+			const btVector3& n = c.m_normal;
+			btVector3 l = node->m_x - BaryEval(face->m_n[0]->m_x, face->m_n[1]->m_x, face->m_n[2]->m_x, w);
+			btScalar d = c.m_margin - n.dot(l);
+			d = btMax(btScalar(0),d);
+			
+			const btVector3& va = node->m_v;
+			btVector3 vb = BaryEval(face->m_n[0]->m_v, face->m_n[1]->m_v, face->m_n[2]->m_v, w);
+			btVector3 vr = va - vb;
+			const btScalar vn = btDot(vr, n); // dn < 0 <==> opposing
+			if (vn > OVERLAP_REDUCTION_FACTOR * d / timeStep)
+				continue;
+			btVector3 vt = vr - vn*n;
+			btScalar I = 0;
+			btScalar mass = node->m_im == 0 ? 0 : btScalar(1)/node->m_im;
+			if (applySpringForce)
+				I = -btMin(m_repulsionStiffness * timeStep * d, mass * (OVERLAP_REDUCTION_FACTOR * d / timeStep - vn));
+			if (vn < 0)
+				I += 0.5 * mass * vn;
+			btScalar face_penetration = 0, node_penetration = node->m_penetration;
+			for (int i = 0; i < 3; ++i)
+				face_penetration =  btMax(face_penetration, face->m_n[i]->m_penetration);
+			btScalar I_tilde = .5 *I /(1.0+w.length2());
+			
+//             double the impulse if node or face is constrained.
+            if (face_penetration > 0 || node_penetration > 0)
+                I_tilde *= 2.0;
+            if (face_penetration <= node_penetration)
+			{
+				for (int j = 0; j < 3; ++j)
+					face->m_n[j]->m_v += w[j]*n*I_tilde*node->m_im;
+			}
+            if (face_penetration >= node_penetration)
+			{
+				node->m_v -= I_tilde*node->m_im*n;
+			}
+			
+			// apply frictional impulse
+			btScalar vt_norm = vt.safeNorm();
+			if (vt_norm > SIMD_EPSILON)
+			{
+				btScalar delta_vn = -2 * I * node->m_im;
+				btScalar mu = c.m_friction;
+				btScalar vt_new = btMax(btScalar(1) - mu * delta_vn / (vt_norm + SIMD_EPSILON), btScalar(0))*vt_norm;
+				I = 0.5 * mass * (vt_norm-vt_new);
+				vt.safeNormalize();
+				I_tilde = .5 *I /(1.0+w.length2());
+//                 double the impulse if node or face is constrained.
+//                if (face_penetration > 0 || node_penetration > 0)
+//                    I_tilde *= 2.0;
+                if (face_penetration <= node_penetration)
+				{
+					for (int j = 0; j < 3; ++j)
+						face->m_n[j]->m_v += w[j] * vt * I_tilde * (face->m_n[j])->m_im;
+				}
+                if (face_penetration >= node_penetration)
+				{
+					node->m_v -= I_tilde * node->m_im * vt;
+				}
+			}
+		}
+	}
 	virtual int calculateSerializeBufferSize() const;
   
 	///fills the dataBuffer and returns the struct name (and 0 on failure)
diff --git a/thirdparty/bullet/BulletSoftBody/btSoftBodyHelpers.cpp b/thirdparty/bullet/BulletSoftBody/btSoftBodyHelpers.cpp
index 649d6f58cf..c1a87c7d57 100644
--- a/thirdparty/bullet/BulletSoftBody/btSoftBodyHelpers.cpp
+++ b/thirdparty/bullet/BulletSoftBody/btSoftBodyHelpers.cpp
@@ -1300,13 +1300,23 @@ btSoftBody* btSoftBodyHelpers::CreateFromVtkFile(btSoftBodyWorldInfo& worldInfo,
         }
         else if (reading_tets)
         {
+			int d;
+			ss >> d;
+			if (d != 4)
+			{
+				printf("Load deformable failed: Only Tetrahedra are supported in VTK file.\n");
+				fs.close();
+				return 0;
+			}
             ss.ignore(128, ' '); // ignore "4"
             Index tet;
             tet.resize(4);
             for (size_t i = 0; i < 4; i++)
             {
                 ss >> tet[i];
+                printf("%d ", tet[i]);
             }
+            printf("\n");
             indices[indices_count++] = tet;
         }
     }
@@ -1500,10 +1510,27 @@ void btSoftBodyHelpers::getBarycentricWeights(const btVector3& a, const btVector
     bary = btVector4(va6*v6, vb6*v6, vc6*v6, vd6*v6);
 }
 
+// Given a simplex with vertices a,b,c, find the barycentric weights of p in this simplex. bary[3] = 0.
+void btSoftBodyHelpers::getBarycentricWeights(const btVector3& a, const btVector3& b, const btVector3& c, const btVector3& p, btVector4& bary)
+{
+    btVector3 v0 = b - a, v1 = c - a, v2 = p - a;
+    btScalar d00 = btDot(v0, v0);
+    btScalar d01 = btDot(v0, v1);
+    btScalar d11 = btDot(v1, v1);
+    btScalar d20 = btDot(v2, v0);
+    btScalar d21 = btDot(v2, v1);
+    btScalar invDenom = 1.0 / (d00 * d11 - d01 * d01);
+    bary[1] = (d11 * d20 - d01 * d21) * invDenom;
+    bary[2] = (d00 * d21 - d01 * d20) * invDenom;
+    bary[0] = 1.0 - bary[1] - bary[2];
+    bary[3] = 0;
+}
+
 // Iterate through all render nodes to find the simulation tetrahedron that contains the render node and record the barycentric weights
 // If the node is not inside any tetrahedron, assign it to the tetrahedron in which the node has the least negative barycentric weight
 void btSoftBodyHelpers::interpolateBarycentricWeights(btSoftBody* psb)
 {
+    psb->m_z.resize(0);
     psb->m_renderNodesInterpolationWeights.resize(psb->m_renderNodes.size());
     psb->m_renderNodesParents.resize(psb->m_renderNodes.size());
     for (int i = 0; i < psb->m_renderNodes.size(); ++i)
@@ -1513,7 +1540,6 @@ void btSoftBodyHelpers::interpolateBarycentricWeights(btSoftBody* psb)
         btVector4 optimal_bary;
         btScalar min_bary_weight = -1e3;
         btAlignedObjectArray<const btSoftBody::Node*> optimal_parents;
-        bool found = false;
         for (int j = 0; j < psb->m_tetras.size(); ++j)
         {
             const btSoftBody::Tetra& t = psb->m_tetras[j];
@@ -1544,3 +1570,55 @@ void btSoftBodyHelpers::interpolateBarycentricWeights(btSoftBody* psb)
         psb->m_renderNodesParents[i] = optimal_parents;
     }
 }
+
+
+// Iterate through all render nodes to find the simulation triangle that's closest to the node in the barycentric sense.
+void btSoftBodyHelpers::extrapolateBarycentricWeights(btSoftBody* psb)
+{
+    psb->m_renderNodesInterpolationWeights.resize(psb->m_renderNodes.size());
+    psb->m_renderNodesParents.resize(psb->m_renderNodes.size());
+    psb->m_z.resize(psb->m_renderNodes.size());
+    for (int i = 0; i < psb->m_renderNodes.size(); ++i)
+    {
+        const btVector3& p = psb->m_renderNodes[i].m_x;
+        btVector4 bary;
+        btVector4 optimal_bary;
+        btScalar min_bary_weight = -SIMD_INFINITY;
+        btAlignedObjectArray<const btSoftBody::Node*> optimal_parents;
+        btScalar dist = 0, optimal_dist = 0;
+        for (int j = 0; j < psb->m_faces.size(); ++j)
+        {
+            const btSoftBody::Face& f = psb->m_faces[j];
+            btVector3 n = btCross(f.m_n[1]->m_x - f.m_n[0]->m_x,  f.m_n[2]->m_x - f.m_n[0]->m_x);
+            btVector3 unit_n = n.normalized();
+            dist = (p-f.m_n[0]->m_x).dot(unit_n);
+            btVector3 proj_p = p - dist*unit_n;
+            getBarycentricWeights(f.m_n[0]->m_x, f.m_n[1]->m_x, f.m_n[2]->m_x, proj_p, bary);
+            btScalar new_min_bary_weight = bary[0];
+            for (int k = 1; k < 3; ++k)
+            {
+                new_min_bary_weight = btMin(new_min_bary_weight, bary[k]);
+            }
+
+            // p is out of the current best triangle, we found a traingle that's better
+            bool better_than_closest_outisde = (new_min_bary_weight > min_bary_weight && min_bary_weight<0.);
+            // p is inside of the current best triangle, we found a triangle that's better
+            bool better_than_best_inside = (new_min_bary_weight>=0 &&  min_bary_weight>=0 && btFabs(dist)<btFabs(optimal_dist));
+
+            if (better_than_closest_outisde || better_than_best_inside)
+            {
+                btAlignedObjectArray<const btSoftBody::Node*> parents;
+                parents.push_back(f.m_n[0]);
+                parents.push_back(f.m_n[1]);
+                parents.push_back(f.m_n[2]);
+                optimal_parents = parents;
+                optimal_bary = bary;
+                optimal_dist = dist;
+                min_bary_weight = new_min_bary_weight;
+            }
+        }
+        psb->m_renderNodesInterpolationWeights[i] = optimal_bary;
+        psb->m_renderNodesParents[i] = optimal_parents;
+        psb->m_z[i] = optimal_dist;
+    }
+}
diff --git a/thirdparty/bullet/BulletSoftBody/btSoftBodyHelpers.h b/thirdparty/bullet/BulletSoftBody/btSoftBodyHelpers.h
index b20f2f6d62..abe1870890 100644
--- a/thirdparty/bullet/BulletSoftBody/btSoftBodyHelpers.h
+++ b/thirdparty/bullet/BulletSoftBody/btSoftBodyHelpers.h
@@ -148,8 +148,12 @@ struct btSoftBodyHelpers
     
     static void getBarycentricWeights(const btVector3& a, const btVector3& b, const btVector3& c, const btVector3& d, const btVector3& p, btVector4& bary);
     
+    static void getBarycentricWeights(const btVector3& a, const btVector3& b, const btVector3& c, const btVector3& p, btVector4& bary);
+    
     static void interpolateBarycentricWeights(btSoftBody* psb);
     
+    static void extrapolateBarycentricWeights(btSoftBody* psb);
+    
     static void generateBoundaryFaces(btSoftBody* psb);
     
     static void duplicateFaces(const char* filename, const btSoftBody* psb);
diff --git a/thirdparty/bullet/BulletSoftBody/btSoftBodyInternals.h b/thirdparty/bullet/BulletSoftBody/btSoftBodyInternals.h
index cde4746d58..b9ebc95b6b 100644
--- a/thirdparty/bullet/BulletSoftBody/btSoftBodyInternals.h
+++ b/thirdparty/bullet/BulletSoftBody/btSoftBodyInternals.h
@@ -18,7 +18,6 @@ subject to the following restrictions:
 #define _BT_SOFT_BODY_INTERNALS_H
 
 #include "btSoftBody.h"
-
 #include "LinearMath/btQuickprof.h"
 #include "LinearMath/btPolarDecomposition.h"
 #include "BulletCollision/BroadphaseCollision/btBroadphaseInterface.h"
@@ -29,9 +28,10 @@ subject to the following restrictions:
 #include "BulletDynamics/Featherstone/btMultiBodyConstraint.h"
 #include <string.h>  //for memset
 #include <cmath>
+#include "poly34.h"
 
 // Given a multibody link, a contact point and a contact direction, fill in the jacobian data needed to calculate the velocity change given an impulse in the contact direction
-static void findJacobian(const btMultiBodyLinkCollider* multibodyLinkCol,
+static SIMD_FORCE_INLINE void findJacobian(const btMultiBodyLinkCollider* multibodyLinkCol,
                          btMultiBodyJacobianData& jacobianData,
                          const btVector3& contact_point,
                          const btVector3& dir)
@@ -44,7 +44,7 @@ static void findJacobian(const btMultiBodyLinkCollider* multibodyLinkCol,
     multibodyLinkCol->m_multiBody->fillContactJacobianMultiDof(multibodyLinkCol->m_link, contact_point, dir, jac, jacobianData.scratch_r, jacobianData.scratch_v, jacobianData.scratch_m);
     multibodyLinkCol->m_multiBody->calcAccelerationDeltasMultiDof(&jacobianData.m_jacobians[0], &jacobianData.m_deltaVelocitiesUnitImpulse[0], jacobianData.scratch_r, jacobianData.scratch_v);
 }
-static btVector3 generateUnitOrthogonalVector(const btVector3& u)
+static SIMD_FORCE_INLINE btVector3 generateUnitOrthogonalVector(const btVector3& u)
 {
     btScalar ux = u.getX();
     btScalar uy = u.getY();
@@ -62,6 +62,571 @@ static btVector3 generateUnitOrthogonalVector(const btVector3& u)
     v.normalize();
     return v;
 }
+
+static SIMD_FORCE_INLINE bool proximityTest(const btVector3& x1, const btVector3& x2, const btVector3& x3, const btVector3& x4, const btVector3& normal, const btScalar& mrg, btVector3& bary)
+{
+    btVector3 x43 = x4-x3;
+    if (std::abs(x43.dot(normal)) > mrg)
+        return false;
+    btVector3 x13 = x1-x3;
+    btVector3 x23 = x2-x3;
+    btScalar a11 = x13.length2();
+    btScalar a22 = x23.length2();
+    btScalar a12 = x13.dot(x23);
+    btScalar b1 = x13.dot(x43);
+    btScalar b2 = x23.dot(x43);
+    btScalar det = a11*a22 - a12*a12;
+    if (det < SIMD_EPSILON)
+        return false;
+    btScalar w1 = (b1*a22-b2*a12)/det;
+    btScalar w2 = (b2*a11-b1*a12)/det;
+    btScalar w3 = 1-w1-w2;
+    btScalar delta = mrg / std::sqrt(0.5*std::abs(x13.cross(x23).safeNorm()));
+    bary = btVector3(w1,w2,w3);
+    for (int i = 0; i < 3; ++i)
+    {
+        if (bary[i] < -delta || bary[i] > 1+delta)
+            return false;
+    }
+    return true;
+}
+static const int KDOP_COUNT = 13;
+static btVector3 dop[KDOP_COUNT]={btVector3(1,0,0),
+	btVector3(0,1,0),
+	btVector3(0,0,1),
+	btVector3(1,1,0),
+	btVector3(1,0,1),
+	btVector3(0,1,1),
+	btVector3(1,-1,0),
+	btVector3(1,0,-1),
+	btVector3(0,1,-1),
+	btVector3(1,1,1),
+	btVector3(1,-1,1),
+	btVector3(1,1,-1),
+	btVector3(1,-1,-1)
+};
+
+static inline int getSign(const btVector3& n, const btVector3& x)
+{
+	btScalar d = n.dot(x);
+	if (d>SIMD_EPSILON)
+		return 1;
+	if (d<-SIMD_EPSILON)
+		return -1;
+	return 0;
+}
+
+static SIMD_FORCE_INLINE bool hasSeparatingPlane(const btSoftBody::Face* face, const btSoftBody::Node* node, const btScalar& dt)
+{
+	btVector3 hex[6] = {face->m_n[0]->m_x - node->m_x,
+		face->m_n[1]->m_x - node->m_x,
+		face->m_n[2]->m_x - node->m_x,
+		face->m_n[0]->m_x + dt*face->m_n[0]->m_v - node->m_x,
+		face->m_n[1]->m_x + dt*face->m_n[1]->m_v - node->m_x,
+		face->m_n[2]->m_x + dt*face->m_n[2]->m_v - node->m_x
+	};
+	btVector3 segment = dt*node->m_v;
+	for (int i = 0; i < KDOP_COUNT; ++i)
+	{
+		int s = getSign(dop[i], segment);
+		int j = 0;
+		for (; j < 6; ++j)
+		{
+			if (getSign(dop[i], hex[j]) == s)
+				break;
+		}
+		if (j == 6)
+			return true;
+	}
+	return false;
+}
+
+static SIMD_FORCE_INLINE bool nearZero(const btScalar& a)
+{
+    return (a>-SAFE_EPSILON && a<SAFE_EPSILON);
+}
+static SIMD_FORCE_INLINE bool sameSign(const btScalar& a, const btScalar& b)
+{
+    return (nearZero(a) || nearZero(b) || (a>SAFE_EPSILON && b>SAFE_EPSILON) || (a<-SAFE_EPSILON && b<-SAFE_EPSILON));
+}
+static SIMD_FORCE_INLINE bool diffSign(const btScalar& a, const btScalar& b)
+{
+    return !sameSign(a, b);
+}
+inline btScalar evaluateBezier2(const btScalar &p0, const btScalar &p1, const btScalar &p2, const btScalar &t, const btScalar &s)
+{
+    btScalar s2 = s*s;
+    btScalar t2 = t*t;
+    
+    return p0*s2+p1*btScalar(2.0)*s*t+p2*t2;
+}
+inline btScalar evaluateBezier(const btScalar &p0, const btScalar &p1, const btScalar &p2, const btScalar &p3, const btScalar &t, const btScalar &s)
+{
+    btScalar s2 = s*s;
+    btScalar s3 = s2*s;
+    btScalar t2 = t*t;
+    btScalar t3 = t2*t;
+    
+    return p0*s3+p1*btScalar(3.0)*s2*t+p2*btScalar(3.0)*s*t2+p3*t3;
+}
+static SIMD_FORCE_INLINE bool getSigns(bool type_c, const btScalar& k0, const btScalar& k1, const btScalar& k2, const btScalar& k3, const btScalar& t0, const btScalar& t1, btScalar &lt0, btScalar &lt1)
+{
+    if (sameSign(t0, t1)) {
+        lt0 = t0;
+        lt1 = t0;
+        return true;
+    }
+    
+    if (type_c || diffSign(k0, k3)) {
+        btScalar ft = evaluateBezier(k0, k1, k2, k3, t0, -t1);
+        if (t0<-0)
+            ft = -ft;
+        
+        if (sameSign(ft, k0)) {
+            lt0 = t1;
+            lt1 = t1;
+        }
+        else {
+            lt0 = t0;
+            lt1 = t0;
+        }
+        return true;
+    }
+    
+    if (!type_c) {
+        btScalar ft = evaluateBezier(k0, k1, k2, k3, t0, -t1);
+        if (t0<-0)
+            ft = -ft;
+        
+        if (diffSign(ft, k0)) {
+            lt0 = t0;
+            lt1 = t1;
+            return true;
+        }
+        
+        btScalar fk = evaluateBezier2(k1-k0, k2-k1, k3-k2, t0, -t1);
+        
+        if (sameSign(fk, k1-k0))
+            lt0 = lt1 = t1;
+        else
+            lt0 = lt1 = t0;
+        
+        return true;
+    }
+    return false;
+}
+
+static SIMD_FORCE_INLINE void getBernsteinCoeff(const btSoftBody::Face* face, const btSoftBody::Node* node, const btScalar& dt, btScalar& k0, btScalar& k1, btScalar& k2, btScalar& k3)
+{
+    const btVector3& n0 = face->m_n0;
+    const btVector3& n1 = face->m_n1;
+    btVector3 n_hat = n0 + n1 - face->m_vn;
+    btVector3 p0ma0 = node->m_x - face->m_n[0]->m_x;
+    btVector3 p1ma1 = node->m_q - face->m_n[0]->m_q;
+    k0 = (p0ma0).dot(n0) * 3.0;
+    k1 = (p0ma0).dot(n_hat) + (p1ma1).dot(n0);
+    k2 = (p1ma1).dot(n_hat) + (p0ma0).dot(n1);
+    k3 = (p1ma1).dot(n1) * 3.0;
+}
+
+static SIMD_FORCE_INLINE void polyDecomposition(const btScalar& k0, const btScalar& k1, const btScalar& k2, const btScalar& k3, const btScalar& j0, const btScalar& j1, const btScalar& j2, btScalar& u0, btScalar& u1, btScalar& v0, btScalar& v1)
+{
+    btScalar denom = 4.0 * (j1-j2) * (j1-j0) + (j2-j0) * (j2-j0);
+    u0 = (2.0*(j1-j2)*(3.0*k1-2.0*k0-k3) - (j0-j2)*(3.0*k2-2.0*k3-k0)) / denom;
+    u1 = (2.0*(j1-j0)*(3.0*k2-2.0*k3-k0) - (j2-j0)*(3.0*k1-2.0*k0-k3)) / denom;
+    v0 = k0-u0*j0;
+    v1 = k3-u1*j2;
+}
+
+static SIMD_FORCE_INLINE bool rootFindingLemma(const btScalar& k0, const btScalar& k1, const btScalar& k2, const btScalar& k3)
+{
+    btScalar u0, u1, v0, v1;
+    btScalar j0 = 3.0*(k1-k0);
+    btScalar j1 = 3.0*(k2-k1);
+    btScalar j2 = 3.0*(k3-k2);
+    polyDecomposition(k0,k1,k2,k3,j0,j1,j2,u0,u1,v0,v1);
+    if (sameSign(v0, v1))
+    {
+        btScalar Ypa = j0*(1.0-v0)*(1.0-v0) + 2.0*j1*v0*(1.0-v0) + j2*v0*v0; // Y'(v0)
+        if (sameSign(Ypa, j0))
+        {
+            return (diffSign(k0,v1));
+        }
+    }
+    return diffSign(k0,v0);
+}
+
+static SIMD_FORCE_INLINE void getJs(const btScalar& k0, const btScalar& k1, const btScalar& k2, const btScalar& k3, const btSoftBody::Node* a, const btSoftBody::Node* b, const btSoftBody::Node* c, const btSoftBody::Node* p, const btScalar& dt, btScalar& j0,  btScalar& j1, btScalar& j2)
+{
+    const btVector3& a0 = a->m_x;
+    const btVector3& b0 = b->m_x;
+    const btVector3& c0 = c->m_x;
+    const btVector3& va = a->m_v;
+    const btVector3& vb = b->m_v;
+    const btVector3& vc = c->m_v;
+    const btVector3 a1 = a0 + dt*va;
+    const btVector3 b1 = b0 + dt*vb;
+    const btVector3 c1 = c0 + dt*vc;
+    btVector3 n0 = (b0-a0).cross(c0-a0);
+    btVector3 n1 = (b1-a1).cross(c1-a1);
+    btVector3 n_hat = n0+n1 - dt*dt*(vb-va).cross(vc-va);
+    const btVector3& p0 = p->m_x;
+    const btVector3& vp = p->m_v;
+    btVector3 p1 = p0 + dt*vp;
+    btVector3 m0 = (b0-p0).cross(c0-p0);
+    btVector3 m1 = (b1-p1).cross(c1-p1);
+    btVector3 m_hat = m0+m1 - dt*dt*(vb-vp).cross(vc-vp);
+    btScalar l0 = m0.dot(n0);
+    btScalar l1 = 0.25 * (m0.dot(n_hat) + m_hat.dot(n0));
+    btScalar l2 = btScalar(1)/btScalar(6)*(m0.dot(n1) + m_hat.dot(n_hat) + m1.dot(n0));
+    btScalar l3 = 0.25 * (m_hat.dot(n1) + m1.dot(n_hat));
+    btScalar l4 = m1.dot(n1);
+    
+    btScalar k1p = 0.25 * k0 + 0.75 * k1;
+    btScalar k2p = 0.5 * k1 + 0.5 * k2;
+    btScalar k3p = 0.75 * k2 + 0.25 * k3;
+    
+    btScalar s0 = (l1 * k0 - l0 * k1p)*4.0;
+    btScalar s1 = (l2 * k0 - l0 * k2p)*2.0;
+    btScalar s2 = (l3 * k0 - l0 * k3p)*btScalar(4)/btScalar(3);
+    btScalar s3 = l4 * k0 - l0 * k3;
+    
+    j0 = (s1*k0 - s0*k1) * 3.0;
+    j1 = (s2*k0 - s0*k2) * 1.5;
+    j2 = (s3*k0 - s0*k3);
+}
+
+static SIMD_FORCE_INLINE bool signDetermination1Internal(const btScalar& k0, const btScalar& k1, const btScalar& k2, const btScalar& k3, const btScalar& u0, const btScalar& u1, const btScalar& v0, const btScalar& v1)
+{
+    btScalar Yu0 = k0*(1.0-u0)*(1.0-u0)*(1.0-u0) + 3.0*k1*u0*(1.0-u0)*(1.0-u0) + 3.0*k2*u0*u0*(1.0-u0) + k3*u0*u0*u0; // Y(u0)
+    btScalar Yv0 = k0*(1.0-v0)*(1.0-v0)*(1.0-v0) + 3.0*k1*v0*(1.0-v0)*(1.0-v0) + 3.0*k2*v0*v0*(1.0-v0) + k3*v0*v0*v0; // Y(v0)
+
+    btScalar sign_Ytp = (u0 > u1) ? Yu0 : -Yu0;
+    btScalar L = sameSign(sign_Ytp, k0) ? u1 : u0;
+    sign_Ytp = (v0 > v1) ? Yv0 : -Yv0;
+    btScalar K = (sameSign(sign_Ytp,k0)) ? v1 : v0;
+    return diffSign(L,K);
+}
+
+static SIMD_FORCE_INLINE bool signDetermination2Internal(const btScalar& k0, const btScalar& k1, const btScalar& k2, const btScalar& k3, const btScalar& j0, const btScalar& j1, const btScalar& j2, const btScalar& u0, const btScalar& u1, const btScalar& v0, const btScalar& v1)
+{
+    btScalar Yu0 = k0*(1.0-u0)*(1.0-u0)*(1.0-u0) + 3.0*k1*u0*(1.0-u0)*(1.0-u0) + 3.0*k2*u0*u0*(1.0-u0) + k3*u0*u0*u0; // Y(u0)
+    btScalar sign_Ytp = (u0 > u1) ? Yu0 : -Yu0, L1, L2;
+    if (diffSign(sign_Ytp,k0))
+    {
+        L1 = u0;
+        L2 = u1;
+    }
+    else
+    {
+        btScalar Yp_u0 = j0*(1.0-u0)*(1.0-u0) + 2.0*j1*(1.0-u0)*u0 + j2*u0*u0;
+        if (sameSign(Yp_u0,j0))
+        {
+            L1 = u1;
+            L2 = u1;
+        }
+        else
+        {
+            L1 = u0;
+            L2 = u0;
+        }
+    }
+    btScalar Yv0 = k0*(1.0-v0)*(1.0-v0)*(1.0-v0) + 3.0*k1*v0*(1.0-v0)*(1.0-v0) + 3.0*k2*v0*v0*(1.0-v0) + k3*v0*v0*v0; // Y(uv0)
+    sign_Ytp = (v0 > v1) ? Yv0 : -Yv0;
+    btScalar K1, K2;
+    if (diffSign(sign_Ytp,k0))
+    {
+        K1 = v0;
+        K2 = v1;
+    }
+    else
+    {
+        btScalar Yp_v0 = j0*(1.0-v0)*(1.0-v0) + 2.0*j1*(1.0-v0)*v0 + j2*v0*v0;
+        if (sameSign(Yp_v0,j0))
+        {
+            K1 = v1;
+            K2 = v1;
+        }
+        else
+        {
+            K1 = v0;
+            K2 = v0;
+        }
+    }
+    return (diffSign(K1, L1) || diffSign(L2, K2));
+}
+
+static SIMD_FORCE_INLINE bool signDetermination1(const btScalar& k0, const btScalar& k1, const btScalar& k2, const btScalar& k3, const btSoftBody::Face* face, const btSoftBody::Node* node, const btScalar& dt)
+{
+    btScalar j0, j1, j2, u0, u1, v0, v1;
+    // p1
+    getJs(k0,k1,k2,k3,face->m_n[0], face->m_n[1], face->m_n[2], node, dt, j0, j1, j2);
+    if (nearZero(j0+j2-j1*2.0))
+    {
+        btScalar lt0, lt1;
+        getSigns(true, k0, k1, k2, k3, j0, j2, lt0, lt1);
+        if (lt0 < -SAFE_EPSILON)
+            return false;
+    }
+    else
+    {
+        polyDecomposition(k0,k1,k2,k3,j0,j1,j2,u0,u1,v0,v1);
+        if (!signDetermination1Internal(k0,k1,k2,k3,u0,u1,v0,v1))
+            return false;
+    }
+    // p2
+    getJs(k0,k1,k2,k3,face->m_n[1], face->m_n[2], face->m_n[0], node, dt, j0, j1, j2);
+    if (nearZero(j0+j2-j1*2.0))
+    {
+        btScalar lt0, lt1;
+        getSigns(true, k0, k1, k2, k3, j0, j2, lt0, lt1);
+        if (lt0 < -SAFE_EPSILON)
+            return false;
+    }
+    else
+    {
+        polyDecomposition(k0,k1,k2,k3,j0,j1,j2,u0,u1,v0,v1);
+        if (!signDetermination1Internal(k0,k1,k2,k3,u0,u1,v0,v1))
+            return false;
+    }
+    // p3
+    getJs(k0,k1,k2,k3,face->m_n[2], face->m_n[0], face->m_n[1], node, dt, j0, j1, j2);
+    if (nearZero(j0+j2-j1*2.0))
+    {
+        btScalar lt0, lt1;
+        getSigns(true, k0, k1, k2, k3, j0, j2, lt0, lt1);
+        if (lt0 < -SAFE_EPSILON)
+            return false;
+    }
+    else
+    {
+        polyDecomposition(k0,k1,k2,k3,j0,j1,j2,u0,u1,v0,v1);
+        if (!signDetermination1Internal(k0,k1,k2,k3,u0,u1,v0,v1))
+            return false;
+    }
+    return true;
+}
+
+static SIMD_FORCE_INLINE bool signDetermination2(const btScalar& k0, const btScalar& k1, const btScalar& k2, const btScalar& k3, const btSoftBody::Face* face, const btSoftBody::Node* node, const btScalar& dt)
+{
+    btScalar j0, j1, j2, u0, u1, v0, v1;
+    // p1
+    getJs(k0,k1,k2,k3,face->m_n[0], face->m_n[1], face->m_n[2], node, dt, j0, j1, j2);
+    if (nearZero(j0+j2-j1*2.0))
+    {
+        btScalar lt0, lt1;
+        bool bt0 = true, bt1=true;
+        getSigns(false, k0, k1, k2, k3, j0, j2, lt0, lt1);
+        if (lt0 < -SAFE_EPSILON)
+            bt0 = false;
+        if (lt1 < -SAFE_EPSILON)
+            bt1 = false;
+        if (!bt0 && !bt1)
+            return false;
+    }
+    else
+    {
+        polyDecomposition(k0,k1,k2,k3,j0,j1,j2,u0,u1,v0,v1);
+        if (!signDetermination2Internal(k0,k1,k2,k3,j0,j1,j2,u0,u1,v0,v1))
+            return false;
+    }
+    // p2
+    getJs(k0,k1,k2,k3,face->m_n[1], face->m_n[2], face->m_n[0], node, dt, j0, j1, j2);
+    if (nearZero(j0+j2-j1*2.0))
+    {
+        btScalar lt0, lt1;
+        bool bt0=true, bt1=true;
+        getSigns(false, k0, k1, k2, k3, j0, j2, lt0, lt1);
+        if (lt0 < -SAFE_EPSILON)
+            bt0 = false;
+        if (lt1 < -SAFE_EPSILON)
+            bt1 = false;
+        if (!bt0 && !bt1)
+            return false;
+    }
+    else
+    {
+        polyDecomposition(k0,k1,k2,k3,j0,j1,j2,u0,u1,v0,v1);
+        if (!signDetermination2Internal(k0,k1,k2,k3,j0,j1,j2,u0,u1,v0,v1))
+            return false;
+    }
+    // p3
+    getJs(k0,k1,k2,k3,face->m_n[2], face->m_n[0], face->m_n[1], node, dt, j0, j1, j2);
+    if (nearZero(j0+j2-j1*2.0))
+    {
+        btScalar lt0, lt1;
+        bool bt0=true, bt1=true;
+        getSigns(false, k0, k1, k2, k3, j0, j2, lt0, lt1);
+        if (lt0 < -SAFE_EPSILON)
+            bt0 = false;
+        if (lt1 < -SAFE_EPSILON)
+            bt1 = false;
+        if (!bt0 && !bt1)
+            return false;
+    }
+    else
+    {
+        polyDecomposition(k0,k1,k2,k3,j0,j1,j2,u0,u1,v0,v1);
+        if (!signDetermination2Internal(k0,k1,k2,k3,j0,j1,j2,u0,u1,v0,v1))
+            return false;
+    }
+    return true;
+}
+
+static SIMD_FORCE_INLINE bool coplanarAndInsideTest(const btScalar& k0, const btScalar& k1, const btScalar& k2, const btScalar& k3, const btSoftBody::Face* face, const btSoftBody::Node* node, const btScalar& dt)
+{
+    // Coplanar test
+    if (diffSign(k1-k0, k3-k2))
+    {
+        // Case b:
+        if (sameSign(k0, k3) && !rootFindingLemma(k0,k1,k2,k3))
+            return false;
+        // inside test
+        return signDetermination2(k0, k1, k2, k3, face, node, dt);
+    }
+    else
+    {
+        // Case c:
+        if (sameSign(k0, k3))
+            return false;
+        // inside test
+        return signDetermination1(k0, k1, k2, k3, face, node, dt);
+    }
+    return false;
+}
+static SIMD_FORCE_INLINE bool conservativeCulling(const btScalar& k0, const btScalar& k1, const btScalar& k2, const btScalar& k3, const btScalar& mrg)
+{
+    if (k0 > mrg && k1 > mrg && k2 > mrg && k3 > mrg)
+        return true;
+    if (k0 < -mrg && k1 < -mrg && k2 < -mrg && k3 < -mrg)
+        return true;
+    return false;
+}
+
+static SIMD_FORCE_INLINE bool bernsteinVFTest(const btScalar& k0, const btScalar& k1, const btScalar& k2, const btScalar& k3, const btScalar& mrg, const btSoftBody::Face* face, const btSoftBody::Node* node, const btScalar& dt)
+{
+    if (conservativeCulling(k0, k1, k2, k3, mrg))
+        return false;
+    return coplanarAndInsideTest(k0, k1, k2, k3, face, node, dt);
+}
+
+static SIMD_FORCE_INLINE void deCasteljau(const btScalar& k0, const btScalar& k1, const btScalar& k2, const btScalar& k3, const btScalar& t0, btScalar& k10, btScalar& k20, btScalar& k30, btScalar& k21, btScalar& k12)
+{
+    k10 = k0*(1.0-t0) + k1*t0;
+    btScalar k11 = k1*(1.0-t0) + k2*t0;
+    k12 = k2*(1.0-t0) + k3*t0;
+    k20 = k10*(1.0-t0) + k11*t0;
+    k21 = k11*(1.0-t0) + k12*t0;
+    k30 = k20*(1.0-t0) + k21*t0;
+}
+static SIMD_FORCE_INLINE bool bernsteinVFTest(const btSoftBody::Face* face, const btSoftBody::Node* node, const btScalar& dt, const btScalar& mrg)
+{
+    btScalar k0, k1, k2, k3;
+    getBernsteinCoeff(face, node, dt, k0, k1, k2, k3);
+    if (conservativeCulling(k0, k1, k2, k3, mrg))
+        return false;
+    return true;
+    if (diffSign(k2-2.0*k1+k0, k3-2.0*k2+k1))
+    {
+        btScalar k10, k20, k30, k21, k12;
+        btScalar t0 = (k2-2.0*k1+k0)/(k0-3.0*k1+3.0*k2-k3);
+        deCasteljau(k0, k1, k2, k3, t0, k10, k20, k30, k21, k12);
+        return bernsteinVFTest(k0, k10, k20, k30, mrg, face, node, dt) || bernsteinVFTest(k30, k21, k12, k3, mrg, face, node, dt);
+    }
+    return coplanarAndInsideTest(k0, k1, k2, k3, face, node, dt);
+}
+
+static SIMD_FORCE_INLINE bool continuousCollisionDetection(const btSoftBody::Face* face, const btSoftBody::Node* node, const btScalar& dt, const btScalar& mrg, btVector3& bary)
+{
+    if (hasSeparatingPlane(face, node, dt))
+        return false;
+    btVector3 x21 = face->m_n[1]->m_x - face->m_n[0]->m_x;
+    btVector3 x31 = face->m_n[2]->m_x - face->m_n[0]->m_x;
+    btVector3 x41 = node->m_x - face->m_n[0]->m_x;
+    btVector3 v21 = face->m_n[1]->m_v - face->m_n[0]->m_v;
+    btVector3 v31 = face->m_n[2]->m_v - face->m_n[0]->m_v;
+    btVector3 v41 = node->m_v - face->m_n[0]->m_v;
+    btVector3 a = x21.cross(x31);
+    btVector3 b = x21.cross(v31) + v21.cross(x31);
+    btVector3 c = v21.cross(v31);
+    btVector3 d = x41;
+    btVector3 e = v41;
+    btScalar a0 = a.dot(d);
+    btScalar a1 = a.dot(e) + b.dot(d);
+    btScalar a2 = c.dot(d) + b.dot(e);
+    btScalar a3 = c.dot(e);
+    btScalar eps = SAFE_EPSILON;
+    int num_roots = 0;
+    btScalar roots[3];
+    if (std::abs(a3) < eps)
+    {
+        // cubic term is zero
+        if (std::abs(a2) < eps)
+        {
+            if (std::abs(a1) < eps)
+            {
+                if (std::abs(a0) < eps)
+                {
+                    num_roots = 2;
+                    roots[0] = 0;
+                    roots[1] = dt;
+                }
+            }
+            else
+            {
+                num_roots = 1;
+                roots[0] = -a0/a1;
+            }
+        }
+        else
+        {
+            num_roots = SolveP2(roots, a1/a2, a0/a2);
+        }
+    }
+    else
+    {
+        num_roots = SolveP3(roots, a2/a3, a1/a3, a0/a3);
+    }
+//    std::sort(roots, roots+num_roots);
+    if (num_roots > 1)
+    {
+        if (roots[0] > roots[1])
+            btSwap(roots[0], roots[1]);
+    }
+    if (num_roots > 2)
+    {
+        if (roots[0] > roots[2])
+            btSwap(roots[0], roots[2]);
+        if (roots[1] > roots[2])
+            btSwap(roots[1], roots[2]);
+    }
+    for (int r = 0; r < num_roots; ++r)
+    {
+        double root = roots[r];
+        if (root <= 0)
+            continue;
+        if (root > dt + SIMD_EPSILON)
+            return false;
+        btVector3 x1 = face->m_n[0]->m_x + root * face->m_n[0]->m_v;
+        btVector3 x2 = face->m_n[1]->m_x + root * face->m_n[1]->m_v;
+        btVector3 x3 = face->m_n[2]->m_x + root * face->m_n[2]->m_v;
+        btVector3 x4 = node->m_x + root * node->m_v;
+        btVector3 normal = (x2-x1).cross(x3-x1);
+        normal.safeNormalize();
+        if (proximityTest(x1, x2, x3, x4, normal, mrg, bary))
+            return true;
+    }
+    return false;
+}
+static SIMD_FORCE_INLINE bool bernsteinCCD(const btSoftBody::Face* face, const btSoftBody::Node* node, const btScalar& dt, const btScalar& mrg, btVector3& bary)
+{
+    if (!bernsteinVFTest(face, node, dt, mrg))
+        return false;
+    if (!continuousCollisionDetection(face, node, dt, 1e-6, bary))
+        return false;
+    return true;
+}
+
 //
 // btSymMatrix
 //
@@ -373,6 +938,26 @@ static inline btMatrix3x3 OuterProduct(const btScalar* v1,const btScalar* v2,con
     return (m);
 }
 
+static inline btMatrix3x3 OuterProduct(const btVector3& v1,const btVector3& v2)
+{
+    btMatrix3x3 m;
+    btScalar a11 = v1[0] * v2[0];
+    btScalar a12 = v1[0] * v2[1];
+    btScalar a13 = v1[0] * v2[2];
+    
+    btScalar a21 = v1[1] * v2[0];
+    btScalar a22 = v1[1] * v2[1];
+    btScalar a23 = v1[1] * v2[2];
+    
+    btScalar a31 = v1[2] * v2[0];
+    btScalar a32 = v1[2] * v2[1];
+    btScalar a33 = v1[2] * v2[2];
+    m[0] = btVector3(a11, a12, a13);
+    m[1] = btVector3(a21, a22, a23);
+    m[2] = btVector3(a31, a32, a33);
+    return (m);
+}
+
 
 //
 static inline btMatrix3x3 Add(const btMatrix3x3& a,
@@ -1070,8 +1655,8 @@ struct btSoftColliders
 
 			if (!n.m_battach)
             {
-                // check for collision at x_{n+1}^* as well at x_n
-                if (psb->checkDeformableContact(m_colObj1Wrap, n.m_x, m, c.m_cti, /*predict = */ true) || psb->checkDeformableContact(m_colObj1Wrap, n.m_q, m, c.m_cti, /*predict = */ true))
+				// check for collision at x_{n+1}^*
+				if (psb->checkDeformableContact(m_colObj1Wrap, n.m_q, m, c.m_cti, /*predict = */ true))
                 {
                     const btScalar ima = n.m_im;
                     // todo: collision between multibody and fixed deformable node will be missed.
@@ -1159,7 +1744,6 @@ struct btSoftColliders
             btSoftBody::Node* n0 = f.m_n[0];
             btSoftBody::Node* n1 = f.m_n[1];
             btSoftBody::Node* n2 = f.m_n[2];
-            
             const btScalar m = (n0->m_im > 0 && n1->m_im > 0 && n2->m_im > 0 )? dynmargin : stamargin;
             btSoftBody::DeformableFaceRigidContact c;
             btVector3 contact_point;
@@ -1174,18 +1758,19 @@ struct btSoftColliders
                 if (ms > 0)
                 {
                     // resolve contact at x_n
-                    psb->checkDeformableFaceContact(m_colObj1Wrap, f, contact_point, bary, m, c.m_cti, /*predict = */ false);
+//                    psb->checkDeformableFaceContact(m_colObj1Wrap, f, contact_point, bary, m, c.m_cti, /*predict = */ false);
                     btSoftBody::sCti& cti = c.m_cti;
                     c.m_contactPoint = contact_point;
                     c.m_bary = bary;
                     // todo xuchenhan@: this is assuming mass of all vertices are the same. Need to modify if mass are different for distinct vertices
                     c.m_weights = btScalar(2)/(btScalar(1) + bary.length2()) * bary;
                     c.m_face = &f;
+					// friction is handled by the nodes to prevent sticking
+//                    const btScalar fc = 0;
                     const btScalar fc = psb->m_cfg.kDF * m_colObj1Wrap->getCollisionObject()->getFriction();
                     
                     // the effective inverse mass of the face as in https://graphics.stanford.edu/papers/cloth-sig02/cloth.pdf
                     ima = bary.getX()*c.m_weights.getX() * n0->m_im + bary.getY()*c.m_weights.getY() * n1->m_im + bary.getZ()*c.m_weights.getZ() * n2->m_im;
-                    
                     c.m_c2 = ima;
                     c.m_c3 = fc;
                     c.m_c4 = m_colObj1Wrap->getCollisionObject()->isStaticOrKinematicObject() ? psb->m_cfg.kKHR : psb->m_cfg.kCHR;
@@ -1316,19 +1901,11 @@ struct btSoftColliders
         {
             btSoftBody::Node* node = (btSoftBody::Node*)lnode->data;
             btSoftBody::Face* face = (btSoftBody::Face*)lface->data;
-            
-            btVector3 o = node->m_x;
-            btVector3 p;
-            btScalar d = SIMD_INFINITY;
-            ProjectOrigin(face->m_n[0]->m_x - o,
-                          face->m_n[1]->m_x - o,
-                          face->m_n[2]->m_x - o,
-                          p, d);
-            const btScalar m = mrg + (o - node->m_q).safeNorm() * 2;
-            if (d < (m * m))
+            btVector3 bary;
+            if (proximityTest(face->m_n[0]->m_x, face->m_n[1]->m_x, face->m_n[2]->m_x, node->m_x, face->m_normal, mrg, bary))
             {
                 const btSoftBody::Node* n[] = {face->m_n[0], face->m_n[1], face->m_n[2]};
-                const btVector3 w = BaryCoord(n[0]->m_x, n[1]->m_x, n[2]->m_x, p + o);
+                const btVector3 w = bary;
                 const btScalar ma = node->m_im;
                 btScalar mb = BaryEval(n[0]->m_im, n[1]->m_im, n[2]->m_im, w);
                 if ((n[0]->m_im <= 0) ||
@@ -1341,20 +1918,14 @@ struct btSoftColliders
                 if (ms > 0)
                 {
                     btSoftBody::DeformableFaceNodeContact c;
-                    if (useFaceNormal)
-                        c.m_normal = face->m_normal;
-                    else
-                        c.m_normal = p / -btSqrt(d);
+                    c.m_normal = face->m_normal;
+                    if (!useFaceNormal && c.m_normal.dot(node->m_x - face->m_n[2]->m_x) < 0)
+                        c.m_normal  = -face->m_normal;
                     c.m_margin = mrg;
                     c.m_node = node;
                     c.m_face = face;
                     c.m_bary = w;
-                    // todo xuchenhan@: this is assuming mass of all vertices are the same. Need to modify if mass are different for distinct vertices
-                    c.m_weights = btScalar(2)/(btScalar(1) + w.length2()) * w;
                     c.m_friction = psb[0]->m_cfg.kDF * psb[1]->m_cfg.kDF;
-                    // the effective inverse mass of the face as in https://graphics.stanford.edu/papers/cloth-sig02/cloth.pdf
-                    c.m_imf = c.m_bary[0]*c.m_weights[0] * n[0]->m_im + c.m_bary[1]*c.m_weights[1] * n[1]->m_im + c.m_bary[2]*c.m_weights[2] * n[2]->m_im;
-                    c.m_c0 = btScalar(1)/(ma + c.m_imf);
                     psb[0]->m_faceNodeContacts.push_back(c);
                 }
             }
@@ -1372,69 +1943,152 @@ struct btSoftColliders
         void Process(const btDbvntNode* lface1,
                      const btDbvntNode* lface2)
         {
-            btSoftBody::Face* f = (btSoftBody::Face*)lface1->data;
-            btSoftBody::Face* face = (btSoftBody::Face*)lface2->data;
+            btSoftBody::Face* f1 = (btSoftBody::Face*)lface1->data;
+            btSoftBody::Face* f2 = (btSoftBody::Face*)lface2->data;
+            if (f1 != f2)
+            {
+                Repel(f1, f2);
+                Repel(f2, f1);
+            }
+        }
+        void Repel(btSoftBody::Face* f1, btSoftBody::Face* f2)
+        {
+            //#define REPEL_NEIGHBOR 1
+#ifndef REPEL_NEIGHBOR
             for (int node_id = 0; node_id < 3; ++node_id)
             {
-                btSoftBody::Node* node = f->m_n[node_id];
-                bool skip = false;
+                btSoftBody::Node* node = f1->m_n[node_id];
                 for (int i = 0; i < 3; ++i)
                 {
-                    if (face->m_n[i] == node)
+                    if (f2->m_n[i] == node)
+                        return;
+                }
+            }
+#endif
+            bool skip = false;
+            for (int node_id = 0; node_id < 3; ++node_id)
+            {
+                btSoftBody::Node* node = f1->m_n[node_id];
+#ifdef REPEL_NEIGHBOR
+                for (int i = 0; i < 3; ++i)
+                {
+                    if (f2->m_n[i] == node)
                     {
                         skip = true;
                         break;
                     }
                 }
                 if (skip)
+                {
+                    skip = false;
+                    continue;
+                }
+#endif
+                btSoftBody::Face* face = f2;
+                btVector3 bary;
+                if (!proximityTest(face->m_n[0]->m_x, face->m_n[1]->m_x, face->m_n[2]->m_x, node->m_x, face->m_normal, mrg, bary))
                     continue;
-                btVector3 o = node->m_x;
-                btVector3 p;
-                btScalar d = SIMD_INFINITY;
-                ProjectOrigin(face->m_n[0]->m_x - o,
-                              face->m_n[1]->m_x - o,
-                              face->m_n[2]->m_x - o,
-                              p, d);
-                const btScalar m = mrg + (o - node->m_q).safeNorm() * 2;
-                if (d < (m * m))
+                btSoftBody::DeformableFaceNodeContact c;
+                c.m_normal = face->m_normal;
+                if (!useFaceNormal && c.m_normal.dot(node->m_x - face->m_n[2]->m_x) < 0)
+                    c.m_normal  = -face->m_normal;
+                c.m_margin = mrg;
+                c.m_node = node;
+                c.m_face = face;
+                c.m_bary = bary;
+                c.m_friction = psb[0]->m_cfg.kDF * psb[1]->m_cfg.kDF;
+                psb[0]->m_faceNodeContacts.push_back(c);
+            }
+        }
+        btSoftBody* psb[2];
+        btScalar mrg;
+        bool useFaceNormal;
+    };
+
+    struct CollideCCD : btDbvt::ICollide
+    {
+        void Process(const btDbvtNode* lnode,
+                     const btDbvtNode* lface)
+        {
+            btSoftBody::Node* node = (btSoftBody::Node*)lnode->data;
+            btSoftBody::Face* face = (btSoftBody::Face*)lface->data;
+            btVector3 bary;
+            if (bernsteinCCD(face, node, dt, SAFE_EPSILON, bary))
+            {
+                btSoftBody::DeformableFaceNodeContact c;
+                c.m_normal = face->m_normal;
+                if (!useFaceNormal && c.m_normal.dot(node->m_x - face->m_n[2]->m_x) < 0)
+                    c.m_normal  = -face->m_normal;
+                c.m_node = node;
+                c.m_face = face;
+                c.m_bary = bary;
+                c.m_friction = psb[0]->m_cfg.kDF * psb[1]->m_cfg.kDF;
+                psb[0]->m_faceNodeContacts.push_back(c);
+            }
+        }
+        void Process(const btDbvntNode* lface1,
+                     const btDbvntNode* lface2)
+        {
+            btSoftBody::Face* f1 = (btSoftBody::Face*)lface1->data;
+            btSoftBody::Face* f2 = (btSoftBody::Face*)lface2->data;
+            if (f1 != f2)
+            {
+                Repel(f1, f2);
+                Repel(f2, f1);
+            }
+        }
+        void Repel(btSoftBody::Face* f1, btSoftBody::Face* f2)
+        {
+            //#define REPEL_NEIGHBOR 1
+#ifndef REPEL_NEIGHBOR
+            for (int node_id = 0; node_id < 3; ++node_id)
+            {
+                btSoftBody::Node* node = f1->m_n[node_id];
+                for (int i = 0; i < 3; ++i)
                 {
-                    const btSoftBody::Node* n[] = {face->m_n[0], face->m_n[1], face->m_n[2]};
-                    const btVector3 w = BaryCoord(n[0]->m_x, n[1]->m_x, n[2]->m_x, p + o);
-                    const btScalar ma = node->m_im;
-                    btScalar mb = BaryEval(n[0]->m_im, n[1]->m_im, n[2]->m_im, w);
-                    if ((n[0]->m_im <= 0) ||
-                        (n[1]->m_im <= 0) ||
-                        (n[2]->m_im <= 0))
-                    {
-                        mb = 0;
-                    }
-                    const btScalar ms = ma + mb;
-                    if (ms > 0)
+                    if (f2->m_n[i] == node)
+                        return;
+                }
+            }
+#endif
+            bool skip = false;
+            for (int node_id = 0; node_id < 3; ++node_id)
+            {
+                btSoftBody::Node* node = f1->m_n[node_id];
+#ifdef REPEL_NEIGHBOR
+                for (int i = 0; i < 3; ++i)
+                {
+                    if (f2->m_n[i] == node)
                     {
-                        btSoftBody::DeformableFaceNodeContact c;
-                        if (useFaceNormal)
-                            c.m_normal = face->m_normal;
-                        else
-                            c.m_normal = p / -btSqrt(d);
-                        c.m_margin = mrg;
-                        c.m_node = node;
-                        c.m_face = face;
-                        c.m_bary = w;
-                        // todo xuchenhan@: this is assuming mass of all vertices are the same. Need to modify if mass are different for distinct vertices
-                        c.m_weights = btScalar(2)/(btScalar(1) + w.length2()) * w;
-                        c.m_friction = psb[0]->m_cfg.kDF * psb[1]->m_cfg.kDF;
-                        // the effective inverse mass of the face as in https://graphics.stanford.edu/papers/cloth-sig02/cloth.pdf
-                        c.m_imf = c.m_bary[0]*c.m_weights[0] * n[0]->m_im + c.m_bary[1]*c.m_weights[1] * n[1]->m_im + c.m_bary[2]*c.m_weights[2] * n[2]->m_im;
-                        c.m_c0 = btScalar(1)/(ma + c.m_imf);
-                        psb[0]->m_faceNodeContacts.push_back(c);
+                        skip = true;
+                        break;
                     }
                 }
+                if (skip)
+                {
+                    skip = false;
+                    continue;
+                }
+#endif
+                btSoftBody::Face* face = f2;
+                btVector3 bary;
+				if (bernsteinCCD(face, node, dt, SAFE_EPSILON, bary))
+                {
+                    btSoftBody::DeformableFaceNodeContact c;
+                    c.m_normal = face->m_normal;
+                    if (!useFaceNormal && c.m_normal.dot(node->m_x - face->m_n[2]->m_x) < 0)
+                        c.m_normal  = -face->m_normal;
+                    c.m_node = node;
+                    c.m_face = face;
+                    c.m_bary = bary;
+                    c.m_friction = psb[0]->m_cfg.kDF * psb[1]->m_cfg.kDF;
+                    psb[0]->m_faceNodeContacts.push_back(c);
+                }
             }
         }
         btSoftBody* psb[2];
-        btScalar mrg;
+        btScalar dt, mrg;
         bool useFaceNormal;
     };
 };
-
 #endif  //_BT_SOFT_BODY_INTERNALS_H
diff --git a/thirdparty/bullet/BulletSoftBody/btSoftRigidCollisionAlgorithm.cpp b/thirdparty/bullet/BulletSoftBody/btSoftRigidCollisionAlgorithm.cpp
index 56d8083f22..5b65216e4b 100644
--- a/thirdparty/bullet/BulletSoftBody/btSoftRigidCollisionAlgorithm.cpp
+++ b/thirdparty/bullet/BulletSoftBody/btSoftRigidCollisionAlgorithm.cpp
@@ -48,9 +48,10 @@ btSoftRigidCollisionAlgorithm::~btSoftRigidCollisionAlgorithm()
 }
 
 #include <stdio.h>
-
+#include "LinearMath/btQuickprof.h"
 void btSoftRigidCollisionAlgorithm::processCollision(const btCollisionObjectWrapper* body0Wrap, const btCollisionObjectWrapper* body1Wrap, const btDispatcherInfo& dispatchInfo, btManifoldResult* resultOut)
 {
+	BT_PROFILE("btSoftRigidCollisionAlgorithm::processCollision");
 	(void)dispatchInfo;
 	(void)resultOut;
 	//printf("btSoftRigidCollisionAlgorithm\n");
diff --git a/thirdparty/bullet/BulletSoftBody/poly34.cpp b/thirdparty/bullet/BulletSoftBody/poly34.cpp
new file mode 100644
index 0000000000..819d0c79f7
--- /dev/null
+++ b/thirdparty/bullet/BulletSoftBody/poly34.cpp
@@ -0,0 +1,419 @@
+// poly34.cpp : solution of cubic and quartic equation
+// (c) Khashin S.I. http://math.ivanovo.ac.ru/dalgebra/Khashin/index.html
+// khash2 (at) gmail.com
+// Thanks to Alexandr Rakhmanin <rakhmanin (at) gmail.com>
+// public domain
+//
+#include <math.h>
+
+#include "poly34.h" // solution of cubic and quartic equation
+#define TwoPi 6.28318530717958648
+const btScalar eps = SIMD_EPSILON;
+
+//=============================================================================
+// _root3, root3 from http://prografix.narod.ru
+//=============================================================================
+static SIMD_FORCE_INLINE btScalar _root3(btScalar x)
+{
+    btScalar s = 1.;
+    while (x < 1.) {
+        x *= 8.;
+        s *= 0.5;
+    }
+    while (x > 8.) {
+        x *= 0.125;
+        s *= 2.;
+    }
+    btScalar r = 1.5;
+    r -= 1. / 3. * (r - x / (r * r));
+    r -= 1. / 3. * (r - x / (r * r));
+    r -= 1. / 3. * (r - x / (r * r));
+    r -= 1. / 3. * (r - x / (r * r));
+    r -= 1. / 3. * (r - x / (r * r));
+    r -= 1. / 3. * (r - x / (r * r));
+    return r * s;
+}
+
+btScalar SIMD_FORCE_INLINE root3(btScalar x)
+{
+    if (x > 0)
+        return _root3(x);
+    else if (x < 0)
+        return -_root3(-x);
+    else
+        return 0.;
+}
+
+// x - array of size 2
+// return 2: 2 real roots x[0], x[1]
+// return 0: pair of complex roots: x[0]i*x[1]
+int SolveP2(btScalar* x, btScalar a, btScalar b)
+{ // solve equation x^2 + a*x + b = 0
+    btScalar D = 0.25 * a * a - b;
+    if (D >= 0) {
+        D = sqrt(D);
+        x[0] = -0.5 * a + D;
+        x[1] = -0.5 * a - D;
+        return 2;
+    }
+    x[0] = -0.5 * a;
+    x[1] = sqrt(-D);
+    return 0;
+}
+//---------------------------------------------------------------------------
+// x - array of size 3
+// In case 3 real roots: => x[0], x[1], x[2], return 3
+//         2 real roots: x[0], x[1],          return 2
+//         1 real root : x[0], x[1]  i*x[2], return 1
+int SolveP3(btScalar* x, btScalar a, btScalar b, btScalar c)
+{ // solve cubic equation x^3 + a*x^2 + b*x + c = 0
+    btScalar a2 = a * a;
+    btScalar q = (a2 - 3 * b) / 9;
+    if (q < 0)
+        q = eps;
+    btScalar r = (a * (2 * a2 - 9 * b) + 27 * c) / 54;
+    // equation x^3 + q*x + r = 0
+    btScalar r2 = r * r;
+    btScalar q3 = q * q * q;
+    btScalar A, B;
+    if (r2 <= (q3 + eps)) { //<<-- FIXED!
+        btScalar t = r / sqrt(q3);
+        if (t < -1)
+            t = -1;
+        if (t > 1)
+            t = 1;
+        t = acos(t);
+        a /= 3;
+        q = -2 * sqrt(q);
+        x[0] = q * cos(t / 3) - a;
+        x[1] = q * cos((t + TwoPi) / 3) - a;
+        x[2] = q * cos((t - TwoPi) / 3) - a;
+        return (3);
+    }
+    else {
+        //A =-pow(fabs(r)+sqrt(r2-q3),1./3);
+        A = -root3(fabs(r) + sqrt(r2 - q3));
+        if (r < 0)
+            A = -A;
+        B = (A == 0 ? 0 : q / A);
+        
+        a /= 3;
+        x[0] = (A + B) - a;
+        x[1] = -0.5 * (A + B) - a;
+        x[2] = 0.5 * sqrt(3.) * (A - B);
+        if (fabs(x[2]) < eps) {
+            x[2] = x[1];
+            return (2);
+        }
+        return (1);
+    }
+} // SolveP3(btScalar *x,btScalar a,btScalar b,btScalar c) {
+//---------------------------------------------------------------------------
+// a>=0!
+void CSqrt(btScalar x, btScalar y, btScalar& a, btScalar& b) // returns:  a+i*s = sqrt(x+i*y)
+{
+    btScalar r = sqrt(x * x + y * y);
+    if (y == 0) {
+        r = sqrt(r);
+        if (x >= 0) {
+            a = r;
+            b = 0;
+        }
+        else {
+            a = 0;
+            b = r;
+        }
+    }
+    else { // y != 0
+        a = sqrt(0.5 * (x + r));
+        b = 0.5 * y / a;
+    }
+}
+//---------------------------------------------------------------------------
+int SolveP4Bi(btScalar* x, btScalar b, btScalar d) // solve equation x^4 + b*x^2 + d = 0
+{
+    btScalar D = b * b - 4 * d;
+    if (D >= 0) {
+        btScalar sD = sqrt(D);
+        btScalar x1 = (-b + sD) / 2;
+        btScalar x2 = (-b - sD) / 2; // x2 <= x1
+        if (x2 >= 0) // 0 <= x2 <= x1, 4 real roots
+        {
+            btScalar sx1 = sqrt(x1);
+            btScalar sx2 = sqrt(x2);
+            x[0] = -sx1;
+            x[1] = sx1;
+            x[2] = -sx2;
+            x[3] = sx2;
+            return 4;
+        }
+        if (x1 < 0) // x2 <= x1 < 0, two pair of imaginary roots
+        {
+            btScalar sx1 = sqrt(-x1);
+            btScalar sx2 = sqrt(-x2);
+            x[0] = 0;
+            x[1] = sx1;
+            x[2] = 0;
+            x[3] = sx2;
+            return 0;
+        }
+        // now x2 < 0 <= x1 , two real roots and one pair of imginary root
+        btScalar sx1 = sqrt(x1);
+        btScalar sx2 = sqrt(-x2);
+        x[0] = -sx1;
+        x[1] = sx1;
+        x[2] = 0;
+        x[3] = sx2;
+        return 2;
+    }
+    else { // if( D < 0 ), two pair of compex roots
+        btScalar sD2 = 0.5 * sqrt(-D);
+        CSqrt(-0.5 * b, sD2, x[0], x[1]);
+        CSqrt(-0.5 * b, -sD2, x[2], x[3]);
+        return 0;
+    } // if( D>=0 )
+} // SolveP4Bi(btScalar *x, btScalar b, btScalar d)    // solve equation x^4 + b*x^2 d
+//---------------------------------------------------------------------------
+#define SWAP(a, b) \
+{              \
+t = b;     \
+b = a;     \
+a = t;     \
+}
+static void dblSort3(btScalar& a, btScalar& b, btScalar& c) // make: a <= b <= c
+{
+    btScalar t;
+    if (a > b)
+        SWAP(a, b); // now a<=b
+    if (c < b) {
+        SWAP(b, c); // now a<=b, b<=c
+        if (a > b)
+            SWAP(a, b); // now a<=b
+    }
+}
+//---------------------------------------------------------------------------
+int SolveP4De(btScalar* x, btScalar b, btScalar c, btScalar d) // solve equation x^4 + b*x^2 + c*x + d
+{
+    //if( c==0 ) return SolveP4Bi(x,b,d); // After that, c!=0
+    if (fabs(c) < 1e-14 * (fabs(b) + fabs(d)))
+        return SolveP4Bi(x, b, d); // After that, c!=0
+    
+    int res3 = SolveP3(x, 2 * b, b * b - 4 * d, -c * c); // solve resolvent
+    // by Viet theorem:  x1*x2*x3=-c*c not equals to 0, so x1!=0, x2!=0, x3!=0
+    if (res3 > 1) // 3 real roots,
+    {
+        dblSort3(x[0], x[1], x[2]); // sort roots to x[0] <= x[1] <= x[2]
+        // Note: x[0]*x[1]*x[2]= c*c > 0
+        if (x[0] > 0) // all roots are positive
+        {
+            btScalar sz1 = sqrt(x[0]);
+            btScalar sz2 = sqrt(x[1]);
+            btScalar sz3 = sqrt(x[2]);
+            // Note: sz1*sz2*sz3= -c (and not equal to 0)
+            if (c > 0) {
+                x[0] = (-sz1 - sz2 - sz3) / 2;
+                x[1] = (-sz1 + sz2 + sz3) / 2;
+                x[2] = (+sz1 - sz2 + sz3) / 2;
+                x[3] = (+sz1 + sz2 - sz3) / 2;
+                return 4;
+            }
+            // now: c<0
+            x[0] = (-sz1 - sz2 + sz3) / 2;
+            x[1] = (-sz1 + sz2 - sz3) / 2;
+            x[2] = (+sz1 - sz2 - sz3) / 2;
+            x[3] = (+sz1 + sz2 + sz3) / 2;
+            return 4;
+        } // if( x[0] > 0) // all roots are positive
+        // now x[0] <= x[1] < 0, x[2] > 0
+        // two pair of comlex roots
+        btScalar sz1 = sqrt(-x[0]);
+        btScalar sz2 = sqrt(-x[1]);
+        btScalar sz3 = sqrt(x[2]);
+        
+        if (c > 0) // sign = -1
+        {
+            x[0] = -sz3 / 2;
+            x[1] = (sz1 - sz2) / 2; // x[0]i*x[1]
+            x[2] = sz3 / 2;
+            x[3] = (-sz1 - sz2) / 2; // x[2]i*x[3]
+            return 0;
+        }
+        // now: c<0 , sign = +1
+        x[0] = sz3 / 2;
+        x[1] = (-sz1 + sz2) / 2;
+        x[2] = -sz3 / 2;
+        x[3] = (sz1 + sz2) / 2;
+        return 0;
+    } // if( res3>1 )    // 3 real roots,
+    // now resoventa have 1 real and pair of compex roots
+    // x[0] - real root, and x[0]>0,
+    // x[1]i*x[2] - complex roots,
+    // x[0] must be >=0. But one times x[0]=~ 1e-17, so:
+    if (x[0] < 0)
+        x[0] = 0;
+    btScalar sz1 = sqrt(x[0]);
+    btScalar szr, szi;
+    CSqrt(x[1], x[2], szr, szi); // (szr+i*szi)^2 = x[1]+i*x[2]
+    if (c > 0) // sign = -1
+    {
+        x[0] = -sz1 / 2 - szr; // 1st real root
+        x[1] = -sz1 / 2 + szr; // 2nd real root
+        x[2] = sz1 / 2;
+        x[3] = szi;
+        return 2;
+    }
+    // now: c<0 , sign = +1
+    x[0] = sz1 / 2 - szr; // 1st real root
+    x[1] = sz1 / 2 + szr; // 2nd real root
+    x[2] = -sz1 / 2;
+    x[3] = szi;
+    return 2;
+} // SolveP4De(btScalar *x, btScalar b, btScalar c, btScalar d)    // solve equation x^4 + b*x^2 + c*x + d
+//-----------------------------------------------------------------------------
+btScalar N4Step(btScalar x, btScalar a, btScalar b, btScalar c, btScalar d) // one Newton step for x^4 + a*x^3 + b*x^2 + c*x + d
+{
+    btScalar fxs = ((4 * x + 3 * a) * x + 2 * b) * x + c; // f'(x)
+    if (fxs == 0)
+        return x; //return 1e99; <<-- FIXED!
+    btScalar fx = (((x + a) * x + b) * x + c) * x + d; // f(x)
+    return x - fx / fxs;
+}
+//-----------------------------------------------------------------------------
+// x - array of size 4
+// return 4: 4 real roots x[0], x[1], x[2], x[3], possible multiple roots
+// return 2: 2 real roots x[0], x[1] and complex x[2]i*x[3],
+// return 0: two pair of complex roots: x[0]i*x[1],  x[2]i*x[3],
+int SolveP4(btScalar* x, btScalar a, btScalar b, btScalar c, btScalar d)
+{ // solve equation x^4 + a*x^3 + b*x^2 + c*x + d by Dekart-Euler method
+    // move to a=0:
+    btScalar d1 = d + 0.25 * a * (0.25 * b * a - 3. / 64 * a * a * a - c);
+    btScalar c1 = c + 0.5 * a * (0.25 * a * a - b);
+    btScalar b1 = b - 0.375 * a * a;
+    int res = SolveP4De(x, b1, c1, d1);
+    if (res == 4) {
+        x[0] -= a / 4;
+        x[1] -= a / 4;
+        x[2] -= a / 4;
+        x[3] -= a / 4;
+    }
+    else if (res == 2) {
+        x[0] -= a / 4;
+        x[1] -= a / 4;
+        x[2] -= a / 4;
+    }
+    else {
+        x[0] -= a / 4;
+        x[2] -= a / 4;
+    }
+    // one Newton step for each real root:
+    if (res > 0) {
+        x[0] = N4Step(x[0], a, b, c, d);
+        x[1] = N4Step(x[1], a, b, c, d);
+    }
+    if (res > 2) {
+        x[2] = N4Step(x[2], a, b, c, d);
+        x[3] = N4Step(x[3], a, b, c, d);
+    }
+    return res;
+}
+//-----------------------------------------------------------------------------
+#define F5(t) (((((t + a) * t + b) * t + c) * t + d) * t + e)
+//-----------------------------------------------------------------------------
+btScalar SolveP5_1(btScalar a, btScalar b, btScalar c, btScalar d, btScalar e) // return real root of x^5 + a*x^4 + b*x^3 + c*x^2 + d*x + e = 0
+{
+    int cnt;
+    if (fabs(e) < eps)
+        return 0;
+    
+    btScalar brd = fabs(a); // brd - border of real roots
+    if (fabs(b) > brd)
+        brd = fabs(b);
+    if (fabs(c) > brd)
+        brd = fabs(c);
+    if (fabs(d) > brd)
+        brd = fabs(d);
+    if (fabs(e) > brd)
+        brd = fabs(e);
+    brd++; // brd - border of real roots
+    
+    btScalar x0, f0; // less than root
+    btScalar x1, f1; // greater than root
+    btScalar x2, f2, f2s; // next values, f(x2), f'(x2)
+    btScalar dx = 0;
+    
+    if (e < 0) {
+        x0 = 0;
+        x1 = brd;
+        f0 = e;
+        f1 = F5(x1);
+        x2 = 0.01 * brd;
+    } // positive root
+    else {
+        x0 = -brd;
+        x1 = 0;
+        f0 = F5(x0);
+        f1 = e;
+        x2 = -0.01 * brd;
+    } // negative root
+    
+    if (fabs(f0) < eps)
+        return x0;
+    if (fabs(f1) < eps)
+        return x1;
+    
+    // now x0<x1, f(x0)<0, f(x1)>0
+    // Firstly 10 bisections
+    for (cnt = 0; cnt < 10; cnt++) {
+        x2 = (x0 + x1) / 2; // next point
+        //x2 = x0 - f0*(x1 - x0) / (f1 - f0);        // next point
+        f2 = F5(x2); // f(x2)
+        if (fabs(f2) < eps)
+            return x2;
+        if (f2 > 0) {
+            x1 = x2;
+            f1 = f2;
+        }
+        else {
+            x0 = x2;
+            f0 = f2;
+        }
+    }
+    
+    // At each step:
+    // x0<x1, f(x0)<0, f(x1)>0.
+    // x2 - next value
+    // we hope that x0 < x2 < x1, but not necessarily
+    do {
+        if (cnt++ > 50)
+            break;
+        if (x2 <= x0 || x2 >= x1)
+            x2 = (x0 + x1) / 2; // now  x0 < x2 < x1
+        f2 = F5(x2); // f(x2)
+        if (fabs(f2) < eps)
+            return x2;
+        if (f2 > 0) {
+            x1 = x2;
+            f1 = f2;
+        }
+        else {
+            x0 = x2;
+            f0 = f2;
+        }
+        f2s = (((5 * x2 + 4 * a) * x2 + 3 * b) * x2 + 2 * c) * x2 + d; // f'(x2)
+        if (fabs(f2s) < eps) {
+            x2 = 1e99;
+            continue;
+        }
+        dx = f2 / f2s;
+        x2 -= dx;
+    } while (fabs(dx) > eps);
+    return x2;
+} // SolveP5_1(btScalar a,btScalar b,btScalar c,btScalar d,btScalar e)    // return real root of x^5 + a*x^4 + b*x^3 + c*x^2 + d*x + e = 0
+//-----------------------------------------------------------------------------
+int SolveP5(btScalar* x, btScalar a, btScalar b, btScalar c, btScalar d, btScalar e) // solve equation x^5 + a*x^4 + b*x^3 + c*x^2 + d*x + e = 0
+{
+    btScalar r = x[0] = SolveP5_1(a, b, c, d, e);
+    btScalar a1 = a + r, b1 = b + r * a1, c1 = c + r * b1, d1 = d + r * c1;
+    return 1 + SolveP4(x + 1, a1, b1, c1, d1);
+} // SolveP5(btScalar *x,btScalar a,btScalar b,btScalar c,btScalar d,btScalar e)    // solve equation x^5 + a*x^4 + b*x^3 + c*x^2 + d*x + e = 0
+//-----------------------------------------------------------------------------
diff --git a/thirdparty/bullet/BulletSoftBody/poly34.h b/thirdparty/bullet/BulletSoftBody/poly34.h
new file mode 100644
index 0000000000..32ad5d7da5
--- /dev/null
+++ b/thirdparty/bullet/BulletSoftBody/poly34.h
@@ -0,0 +1,38 @@
+// poly34.h : solution of cubic and quartic equation
+// (c) Khashin S.I. http://math.ivanovo.ac.ru/dalgebra/Khashin/index.html
+// khash2 (at) gmail.com
+
+#ifndef POLY_34
+#define POLY_34
+#include "LinearMath/btScalar.h"
+// x - array of size 2
+// return 2: 2 real roots x[0], x[1]
+// return 0: pair of complex roots: x[0]i*x[1]
+int SolveP2(btScalar* x, btScalar a, btScalar b); // solve equation x^2 + a*x + b = 0
+
+// x - array of size 3
+// return 3: 3 real roots x[0], x[1], x[2]
+// return 1: 1 real root x[0] and pair of complex roots: x[1]i*x[2]
+int SolveP3(btScalar* x, btScalar a, btScalar b, btScalar c); // solve cubic equation x^3 + a*x^2 + b*x + c = 0
+
+// x - array of size 4
+// return 4: 4 real roots x[0], x[1], x[2], x[3], possible multiple roots
+// return 2: 2 real roots x[0], x[1] and complex x[2]i*x[3],
+// return 0: two pair of complex roots: x[0]i*x[1],  x[2]i*x[3],
+int SolveP4(btScalar* x, btScalar a, btScalar b, btScalar c, btScalar d); // solve equation x^4 + a*x^3 + b*x^2 + c*x + d = 0 by Dekart-Euler method
+
+// x - array of size 5
+// return 5: 5 real roots x[0], x[1], x[2], x[3], x[4], possible multiple roots
+// return 3: 3 real roots x[0], x[1], x[2] and complex x[3]i*x[4],
+// return 1: 1 real root x[0] and two pair of complex roots: x[1]i*x[2],  x[3]i*x[4],
+int SolveP5(btScalar* x, btScalar a, btScalar b, btScalar c, btScalar d, btScalar e); // solve equation x^5 + a*x^4 + b*x^3 + c*x^2 + d*x + e = 0
+
+//-----------------------------------------------------------------------------
+// And some additional functions for internal use.
+// Your may remove this definitions from here
+int SolveP4Bi(btScalar* x, btScalar b, btScalar d); // solve equation x^4 + b*x^2 + d = 0
+int SolveP4De(btScalar* x, btScalar b, btScalar c, btScalar d); // solve equation x^4 + b*x^2 + c*x + d = 0
+void CSqrt(btScalar x, btScalar y, btScalar& a, btScalar& b); // returns as a+i*s,  sqrt(x+i*y)
+btScalar N4Step(btScalar x, btScalar a, btScalar b, btScalar c, btScalar d); // one Newton step for x^4 + a*x^3 + b*x^2 + c*x + d
+btScalar SolveP5_1(btScalar a, btScalar b, btScalar c, btScalar d, btScalar e); // return real root of x^5 + a*x^4 + b*x^3 + c*x^2 + d*x + e = 0
+#endif