diff --git a/examples/DeformableDemo/GraspDeformable.cpp b/examples/DeformableDemo/GraspDeformable.cpp index 970fde5f4..98fb63a8d 100644 --- a/examples/DeformableDemo/GraspDeformable.cpp +++ b/examples/DeformableDemo/GraspDeformable.cpp @@ -190,6 +190,7 @@ void GraspDeformable::initPhysics() m_guiHelper->createPhysicsDebugDrawer(m_dynamicsWorld); // build a gripper + if(1) { bool damping = true; bool gyro = false; @@ -269,18 +270,18 @@ void GraspDeformable::initPhysics() } // create a soft block - if(0) + if(1) { char relative_path[1024]; // b3FileUtils::findFile("banana.vtk", relative_path, 1024); // b3FileUtils::findFile("ball.vtk", relative_path, 1024); -// b3FileUtils::findFile("paper_collision.vtk", relative_path, 1024); +// b3FileUtils::findFile("deformable_crumpled_napkin_sim.vtk", relative_path, 1024); // b3FileUtils::findFile("single_tet.vtk", relative_path, 1024); - b3FileUtils::findFile("tube.vtk", relative_path, 1024); +// b3FileUtils::findFile("tube.vtk", relative_path, 1024); // b3FileUtils::findFile("torus.vtk", relative_path, 1024); // b3FileUtils::findFile("paper_roll.vtk", relative_path, 1024); // b3FileUtils::findFile("bread.vtk", relative_path, 1024); -// b3FileUtils::findFile("ditto.vtk", relative_path, 1024); + b3FileUtils::findFile("ditto.vtk", relative_path, 1024); // b3FileUtils::findFile("boot.vtk", relative_path, 1024); // btSoftBody* psb = btSoftBodyHelpers::CreateFromTetGenData(getDeformableDynamicsWorld()->getWorldInfo(), // TetraCube::getElements(), @@ -290,12 +291,13 @@ void GraspDeformable::initPhysics() btSoftBody* psb = btSoftBodyHelpers::CreateFromVtkFile(getDeformableDynamicsWorld()->getWorldInfo(), relative_path); // psb->scale(btVector3(30, 30, 30)); // for banana -// psb->scale(btVector3(.2, .2, .2)); + psb->scale(btVector3(.7, .7, .7)); // psb->scale(btVector3(2, 2, 2)); - psb->scale(btVector3(.1, .1, .1)); // for tube, torus, boot -// psb->scale(btVector3(.1, .1, .1)); // for ditto -// psb->translate(btVector3(.25, 2, 0.4)); - psb->getCollisionShape()->setMargin(0.01); +// psb->scale(btVector3(.3, .3, .3)); // for tube, torus, boot + psb->scale(btVector3(.1, .1, .1)); // for ditto +// psb->translate(btVector3(.25, 10, 0.4)); + psb->getCollisionShape()->setMargin(0.0005); + psb->setMaxStress(50); psb->setTotalMass(.01); psb->m_cfg.kKHR = 1; // collision hardness with kinematic objects psb->m_cfg.kCHR = 1; // collision hardness with rigid body @@ -307,20 +309,22 @@ void GraspDeformable::initPhysics() getDeformableDynamicsWorld()->addForce(psb, gravity_force); m_forces.push_back(gravity_force); - btDeformableNeoHookeanForce* neohookean = new btDeformableNeoHookeanForce(100,200, 0.05); + btDeformableNeoHookeanForce* neohookean = new btDeformableNeoHookeanForce(8,32, .05); getDeformableDynamicsWorld()->addForce(psb, neohookean); m_forces.push_back(neohookean); } getDeformableDynamicsWorld()->setImplicit(false); // create a piece of cloth + if(0) { bool onGround = false; - const btScalar s = .4; - btSoftBody* psb = btSoftBodyHelpers::CreatePatch(getDeformableDynamicsWorld()->getWorldInfo(), btVector3(-s, 0, -s), - btVector3(+s, 0, -s), - btVector3(-s, 0, +s), - btVector3(+s, 0, +s), + const btScalar s = .1; + const btScalar h = 1; + btSoftBody* psb = btSoftBodyHelpers::CreatePatch(getDeformableDynamicsWorld()->getWorldInfo(), btVector3(-s, h, -s), + btVector3(+s, h, -s), + btVector3(-s, h, +s), + btVector3(+s, h, +s), 20,20, 0, true); @@ -344,7 +348,8 @@ void GraspDeformable::initPhysics() psb->m_cfg.collisions = btSoftBody::fCollision::SDF_RD; psb->m_cfg.collisions |= btSoftBody::fCollision::VF_DD; getDeformableDynamicsWorld()->addSoftBody(psb); - getDeformableDynamicsWorld()->addForce(psb, new btDeformableMassSpringForce(.2,0.02, true)); +// getDeformableDynamicsWorld()->addForce(psb, new btDeformableMassSpringForce(.0,0.0, true)); + getDeformableDynamicsWorld()->addForce(psb, new btDeformableMassSpringForce(1,0.05, false)); getDeformableDynamicsWorld()->addForce(psb, new btDeformableGravityForce(gravity)); } diff --git a/src/BulletSoftBody/btDeformableNeoHookeanForce.h b/src/BulletSoftBody/btDeformableNeoHookeanForce.h index 589abf2e3..26f738fda 100644 --- a/src/BulletSoftBody/btDeformableNeoHookeanForce.h +++ b/src/BulletSoftBody/btDeformableNeoHookeanForce.h @@ -33,7 +33,7 @@ public: m_lambda_damp = damping * m_lambda; } - btDeformableNeoHookeanForce(btScalar mu, btScalar lambda, btScalar damping = 0): m_mu(mu), m_lambda(lambda) + 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; @@ -73,12 +73,10 @@ public: size_t id2 = node2->index; size_t id3 = node3->index; btMatrix3x3 dF = DsFromVelocity(node0, node1, node2, node3) * tetra.m_Dm_inverse; - btMatrix3x3 C = dF + dF.transpose(); - btMatrix3x3 dP; -// firstPiolaDampingDifferential(psb->m_tetraScratchesTn[j], dF, dP); btMatrix3x3 I; I.setIdentity(); - dP = C * m_mu_damp + I * (dF[0][0]+dF[1][1]+dF[2][2]) * m_lambda_damp; + btMatrix3x3 dP = (dF + dF.transpose()) * m_mu_damp + I * (dF[0][0]+dF[1][1]+dF[2][2]) * m_lambda_damp; +// firstPiolaDampingDifferential(psb->m_tetraScratchesTn[j], dF, dP); btVector3 df_on_node0 = dP * (tetra.m_Dm_inverse.transpose()*grad_N_hat_1st_col); btMatrix3x3 df_on_node123 = dP * tetra.m_Dm_inverse.transpose(); @@ -157,13 +155,17 @@ public: virtual void addScaledElasticForce(btScalar scale, TVStack& force) { -// btScalar max_p = 0; int numNodes = getNumNodes(); btAssert(numNodes <= force.size()); btVector3 grad_N_hat_1st_col = btVector3(-1,-1,-1); for (int i = 0; i < m_softBodies.size(); ++i) { btSoftBody* psb = m_softBodies[i]; + if (!psb->isActive()) + { + continue; + } + btScalar max_p = psb->m_cfg.m_maxStress; for (int j = 0; j < psb->m_tetras.size(); ++j) { btSoftBody::Tetra& tetra = psb->m_tetras[j]; @@ -173,47 +175,21 @@ public: btMatrix3x3 U, V; btVector3 sigma; singularValueDecomposition(P, U, sigma, V); - sigma[0] = btMin(sigma[0], (btScalar)1000); - sigma[1] = btMin(sigma[1], (btScalar)1000); - sigma[2] = btMin(sigma[2], (btScalar)1000); - sigma[0] = btMax(sigma[0], (btScalar)-1000); - sigma[1] = btMax(sigma[1], (btScalar)-1000); - sigma[2] = btMax(sigma[2], (btScalar)-1000); + if (max_p > 0) + { + sigma[0] = btMin(sigma[0], max_p); + sigma[1] = btMin(sigma[1], max_p); + sigma[2] = btMin(sigma[2], max_p); + sigma[0] = btMax(sigma[0], -max_p); + sigma[1] = btMax(sigma[1], -max_p); + sigma[2] = btMax(sigma[2], -max_p); + } btMatrix3x3 Sigma; Sigma.setIdentity(); Sigma[0][0] = sigma[0]; Sigma[1][1] = sigma[1]; Sigma[2][2] = sigma[2]; -// max_p = btMax(max_p, sigma[0]); -// max_p = btMax(max_p, sigma[1]); -// max_p = btMax(max_p, sigma[2]); P = U * Sigma * V.transpose(); - -// Eigen::Matrix eigenP; -// eigenP << P[0][0], P[0][1], P[0][2], -// P[1][0],P[1][1],P[1][2], -// P[2][0],P[2][1],P[2][2]; -// Eigen::JacobiSVD > svd(eigenP, Eigen::ComputeFullU | Eigen::ComputeFullV); -// Eigen::Matrix P_hat = svd.singularValues().asDiagonal(); -// P_hat(0,0) = btMin((btScalar)P_hat(0,0), (btScalar)20); -// P_hat(1,1) = btMin((btScalar)P_hat(1,1), (btScalar)20); -// P_hat(2,2) = btMin((btScalar)P_hat(2,2), (btScalar)20); -// eigenP = svd.matrixU() * P_hat * svd.matrixV().transpose(); -// max_p = btMax(max_p, (btScalar)P_hat(0,0)); -// max_p = btMax(max_p, (btScalar)P_hat(1,1)); -// max_p = btMax(max_p, (btScalar)P_hat(2,2)); -// -// P[0][0] = eigenP(0,0); -// P[0][1] = eigenP(0,1); -// P[0][2] = eigenP(0,2); -// P[1][0] = eigenP(1,0); -// P[1][1] = eigenP(1,1); -// P[1][2] = eigenP(1,2); -// P[2][0] = eigenP(2,0); -// P[2][1] = eigenP(2,1); -// P[2][2] = eigenP(2,2); - - // 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; @@ -235,7 +211,6 @@ public: force[id3] -= scale1 * force_on_node123.getColumn(2); } } -// std::cout << max_p << std::endl; } // The damping matrix is calculated using the time n state as described in https://www.math.ucla.edu/~jteran/papers/GSSJT15.pdf to allow line search @@ -261,12 +236,10 @@ public: size_t id2 = node2->index; size_t id3 = node3->index; btMatrix3x3 dF = Ds(id0, id1, id2, id3, dv) * tetra.m_Dm_inverse; - btMatrix3x3 C = dF + dF.transpose(); - btMatrix3x3 dP; -// firstPiolaDampingDifferential(psb->m_tetraScratchesTn[j], dF, dP); btMatrix3x3 I; I.setIdentity(); - dP = C * m_mu_damp + I * (dF[0][0]+dF[1][1]+dF[2][2]) * m_lambda_damp; + btMatrix3x3 dP = (dF + dF.transpose()) * m_mu_damp + I * (dF[0][0]+dF[1][1]+dF[2][2]) * m_lambda_damp; +// firstPiolaDampingDifferential(psb->m_tetraScratchesTn[j], dF, dP); // btVector3 df_on_node0 = dP * (tetra.m_Dm_inverse.transpose()*grad_N_hat_1st_col); btMatrix3x3 df_on_node123 = dP * tetra.m_Dm_inverse.transpose(); btVector3 df_on_node0 = df_on_node123 * grad_N_hat_1st_col; diff --git a/src/BulletSoftBody/btSoftBody.cpp b/src/BulletSoftBody/btSoftBody.cpp index c3d99590a..253de645d 100644 --- a/src/BulletSoftBody/btSoftBody.cpp +++ b/src/BulletSoftBody/btSoftBody.cpp @@ -90,6 +90,7 @@ void btSoftBody::initDefaults() m_cfg.diterations = 0; m_cfg.citerations = 4; m_cfg.drag = 0; + m_cfg.m_maxStress = 0; m_cfg.collisions = fCollision::Default; m_cfg.collisions |= fCollision::VF_DD; m_pose.m_bvolume = false; @@ -3195,6 +3196,12 @@ void btSoftBody::applyForces() } } +// +void btSoftBody::setMaxStress(btScalar maxStress) +{ + m_cfg.m_maxStress = maxStress; +} + // void btSoftBody::interpolateRenderMesh() { diff --git a/src/BulletSoftBody/btSoftBody.h b/src/BulletSoftBody/btSoftBody.h index 31ddfabb9..4f4f6f22b 100644 --- a/src/BulletSoftBody/btSoftBody.h +++ b/src/BulletSoftBody/btSoftBody.h @@ -703,6 +703,7 @@ public: tPSolverArray m_psequence; // Position solvers sequence tPSolverArray m_dsequence; // Drift solvers sequence btScalar drag; // deformable air drag + btScalar m_maxStress; // Maximum principle first Piola stress }; /* SolverState */ struct SolverState @@ -1107,6 +1108,7 @@ public: void updateDeformation(); void advanceDeformation(); void applyForces(); + void setMaxStress(btScalar maxStress); void interpolateRenderMesh(); static void PSolve_Anchors(btSoftBody* psb, btScalar kst, btScalar ti); static void PSolve_RContacts(btSoftBody* psb, btScalar kst, btScalar ti);