diff --git a/examples/DeformableDemo/DeformableContact.cpp b/examples/DeformableDemo/DeformableContact.cpp index 2432d3e27..8712d6de8 100644 --- a/examples/DeformableDemo/DeformableContact.cpp +++ b/examples/DeformableDemo/DeformableContact.cpp @@ -105,7 +105,8 @@ void DeformableContact::initPhysics() btVector3 gravity = btVector3(0, -10, 0); m_dynamicsWorld->setGravity(gravity); getDeformableDynamicsWorld()->getWorldInfo().m_gravity = gravity; - + getDeformableDynamicsWorld()->getWorldInfo().m_sparsesdf.setDefaultVoxelsz(0.25); + getDeformableDynamicsWorld()->getWorldInfo().m_sparsesdf.Reset(); m_guiHelper->createPhysicsDebugDrawer(m_dynamicsWorld); { @@ -182,7 +183,7 @@ void DeformableContact::initPhysics() psb2->setTotalMass(1); psb2->m_cfg.kKHR = 1; // collision hardness with kinematic objects psb2->m_cfg.kCHR = 1; // collision hardness with rigid body - psb2->m_cfg.kDF = 0; + psb2->m_cfg.kDF = 0.5; psb2->m_cfg.collisions = btSoftBody::fCollision::SDF_RD; psb2->m_cfg.collisions |= btSoftBody::fCollision::VF_DD; psb->translate(btVector3(3.5,0,0)); diff --git a/src/.DS_Store b/src/.DS_Store index ee616d543..26f6e2211 100644 Binary files a/src/.DS_Store and b/src/.DS_Store differ diff --git a/src/BulletSoftBody/btDeformableMultiBodyDynamicsWorld.cpp b/src/BulletSoftBody/btDeformableMultiBodyDynamicsWorld.cpp index 9f8e13e20..03a884139 100644 --- a/src/BulletSoftBody/btDeformableMultiBodyDynamicsWorld.cpp +++ b/src/BulletSoftBody/btDeformableMultiBodyDynamicsWorld.cpp @@ -17,14 +17,19 @@ /* A single step of the deformable body simulation contains the following main components: -1. Update velocity to a temporary state v_{n+1}^* = v_n + explicit_force * dt / mass, where explicit forces include gravity and elastic forces. -2. Detect collisions between rigid and deformable bodies at position x_{n+1}^* = x_n + dt * v_{n+1}^*. -3. Then velocities of deformable bodies v_{n+1} are solved in +Call internalStepSimulation multiple times, to achieve 240Hz (4 steps of 60Hz). +1. Deformable maintaintenance of rest lengths and volume preservation. Forces only depend on position: Update velocity to a temporary state v_{n+1}^* = v_n + explicit_force * dt / mass, where explicit forces include gravity and elastic forces. +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, by a conjugate gradient solver, where the damping force is implicit and depends on v_{n+1}. -4. Contact constraints are solved as projections as in the paper by Baraff and Witkin https://www.cs.cmu.edu/~baraff/papers/sig98.pdf. Dynamic frictions are treated as a force and added to the rhs of the CG solve, whereas static frictions are treated as constraints similar to contact. -5. Position is updated via x_{n+1} = x_n + dt * v_{n+1}. -6. Apply position correction to prevent numerical drift. + Make sure contact constraints are not violated in step b by performing velocity projections as in the paper by Baraff and Witkin https://www.cs.cmu.edu/~baraff/papers/sig98.pdf. Dynamic frictions are treated as a force and added to the rhs of the CG solve, whereas static frictions are treated as constraints similar to contact. +4. Position is updated via x_{n+1} = x_n + dt * v_{n+1}. + The algorithm also closely resembles the one in http://physbam.stanford.edu/~fedkiw/papers/stanford2008-03.pdf */ diff --git a/src/BulletSoftBody/btDeformableNeoHookeanForce.h b/src/BulletSoftBody/btDeformableNeoHookeanForce.h index b113b4546..cbe959bdc 100644 --- a/src/BulletSoftBody/btDeformableNeoHookeanForce.h +++ b/src/BulletSoftBody/btDeformableNeoHookeanForce.h @@ -175,7 +175,7 @@ public: btSoftBody::Tetra& tetra = psb->m_tetras[j]; btMatrix3x3 P; firstPiola(psb->m_tetraScratches[j],P); - +#if USE_SVD btMatrix3x3 U, V; btVector3 sigma; singularValueDecomposition(P, U, sigma, V); @@ -194,6 +194,7 @@ public: Sigma[1][1] = sigma[1]; Sigma[2][2] = sigma[2]; P = U * Sigma * V.transpose(); +#endif // btVector3 force_on_node0 = P * (tetra.m_Dm_inverse.transpose()*grad_N_hat_1st_col); btMatrix3x3 force_on_node123 = P * tetra.m_Dm_inverse.transpose(); btVector3 force_on_node0 = force_on_node123 * grad_N_hat_1st_col; diff --git a/src/BulletSoftBody/btSoftBodyInternals.h b/src/BulletSoftBody/btSoftBodyInternals.h index 6e20b3222..93aa1d3c0 100644 --- a/src/BulletSoftBody/btSoftBodyInternals.h +++ b/src/BulletSoftBody/btSoftBodyInternals.h @@ -1288,7 +1288,7 @@ struct btSoftColliders c.m_node = node; c.m_face = face; c.m_weights = w; - c.m_friction = btMax(psb[0]->m_cfg.kDF, psb[1]->m_cfg.kDF); + c.m_friction = btMax (psb[0]->m_cfg.kDF, psb[1]->m_cfg.kDF); c.m_cfm[0] = ma / ms * psb[0]->m_cfg.kSHR; c.m_cfm[1] = mb / ms * psb[1]->m_cfg.kSHR; psb[0]->m_scontacts.push_back(c); @@ -1346,7 +1346,7 @@ struct btSoftColliders 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 = btMax(psb[0]->m_cfg.kDF, psb[1]->m_cfg.kDF); + 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); @@ -1408,7 +1408,7 @@ struct btSoftColliders 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 = btMax(psb[0]->m_cfg.kDF, psb[1]->m_cfg.kDF); + 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); @@ -1462,7 +1462,7 @@ struct btSoftColliders 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 = btMax(psb[0]->m_cfg.kDF, psb[1]->m_cfg.kDF); + 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);