AuroraMiddleware/bullet3
1
0
Fork 0
mirror of https://github.com/bulletphysics/bullet3 synced 2025-01-19 05:20:06 +00:00
Code Issues Releases Wiki Activity

add mass damping into damping model

Browse Source
This commit is contained in:
Xuchen Han 2020-06-18 16:43:53 -07:00
parent abaf278c2d
commit 92f2f56b3d
6 changed files with 71 additions and 33 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

14
examples/DeformableDemo/LargeDeformation.cpp
View File

@ -29,7 +29,8 @@
///The LargeDeformation shows the contact between volumetric deformable objects and rigid objects.
static btScalar E = 25;
static btScalar nu = 0.3;
static btScalar damping = 0.01;
static btScalar damping_alpha = 0.01;
static btScalar damping_beta = 0.01;
struct TetraCube
{
@ -66,7 +67,7 @@ public:
{
m_linearElasticity->setPoissonRatio(nu);
m_linearElasticity->setYoungsModulus(E);
m_linearElasticity->setDamping(damping);
m_linearElasticity->setDamping(damping_alpha, damping_beta);
float internalTimeStep = 1. / 60.f;
m_dynamicsWorld->stepSimulation(deltaTime, 1, internalTimeStep);
}
@ -164,12 +165,19 @@ void LargeDeformation::initPhysics()
m_guiHelper->getParameterInterface()->registerSliderFloatParameter(slider);
}
{
SliderParams slider("Damping", &damping);
SliderParams slider("Mass Damping", &damping_alpha);
slider.m_minVal = 0.001;
slider.m_maxVal = 0.01;
if (m_guiHelper->getParameterInterface())
m_guiHelper->getParameterInterface()->registerSliderFloatParameter(slider);
}
{
SliderParams slider("Stiffness Damping", &damping_beta);
slider.m_minVal = 0.001;
slider.m_maxVal = 0.01;
if (m_guiHelper->getParameterInterface())
m_guiHelper->getParameterInterface()->registerSliderFloatParameter(slider);
}
}
void LargeDeformation::exitPhysics()

26
examples/DeformableDemo/VolumetricDeformable.cpp
View File

@ -27,9 +27,10 @@
#include "../Utils/b3ResourcePath.h"
///The VolumetricDeformable shows the contact between volumetric deformable objects and rigid objects.
static btScalar E = 25;
static btScalar E = 50;
static btScalar nu = 0.3;
static btScalar damping = 0.01;
static btScalar damping_alpha = 0.1;
static btScalar damping_beta = 0.01;
struct TetraCube
{
@ -66,9 +67,9 @@ public:
{
m_linearElasticity->setPoissonRatio(nu);
m_linearElasticity->setYoungsModulus(E);
m_linearElasticity->setDamping(damping);
m_linearElasticity->setDamping(damping_alpha, damping_beta);
//use a smaller internal timestep, there are stability issues
float internalTimeStep = 1. / 240.f;
float internalTimeStep = 1. / 240;
m_dynamicsWorld->stepSimulation(deltaTime, 4, internalTimeStep);
}
@ -187,7 +188,7 @@ void VolumetricDeformable::initPhysics()
btDefaultMotionState* myMotionState = new btDefaultMotionState(groundTransform);
btRigidBody::btRigidBodyConstructionInfo rbInfo(mass, myMotionState, groundShape, localInertia);
btRigidBody* body = new btRigidBody(rbInfo);
body->setFriction(0.5);
body->setFriction(1);
//add the ground to the dynamics world
m_dynamicsWorld->addRigidBody(body);
@ -212,7 +213,7 @@ void VolumetricDeformable::initPhysics()
psb->setTotalMass(0.5);
psb->m_cfg.kKHR = 1; // collision hardness with kinematic objects
psb->m_cfg.kCHR = 1; // collision hardness with rigid body
psb->m_cfg.kDF = 0.5;
psb->m_cfg.kDF = 2;
psb->m_cfg.collisions = btSoftBody::fCollision::SDF_RD;
psb->m_cfg.collisions |= btSoftBody::fCollision::SDF_RDN;
psb->m_sleepingThreshold = 0;
@ -253,12 +254,19 @@ void VolumetricDeformable::initPhysics()
m_guiHelper->getParameterInterface()->registerSliderFloatParameter(slider);
}
{
SliderParams slider("Damping", &damping);
slider.m_minVal = 0.001;
slider.m_maxVal = 0.01;
SliderParams slider("Mass Damping", &damping_alpha);
slider.m_minVal = 0;
slider.m_maxVal = 1;
if (m_guiHelper->getParameterInterface())
m_guiHelper->getParameterInterface()->registerSliderFloatParameter(slider);
}
{
SliderParams slider("Stiffness Damping", &damping_beta);
slider.m_minVal = 0;
slider.m_maxVal = 0.1;
if (m_guiHelper->getParameterInterface())
m_guiHelper->getParameterInterface()->registerSliderFloatParameter(slider);
}
}
void VolumetricDeformable::exitPhysics()

2
src/BulletSoftBody/btConjugateGradient.h
View File

@ -33,6 +33,7 @@ public:
// return the number of iterations taken
int solve(MatrixX& A, TVStack& x, const TVStack& b, bool verbose = false)
{
verbose = true;
BT_PROFILE("CGSolve");
btAssert(x.size() == b.size());
reinitialize(b);
@ -41,6 +42,7 @@ public:
p = temp;
A.precondition(p,z);
btScalar d0 = this->dot(z,temp);
d0 = btMin(btScalar(1), d0);
// r = b - A * x --with assigned dof zeroed out
A.multiply(x, temp);
r = this->sub(b, temp);

53
src/BulletSoftBody/btDeformableLinearElasticityForce.h
View File

@ -18,26 +18,21 @@
#include "btDeformableLagrangianForce.h"
#include "LinearMath/btQuickprof.h"
#define TETRA_FLAT_THRESHOLD 0.6
#define TETRA_FLAT_THRESHOLD 0.01
class btDeformableLinearElasticityForce : public btDeformableLagrangianForce
{
public:
typedef btAlignedObjectArray<btVector3> TVStack;
btScalar m_mu, m_lambda;
btScalar m_E, m_nu; // Young's modulus and Poisson ratio
btScalar m_mu_damp, m_lambda_damp;
btDeformableLinearElasticityForce(): m_mu(1), m_lambda(1)
btScalar m_damping_alpha, m_damping_beta;
btDeformableLinearElasticityForce(): m_mu(1), m_lambda(1), m_damping_alpha(0.01), m_damping_beta(0.01)
{
btScalar damping = 0.05;
m_mu_damp = damping * m_mu;
m_lambda_damp = damping * m_lambda;
updateYoungsModulusAndPoissonRatio();
}
btDeformableLinearElasticityForce(btScalar mu, btScalar lambda, btScalar damping = 0.05): m_mu(mu), m_lambda(lambda)
btDeformableLinearElasticityForce(btScalar mu, btScalar lambda, btScalar damping_alpha = 0.01, btScalar damping_beta = 0.01): m_mu(mu), m_lambda(lambda), m_damping_alpha(damping_alpha), m_damping_beta(damping_beta)
{
m_mu_damp = damping * m_mu;
m_lambda_damp = damping * m_lambda;
updateYoungsModulusAndPoissonRatio();
}
@ -69,10 +64,10 @@ public:
updateLameParameters();
}
void setDamping(btScalar damping)
void setDamping(btScalar damping_alpha, btScalar damping_beta)
{
m_mu_damp = damping * m_mu;
m_lambda_damp = damping * m_lambda;
m_damping_alpha = damping_alpha;
m_damping_beta = damping_beta;
}
void setLameParameters(btScalar mu, btScalar lambda)
@ -96,8 +91,10 @@ public:
// 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
virtual void addScaledDampingForce(btScalar scale, TVStack& force)
{
if (m_mu_damp == 0 && m_lambda_damp == 0)
if (m_damping_alpha == 0 && m_damping_beta == 0)
return;
btScalar mu_damp = m_damping_beta * m_mu;
btScalar lambda_damp = m_damping_beta * m_lambda;
int numNodes = getNumNodes();
btAssert(numNodes <= force.size());
btVector3 grad_N_hat_1st_col = btVector3(-1,-1,-1);
@ -127,7 +124,7 @@ public:
}
btMatrix3x3 I;
I.setIdentity();
btMatrix3x3 dP = (dF + dF.transpose()) * m_mu_damp + I * ((dF[0][0]+dF[1][1]+dF[2][2]) * m_lambda_damp);
btMatrix3x3 dP = (dF + dF.transpose()) * mu_damp + I * ((dF[0][0]+dF[1][1]+dF[2][2]) * lambda_damp);
btMatrix3x3 df_on_node123 = dP * tetra.m_Dm_inverse.transpose();
if (!close_to_flat)
{
@ -141,6 +138,15 @@ public:
force[id2] -= scale1 * df_on_node123.getColumn(1);
force[id3] -= scale1 * df_on_node123.getColumn(2);
}
for (int j = 0; j < psb->m_nodes.size(); ++j)
{
const btSoftBody::Node& node = psb->m_nodes[j];
size_t id = node.index;
if (node.m_im > 0)
{
force[id] -= scale * node.m_v / node.m_im * m_damping_alpha;
}
}
}
}
@ -280,8 +286,10 @@ public:
// 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
virtual void addScaledDampingForceDifferential(btScalar scale, const TVStack& dv, TVStack& df)
{
if (m_mu_damp == 0 && m_lambda_damp == 0)
if (m_damping_alpha == 0 && m_damping_beta == 0)
return;
btScalar mu_damp = m_damping_beta * m_mu;
btScalar lambda_damp = m_damping_beta * m_lambda;
int numNodes = getNumNodes();
btAssert(numNodes <= df.size());
btVector3 grad_N_hat_1st_col = btVector3(-1,-1,-1);
@ -311,7 +319,7 @@ public:
}
btMatrix3x3 I;
I.setIdentity();
btMatrix3x3 dP = (dF + dF.transpose()) * m_mu_damp + I * ((dF[0][0]+dF[1][1]+dF[2][2]) * m_lambda_damp);
btMatrix3x3 dP = (dF + dF.transpose()) * mu_damp + I * ((dF[0][0]+dF[1][1]+dF[2][2]) * lambda_damp);
btMatrix3x3 df_on_node123 = dP * tetra.m_Dm_inverse.transpose();
if (!close_to_flat)
{
@ -326,6 +334,15 @@ public:
df[id2] -= scale1 * df_on_node123.getColumn(1);
df[id3] -= scale1 * df_on_node123.getColumn(2);
}
for (int j = 0; j < psb->m_nodes.size(); ++j)
{
const btSoftBody::Node& node = psb->m_nodes[j];
size_t id = node.index;
if (node.m_im > 0)
{
df[id] -= scale * dv[id] / node.m_im * m_damping_alpha;
}
}
}
}
@ -390,8 +407,10 @@ public:
// This function calculates the dP = dQ/dF * dF
void firstPiolaDampingDifferential(const btSoftBody::TetraScratch& s, const btMatrix3x3& dF, btMatrix3x3& dP)
{
btScalar mu_damp = m_damping_beta * m_mu;
btScalar lambda_damp = m_damping_beta * m_lambda;
btScalar trace = (dF[0][0] + dF[1][1] + dF[2][2]);
dP = (dF + dF.transpose()) * m_mu_damp + btMatrix3x3::getIdentity() * m_lambda_damp * trace;
dP = (dF + dF.transpose()) * mu_damp + btMatrix3x3::getIdentity() * lambda_damp * trace;
}
virtual btDeformableLagrangianForceType getForceType()

1
src/BulletSoftBody/btDeformableMultiBodyConstraintSolver.cpp
View File

@ -40,6 +40,7 @@ btScalar btDeformableMultiBodyConstraintSolver::solveDeformableGroupIterations(b
if (m_leastSquaresResidual <= infoGlobal.m_leastSquaresResidualThreshold || (iteration >= (maxIterations - 1)))
{
#define VERBOSE_RESIDUAL_PRINTF
#ifdef VERBOSE_RESIDUAL_PRINTF
if (iteration >= (maxIterations - 1))
printf("residual = %f at iteration #%d\n", m_leastSquaresResidual, iteration);

8
src/BulletSoftBody/btSoftBody.cpp
View File

@ -3525,10 +3525,10 @@ void btSoftBody::updateDeformation()
s.m_trace = C[0].getX() + C[1].getY() + C[2].getZ();
s.m_cofF = t.m_F.adjoint().transpose();
btVector3 a = t.m_n[0]->m_x;
btVector3 b = t.m_n[1]->m_x;
btVector3 c = t.m_n[2]->m_x;
btVector3 d = t.m_n[3]->m_x;
btVector3 a = t.m_n[0]->m_q;
btVector3 b = t.m_n[1]->m_q;
btVector3 c = t.m_n[2]->m_q;
btVector3 d = t.m_n[3]->m_q;
btVector4 q1(a[0], b[0], c[0], d[0]);
btVector4 q2(a[1], b[1], c[1], d[1]);
btVector4 q3(a[2], b[2], c[2], d[2]);