1118 lines
36 KiB
C++
1118 lines
36 KiB
C++
|
/*
|
||
|
|
Bullet Continuous Collision Detection and Physics Library
|
||
|
|
btConeTwistConstraint is Copyright (c) 2007 Starbreeze Studios
|
||
|
|
|
||
|
|
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.
|
||
|
|
|
||
|
|
Written by: Marcus Hennix
|
||
|
|
*/
|
||
|
|
|
||
|
|
|
||
|
|
#include "btConeTwistConstraint.h"
|
||
|
|
#include "BulletDynamics/Dynamics/btRigidBody.h"
|
||
|
|
#include "LinearMath/btTransformUtil.h"
|
||
|
|
#include "LinearMath/btMinMax.h"
|
||
|
|
#include <new>
|
||
|
|
|
||
|
|
|
||
|
|
|
||
|
|
//#define CONETWIST_USE_OBSOLETE_SOLVER true
|
||
|
|
#define CONETWIST_USE_OBSOLETE_SOLVER false
|
||
|
|
#define CONETWIST_DEF_FIX_THRESH btScalar(.05f)
|
||
|
|
|
||
|
|
|
||
|
|
SIMD_FORCE_INLINE btScalar computeAngularImpulseDenominator(const btVector3& axis, const btMatrix3x3& invInertiaWorld)
|
||
|
|
{
|
||
|
|
btVector3 vec = axis * invInertiaWorld;
|
||
|
|
return axis.dot(vec);
|
||
|
|
}
|
||
|
|
|
||
|
|
|
||
|
|
|
||
|
|
|
||
|
|
btConeTwistConstraint::btConeTwistConstraint(btRigidBody& rbA,btRigidBody& rbB,
|
||
|
|
const btTransform& rbAFrame,const btTransform& rbBFrame)
|
||
|
|
:btTypedConstraint(CONETWIST_CONSTRAINT_TYPE, rbA,rbB),m_rbAFrame(rbAFrame),m_rbBFrame(rbBFrame),
|
||
|
|
m_angularOnly(false),
|
||
|
|
m_useSolveConstraintObsolete(CONETWIST_USE_OBSOLETE_SOLVER)
|
||
|
|
{
|
||
|
|
init();
|
||
|
|
}
|
||
|
|
|
||
|
|
btConeTwistConstraint::btConeTwistConstraint(btRigidBody& rbA,const btTransform& rbAFrame)
|
||
|
|
:btTypedConstraint(CONETWIST_CONSTRAINT_TYPE,rbA),m_rbAFrame(rbAFrame),
|
||
|
|
m_angularOnly(false),
|
||
|
|
m_useSolveConstraintObsolete(CONETWIST_USE_OBSOLETE_SOLVER)
|
||
|
|
{
|
||
|
|
m_rbBFrame = m_rbAFrame;
|
||
|
|
init();
|
||
|
|
}
|
||
|
|
|
||
|
|
|
||
|
|
void btConeTwistConstraint::init()
|
||
|
|
{
|
||
|
|
m_angularOnly = false;
|
||
|
|
m_solveTwistLimit = false;
|
||
|
|
m_solveSwingLimit = false;
|
||
|
|
m_bMotorEnabled = false;
|
||
|
|
m_maxMotorImpulse = btScalar(-1);
|
||
|
|
|
||
|
|
setLimit(btScalar(BT_LARGE_FLOAT), btScalar(BT_LARGE_FLOAT), btScalar(BT_LARGE_FLOAT));
|
||
|
|
m_damping = btScalar(0.01);
|
||
|
|
m_fixThresh = CONETWIST_DEF_FIX_THRESH;
|
||
|
|
m_flags = 0;
|
||
|
|
m_linCFM = btScalar(0.f);
|
||
|
|
m_linERP = btScalar(0.7f);
|
||
|
|
m_angCFM = btScalar(0.f);
|
||
|
|
}
|
||
|
|
|
||
|
|
|
||
|
|
void btConeTwistConstraint::getInfo1 (btConstraintInfo1* info)
|
||
|
|
{
|
||
|
|
if (m_useSolveConstraintObsolete)
|
||
|
|
{
|
||
|
|
info->m_numConstraintRows = 0;
|
||
|
|
info->nub = 0;
|
||
|
|
}
|
||
|
|
else
|
||
|
|
{
|
||
|
|
info->m_numConstraintRows = 3;
|
||
|
|
info->nub = 3;
|
||
|
|
calcAngleInfo2(m_rbA.getCenterOfMassTransform(),m_rbB.getCenterOfMassTransform(),m_rbA.getInvInertiaTensorWorld(),m_rbB.getInvInertiaTensorWorld());
|
||
|
|
if(m_solveSwingLimit)
|
||
|
|
{
|
||
|
|
info->m_numConstraintRows++;
|
||
|
|
info->nub--;
|
||
|
|
if((m_swingSpan1 < m_fixThresh) && (m_swingSpan2 < m_fixThresh))
|
||
|
|
{
|
||
|
|
info->m_numConstraintRows++;
|
||
|
|
info->nub--;
|
||
|
|
}
|
||
|
|
}
|
||
|
|
if(m_solveTwistLimit)
|
||
|
|
{
|
||
|
|
info->m_numConstraintRows++;
|
||
|
|
info->nub--;
|
||
|
|
}
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
void btConeTwistConstraint::getInfo1NonVirtual (btConstraintInfo1* info)
|
||
|
|
{
|
||
|
|
//always reserve 6 rows: object transform is not available on SPU
|
||
|
|
info->m_numConstraintRows = 6;
|
||
|
|
info->nub = 0;
|
||
|
|
|
||
|
|
}
|
||
|
|
|
||
|
|
|
||
|
|
void btConeTwistConstraint::getInfo2 (btConstraintInfo2* info)
|
||
|
|
{
|
||
|
|
getInfo2NonVirtual(info,m_rbA.getCenterOfMassTransform(),m_rbB.getCenterOfMassTransform(),m_rbA.getInvInertiaTensorWorld(),m_rbB.getInvInertiaTensorWorld());
|
||
|
|
}
|
||
|
|
|
||
|
|
void btConeTwistConstraint::getInfo2NonVirtual (btConstraintInfo2* info,const btTransform& transA,const btTransform& transB,const btMatrix3x3& invInertiaWorldA,const btMatrix3x3& invInertiaWorldB)
|
||
|
|
{
|
||
|
|
calcAngleInfo2(transA,transB,invInertiaWorldA,invInertiaWorldB);
|
||
|
|
|
||
|
|
btAssert(!m_useSolveConstraintObsolete);
|
||
|
|
// set jacobian
|
||
|
|
info->m_J1linearAxis[0] = 1;
|
||
|
|
info->m_J1linearAxis[info->rowskip+1] = 1;
|
||
|
|
info->m_J1linearAxis[2*info->rowskip+2] = 1;
|
||
|
|
btVector3 a1 = transA.getBasis() * m_rbAFrame.getOrigin();
|
||
|
|
{
|
||
|
|
btVector3* angular0 = (btVector3*)(info->m_J1angularAxis);
|
||
|
|
btVector3* angular1 = (btVector3*)(info->m_J1angularAxis+info->rowskip);
|
||
|
|
btVector3* angular2 = (btVector3*)(info->m_J1angularAxis+2*info->rowskip);
|
||
|
|
btVector3 a1neg = -a1;
|
||
|
|
a1neg.getSkewSymmetricMatrix(angular0,angular1,angular2);
|
||
|
|
}
|
||
|
|
btVector3 a2 = transB.getBasis() * m_rbBFrame.getOrigin();
|
||
|
|
{
|
||
|
|
btVector3* angular0 = (btVector3*)(info->m_J2angularAxis);
|
||
|
|
btVector3* angular1 = (btVector3*)(info->m_J2angularAxis+info->rowskip);
|
||
|
|
btVector3* angular2 = (btVector3*)(info->m_J2angularAxis+2*info->rowskip);
|
||
|
|
a2.getSkewSymmetricMatrix(angular0,angular1,angular2);
|
||
|
|
}
|
||
|
|
// set right hand side
|
||
|
|
btScalar linERP = (m_flags & BT_CONETWIST_FLAGS_LIN_ERP) ? m_linERP : info->erp;
|
||
|
|
btScalar k = info->fps * linERP;
|
||
|
|
int j;
|
||
|
|
for (j=0; j<3; j++)
|
||
|
|
{
|
||
|
|
info->m_constraintError[j*info->rowskip] = k * (a2[j] + transB.getOrigin()[j] - a1[j] - transA.getOrigin()[j]);
|
||
|
|
info->m_lowerLimit[j*info->rowskip] = -SIMD_INFINITY;
|
||
|
|
info->m_upperLimit[j*info->rowskip] = SIMD_INFINITY;
|
||
|
|
if(m_flags & BT_CONETWIST_FLAGS_LIN_CFM)
|
||
|
|
{
|
||
|
|
info->cfm[j*info->rowskip] = m_linCFM;
|
||
|
|
}
|
||
|
|
}
|
||
|
|
int row = 3;
|
||
|
|
int srow = row * info->rowskip;
|
||
|
|
btVector3 ax1;
|
||
|
|
// angular limits
|
||
|
|
if(m_solveSwingLimit)
|
||
|
|
{
|
||
|
|
btScalar *J1 = info->m_J1angularAxis;
|
||
|
|
btScalar *J2 = info->m_J2angularAxis;
|
||
|
|
if((m_swingSpan1 < m_fixThresh) && (m_swingSpan2 < m_fixThresh))
|
||
|
|
{
|
||
|
|
btTransform trA = transA*m_rbAFrame;
|
||
|
|
btVector3 p = trA.getBasis().getColumn(1);
|
||
|
|
btVector3 q = trA.getBasis().getColumn(2);
|
||
|
|
int srow1 = srow + info->rowskip;
|
||
|
|
J1[srow+0] = p[0];
|
||
|
|
J1[srow+1] = p[1];
|
||
|
|
J1[srow+2] = p[2];
|
||
|
|
J1[srow1+0] = q[0];
|
||
|
|
J1[srow1+1] = q[1];
|
||
|
|
J1[srow1+2] = q[2];
|
||
|
|
J2[srow+0] = -p[0];
|
||
|
|
J2[srow+1] = -p[1];
|
||
|
|
J2[srow+2] = -p[2];
|
||
|
|
J2[srow1+0] = -q[0];
|
||
|
|
J2[srow1+1] = -q[1];
|
||
|
|
J2[srow1+2] = -q[2];
|
||
|
|
btScalar fact = info->fps * m_relaxationFactor;
|
||
|
|
info->m_constraintError[srow] = fact * m_swingAxis.dot(p);
|
||
|
|
info->m_constraintError[srow1] = fact * m_swingAxis.dot(q);
|
||
|
|
info->m_lowerLimit[srow] = -SIMD_INFINITY;
|
||
|
|
info->m_upperLimit[srow] = SIMD_INFINITY;
|
||
|
|
info->m_lowerLimit[srow1] = -SIMD_INFINITY;
|
||
|
|
info->m_upperLimit[srow1] = SIMD_INFINITY;
|
||
|
|
srow = srow1 + info->rowskip;
|
||
|
|
}
|
||
|
|
else
|
||
|
|
{
|
||
|
|
ax1 = m_swingAxis * m_relaxationFactor * m_relaxationFactor;
|
||
|
|
J1[srow+0] = ax1[0];
|
||
|
|
J1[srow+1] = ax1[1];
|
||
|
|
J1[srow+2] = ax1[2];
|
||
|
|
J2[srow+0] = -ax1[0];
|
||
|
|
J2[srow+1] = -ax1[1];
|
||
|
|
J2[srow+2] = -ax1[2];
|
||
|
|
btScalar k = info->fps * m_biasFactor;
|
||
|
|
|
||
|
|
info->m_constraintError[srow] = k * m_swingCorrection;
|
||
|
|
if(m_flags & BT_CONETWIST_FLAGS_ANG_CFM)
|
||
|
|
{
|
||
|
|
info->cfm[srow] = m_angCFM;
|
||
|
|
}
|
||
|
|
// m_swingCorrection is always positive or 0
|
||
|
|
info->m_lowerLimit[srow] = 0;
|
||
|
|
info->m_upperLimit[srow] = SIMD_INFINITY;
|
||
|
|
srow += info->rowskip;
|
||
|
|
}
|
||
|
|
}
|
||
|
|
if(m_solveTwistLimit)
|
||
|
|
{
|
||
|
|
ax1 = m_twistAxis * m_relaxationFactor * m_relaxationFactor;
|
||
|
|
btScalar *J1 = info->m_J1angularAxis;
|
||
|
|
btScalar *J2 = info->m_J2angularAxis;
|
||
|
|
J1[srow+0] = ax1[0];
|
||
|
|
J1[srow+1] = ax1[1];
|
||
|
|
J1[srow+2] = ax1[2];
|
||
|
|
J2[srow+0] = -ax1[0];
|
||
|
|
J2[srow+1] = -ax1[1];
|
||
|
|
J2[srow+2] = -ax1[2];
|
||
|
|
btScalar k = info->fps * m_biasFactor;
|
||
|
|
info->m_constraintError[srow] = k * m_twistCorrection;
|
||
|
|
if(m_flags & BT_CONETWIST_FLAGS_ANG_CFM)
|
||
|
|
{
|
||
|
|
info->cfm[srow] = m_angCFM;
|
||
|
|
}
|
||
|
|
if(m_twistSpan > 0.0f)
|
||
|
|
{
|
||
|
|
|
||
|
|
if(m_twistCorrection > 0.0f)
|
||
|
|
{
|
||
|
|
info->m_lowerLimit[srow] = 0;
|
||
|
|
info->m_upperLimit[srow] = SIMD_INFINITY;
|
||
|
|
}
|
||
|
|
else
|
||
|
|
{
|
||
|
|
info->m_lowerLimit[srow] = -SIMD_INFINITY;
|
||
|
|
info->m_upperLimit[srow] = 0;
|
||
|
|
}
|
||
|
|
}
|
||
|
|
else
|
||
|
|
{
|
||
|
|
info->m_lowerLimit[srow] = -SIMD_INFINITY;
|
||
|
|
info->m_upperLimit[srow] = SIMD_INFINITY;
|
||
|
|
}
|
||
|
|
srow += info->rowskip;
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
|
||
|
|
|
||
|
|
void btConeTwistConstraint::buildJacobian()
|
||
|
|
{
|
||
|
|
if (m_useSolveConstraintObsolete)
|
||
|
|
{
|
||
|
|
m_appliedImpulse = btScalar(0.);
|
||
|
|
m_accTwistLimitImpulse = btScalar(0.);
|
||
|
|
m_accSwingLimitImpulse = btScalar(0.);
|
||
|
|
m_accMotorImpulse = btVector3(0.,0.,0.);
|
||
|
|
|
||
|
|
if (!m_angularOnly)
|
||
|
|
{
|
||
|
|
btVector3 pivotAInW = m_rbA.getCenterOfMassTransform()*m_rbAFrame.getOrigin();
|
||
|
|
btVector3 pivotBInW = m_rbB.getCenterOfMassTransform()*m_rbBFrame.getOrigin();
|
||
|
|
btVector3 relPos = pivotBInW - pivotAInW;
|
||
|
|
|
||
|
|
btVector3 normal[3];
|
||
|
|
if (relPos.length2() > SIMD_EPSILON)
|
||
|
|
{
|
||
|
|
normal[0] = relPos.normalized();
|
||
|
|
}
|
||
|
|
else
|
||
|
|
{
|
||
|
|
normal[0].setValue(btScalar(1.0),0,0);
|
||
|
|
}
|
||
|
|
|
||
|
|
btPlaneSpace1(normal[0], normal[1], normal[2]);
|
||
|
|
|
||
|
|
for (int i=0;i<3;i++)
|
||
|
|
{
|
||
|
|
new (&m_jac[i]) btJacobianEntry(
|
||
|
|
m_rbA.getCenterOfMassTransform().getBasis().transpose(),
|
||
|
|
m_rbB.getCenterOfMassTransform().getBasis().transpose(),
|
||
|
|
pivotAInW - m_rbA.getCenterOfMassPosition(),
|
||
|
|
pivotBInW - m_rbB.getCenterOfMassPosition(),
|
||
|
|
normal[i],
|
||
|
|
m_rbA.getInvInertiaDiagLocal(),
|
||
|
|
m_rbA.getInvMass(),
|
||
|
|
m_rbB.getInvInertiaDiagLocal(),
|
||
|
|
m_rbB.getInvMass());
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
calcAngleInfo2(m_rbA.getCenterOfMassTransform(),m_rbB.getCenterOfMassTransform(),m_rbA.getInvInertiaTensorWorld(),m_rbB.getInvInertiaTensorWorld());
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
|
||
|
|
|
||
|
|
void btConeTwistConstraint::solveConstraintObsolete(btRigidBody& bodyA,btRigidBody& bodyB,btScalar timeStep)
|
||
|
|
{
|
||
|
|
#ifndef __SPU__
|
||
|
|
if (m_useSolveConstraintObsolete)
|
||
|
|
{
|
||
|
|
btVector3 pivotAInW = m_rbA.getCenterOfMassTransform()*m_rbAFrame.getOrigin();
|
||
|
|
btVector3 pivotBInW = m_rbB.getCenterOfMassTransform()*m_rbBFrame.getOrigin();
|
||
|
|
|
||
|
|
btScalar tau = btScalar(0.3);
|
||
|
|
|
||
|
|
//linear part
|
||
|
|
if (!m_angularOnly)
|
||
|
|
{
|
||
|
|
btVector3 rel_pos1 = pivotAInW - m_rbA.getCenterOfMassPosition();
|
||
|
|
btVector3 rel_pos2 = pivotBInW - m_rbB.getCenterOfMassPosition();
|
||
|
|
|
||
|
|
btVector3 vel1;
|
||
|
|
bodyA.internalGetVelocityInLocalPointObsolete(rel_pos1,vel1);
|
||
|
|
btVector3 vel2;
|
||
|
|
bodyB.internalGetVelocityInLocalPointObsolete(rel_pos2,vel2);
|
||
|
|
btVector3 vel = vel1 - vel2;
|
||
|
|
|
||
|
|
for (int i=0;i<3;i++)
|
||
|
|
{
|
||
|
|
const btVector3& normal = m_jac[i].m_linearJointAxis;
|
||
|
|
btScalar jacDiagABInv = btScalar(1.) / m_jac[i].getDiagonal();
|
||
|
|
|
||
|
|
btScalar rel_vel;
|
||
|
|
rel_vel = normal.dot(vel);
|
||
|
|
//positional error (zeroth order error)
|
||
|
|
btScalar depth = -(pivotAInW - pivotBInW).dot(normal); //this is the error projected on the normal
|
||
|
|
btScalar impulse = depth*tau/timeStep * jacDiagABInv - rel_vel * jacDiagABInv;
|
||
|
|
m_appliedImpulse += impulse;
|
||
|
|
|
||
|
|
btVector3 ftorqueAxis1 = rel_pos1.cross(normal);
|
||
|
|
btVector3 ftorqueAxis2 = rel_pos2.cross(normal);
|
||
|
|
bodyA.internalApplyImpulse(normal*m_rbA.getInvMass(), m_rbA.getInvInertiaTensorWorld()*ftorqueAxis1,impulse);
|
||
|
|
bodyB.internalApplyImpulse(normal*m_rbB.getInvMass(), m_rbB.getInvInertiaTensorWorld()*ftorqueAxis2,-impulse);
|
||
|
|
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
// apply motor
|
||
|
|
if (m_bMotorEnabled)
|
||
|
|
{
|
||
|
|
// compute current and predicted transforms
|
||
|
|
btTransform trACur = m_rbA.getCenterOfMassTransform();
|
||
|
|
btTransform trBCur = m_rbB.getCenterOfMassTransform();
|
||
|
|
btVector3 omegaA; bodyA.internalGetAngularVelocity(omegaA);
|
||
|
|
btVector3 omegaB; bodyB.internalGetAngularVelocity(omegaB);
|
||
|
|
btTransform trAPred; trAPred.setIdentity();
|
||
|
|
btVector3 zerovec(0,0,0);
|
||
|
|
btTransformUtil::integrateTransform(
|
||
|
|
trACur, zerovec, omegaA, timeStep, trAPred);
|
||
|
|
btTransform trBPred; trBPred.setIdentity();
|
||
|
|
btTransformUtil::integrateTransform(
|
||
|
|
trBCur, zerovec, omegaB, timeStep, trBPred);
|
||
|
|
|
||
|
|
// compute desired transforms in world
|
||
|
|
btTransform trPose(m_qTarget);
|
||
|
|
btTransform trABDes = m_rbBFrame * trPose * m_rbAFrame.inverse();
|
||
|
|
btTransform trADes = trBPred * trABDes;
|
||
|
|
btTransform trBDes = trAPred * trABDes.inverse();
|
||
|
|
|
||
|
|
// compute desired omegas in world
|
||
|
|
btVector3 omegaADes, omegaBDes;
|
||
|
|
|
||
|
|
btTransformUtil::calculateVelocity(trACur, trADes, timeStep, zerovec, omegaADes);
|
||
|
|
btTransformUtil::calculateVelocity(trBCur, trBDes, timeStep, zerovec, omegaBDes);
|
||
|
|
|
||
|
|
// compute delta omegas
|
||
|
|
btVector3 dOmegaA = omegaADes - omegaA;
|
||
|
|
btVector3 dOmegaB = omegaBDes - omegaB;
|
||
|
|
|
||
|
|
// compute weighted avg axis of dOmega (weighting based on inertias)
|
||
|
|
btVector3 axisA, axisB;
|
||
|
|
btScalar kAxisAInv = 0, kAxisBInv = 0;
|
||
|
|
|
||
|
|
if (dOmegaA.length2() > SIMD_EPSILON)
|
||
|
|
{
|
||
|
|
axisA = dOmegaA.normalized();
|
||
|
|
kAxisAInv = getRigidBodyA().computeAngularImpulseDenominator(axisA);
|
||
|
|
}
|
||
|
|
|
||
|
|
if (dOmegaB.length2() > SIMD_EPSILON)
|
||
|
|
{
|
||
|
|
axisB = dOmegaB.normalized();
|
||
|
|
kAxisBInv = getRigidBodyB().computeAngularImpulseDenominator(axisB);
|
||
|
|
}
|
||
|
|
|
||
|
|
btVector3 avgAxis = kAxisAInv * axisA + kAxisBInv * axisB;
|
||
|
|
|
||
|
|
static bool bDoTorque = true;
|
||
|
|
if (bDoTorque && avgAxis.length2() > SIMD_EPSILON)
|
||
|
|
{
|
||
|
|
avgAxis.normalize();
|
||
|
|
kAxisAInv = getRigidBodyA().computeAngularImpulseDenominator(avgAxis);
|
||
|
|
kAxisBInv = getRigidBodyB().computeAngularImpulseDenominator(avgAxis);
|
||
|
|
btScalar kInvCombined = kAxisAInv + kAxisBInv;
|
||
|
|
|
||
|
|
btVector3 impulse = (kAxisAInv * dOmegaA - kAxisBInv * dOmegaB) /
|
||
|
|
(kInvCombined * kInvCombined);
|
||
|
|
|
||
|
|
if (m_maxMotorImpulse >= 0)
|
||
|
|
{
|
||
|
|
btScalar fMaxImpulse = m_maxMotorImpulse;
|
||
|
|
if (m_bNormalizedMotorStrength)
|
||
|
|
fMaxImpulse = fMaxImpulse/kAxisAInv;
|
||
|
|
|
||
|
|
btVector3 newUnclampedAccImpulse = m_accMotorImpulse + impulse;
|
||
|
|
btScalar newUnclampedMag = newUnclampedAccImpulse.length();
|
||
|
|
if (newUnclampedMag > fMaxImpulse)
|
||
|
|
{
|
||
|
|
newUnclampedAccImpulse.normalize();
|
||
|
|
newUnclampedAccImpulse *= fMaxImpulse;
|
||
|
|
impulse = newUnclampedAccImpulse - m_accMotorImpulse;
|
||
|
|
}
|
||
|
|
m_accMotorImpulse += impulse;
|
||
|
|
}
|
||
|
|
|
||
|
|
btScalar impulseMag = impulse.length();
|
||
|
|
btVector3 impulseAxis = impulse / impulseMag;
|
||
|
|
|
||
|
|
bodyA.internalApplyImpulse(btVector3(0,0,0), m_rbA.getInvInertiaTensorWorld()*impulseAxis, impulseMag);
|
||
|
|
bodyB.internalApplyImpulse(btVector3(0,0,0), m_rbB.getInvInertiaTensorWorld()*impulseAxis, -impulseMag);
|
||
|
|
|
||
|
|
}
|
||
|
|
}
|
||
|
|
else if (m_damping > SIMD_EPSILON) // no motor: do a little damping
|
||
|
|
{
|
||
|
|
btVector3 angVelA; bodyA.internalGetAngularVelocity(angVelA);
|
||
|
|
btVector3 angVelB; bodyB.internalGetAngularVelocity(angVelB);
|
||
|
|
btVector3 relVel = angVelB - angVelA;
|
||
|
|
if (relVel.length2() > SIMD_EPSILON)
|
||
|
|
{
|
||
|
|
btVector3 relVelAxis = relVel.normalized();
|
||
|
|
btScalar m_kDamping = btScalar(1.) /
|
||
|
|
(getRigidBodyA().computeAngularImpulseDenominator(relVelAxis) +
|
||
|
|
getRigidBodyB().computeAngularImpulseDenominator(relVelAxis));
|
||
|
|
btVector3 impulse = m_damping * m_kDamping * relVel;
|
||
|
|
|
||
|
|
btScalar impulseMag = impulse.length();
|
||
|
|
btVector3 impulseAxis = impulse / impulseMag;
|
||
|
|
bodyA.internalApplyImpulse(btVector3(0,0,0), m_rbA.getInvInertiaTensorWorld()*impulseAxis, impulseMag);
|
||
|
|
bodyB.internalApplyImpulse(btVector3(0,0,0), m_rbB.getInvInertiaTensorWorld()*impulseAxis, -impulseMag);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
// joint limits
|
||
|
|
{
|
||
|
|
///solve angular part
|
||
|
|
btVector3 angVelA;
|
||
|
|
bodyA.internalGetAngularVelocity(angVelA);
|
||
|
|
btVector3 angVelB;
|
||
|
|
bodyB.internalGetAngularVelocity(angVelB);
|
||
|
|
|
||
|
|
// solve swing limit
|
||
|
|
if (m_solveSwingLimit)
|
||
|
|
{
|
||
|
|
btScalar amplitude = m_swingLimitRatio * m_swingCorrection*m_biasFactor/timeStep;
|
||
|
|
btScalar relSwingVel = (angVelB - angVelA).dot(m_swingAxis);
|
||
|
|
if (relSwingVel > 0)
|
||
|
|
amplitude += m_swingLimitRatio * relSwingVel * m_relaxationFactor;
|
||
|
|
btScalar impulseMag = amplitude * m_kSwing;
|
||
|
|
|
||
|
|
// Clamp the accumulated impulse
|
||
|
|
btScalar temp = m_accSwingLimitImpulse;
|
||
|
|
m_accSwingLimitImpulse = btMax(m_accSwingLimitImpulse + impulseMag, btScalar(0.0) );
|
||
|
|
impulseMag = m_accSwingLimitImpulse - temp;
|
||
|
|
|
||
|
|
btVector3 impulse = m_swingAxis * impulseMag;
|
||
|
|
|
||
|
|
// don't let cone response affect twist
|
||
|
|
// (this can happen since body A's twist doesn't match body B's AND we use an elliptical cone limit)
|
||
|
|
{
|
||
|
|
btVector3 impulseTwistCouple = impulse.dot(m_twistAxisA) * m_twistAxisA;
|
||
|
|
btVector3 impulseNoTwistCouple = impulse - impulseTwistCouple;
|
||
|
|
impulse = impulseNoTwistCouple;
|
||
|
|
}
|
||
|
|
|
||
|
|
impulseMag = impulse.length();
|
||
|
|
btVector3 noTwistSwingAxis = impulse / impulseMag;
|
||
|
|
|
||
|
|
bodyA.internalApplyImpulse(btVector3(0,0,0), m_rbA.getInvInertiaTensorWorld()*noTwistSwingAxis, impulseMag);
|
||
|
|
bodyB.internalApplyImpulse(btVector3(0,0,0), m_rbB.getInvInertiaTensorWorld()*noTwistSwingAxis, -impulseMag);
|
||
|
|
}
|
||
|
|
|
||
|
|
|
||
|
|
// solve twist limit
|
||
|
|
if (m_solveTwistLimit)
|
||
|
|
{
|
||
|
|
btScalar amplitude = m_twistLimitRatio * m_twistCorrection*m_biasFactor/timeStep;
|
||
|
|
btScalar relTwistVel = (angVelB - angVelA).dot( m_twistAxis );
|
||
|
|
if (relTwistVel > 0) // only damp when moving towards limit (m_twistAxis flipping is important)
|
||
|
|
amplitude += m_twistLimitRatio * relTwistVel * m_relaxationFactor;
|
||
|
|
btScalar impulseMag = amplitude * m_kTwist;
|
||
|
|
|
||
|
|
// Clamp the accumulated impulse
|
||
|
|
btScalar temp = m_accTwistLimitImpulse;
|
||
|
|
m_accTwistLimitImpulse = btMax(m_accTwistLimitImpulse + impulseMag, btScalar(0.0) );
|
||
|
|
impulseMag = m_accTwistLimitImpulse - temp;
|
||
|
|
|
||
|
|
btVector3 impulse = m_twistAxis * impulseMag;
|
||
|
|
|
||
|
|
bodyA.internalApplyImpulse(btVector3(0,0,0), m_rbA.getInvInertiaTensorWorld()*m_twistAxis,impulseMag);
|
||
|
|
bodyB.internalApplyImpulse(btVector3(0,0,0), m_rbB.getInvInertiaTensorWorld()*m_twistAxis,-impulseMag);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
}
|
||
|
|
#else
|
||
|
|
btAssert(0);
|
||
|
|
#endif //__SPU__
|
||
|
|
}
|
||
|
|
|
||
|
|
|
||
|
|
|
||
|
|
|
||
|
|
void btConeTwistConstraint::updateRHS(btScalar timeStep)
|
||
|
|
{
|
||
|
|
(void)timeStep;
|
||
|
|
|
||
|
|
}
|
||
|
|
|
||
|
|
|
||
|
|
#ifndef __SPU__
|
||
|
|
void btConeTwistConstraint::calcAngleInfo()
|
||
|
|
{
|
||
|
|
m_swingCorrection = btScalar(0.);
|
||
|
|
m_twistLimitSign = btScalar(0.);
|
||
|
|
m_solveTwistLimit = false;
|
||
|
|
m_solveSwingLimit = false;
|
||
|
|
|
||
|
|
btVector3 b1Axis1,b1Axis2,b1Axis3;
|
||
|
|
btVector3 b2Axis1,b2Axis2;
|
||
|
|
|
||
|
|
b1Axis1 = getRigidBodyA().getCenterOfMassTransform().getBasis() * this->m_rbAFrame.getBasis().getColumn(0);
|
||
|
|
b2Axis1 = getRigidBodyB().getCenterOfMassTransform().getBasis() * this->m_rbBFrame.getBasis().getColumn(0);
|
||
|
|
|
||
|
|
btScalar swing1=btScalar(0.),swing2 = btScalar(0.);
|
||
|
|
|
||
|
|
btScalar swx=btScalar(0.),swy = btScalar(0.);
|
||
|
|
btScalar thresh = btScalar(10.);
|
||
|
|
btScalar fact;
|
||
|
|
|
||
|
|
// Get Frame into world space
|
||
|
|
if (m_swingSpan1 >= btScalar(0.05f))
|
||
|
|
{
|
||
|
|
b1Axis2 = getRigidBodyA().getCenterOfMassTransform().getBasis() * this->m_rbAFrame.getBasis().getColumn(1);
|
||
|
|
swx = b2Axis1.dot(b1Axis1);
|
||
|
|
swy = b2Axis1.dot(b1Axis2);
|
||
|
|
swing1 = btAtan2Fast(swy, swx);
|
||
|
|
fact = (swy*swy + swx*swx) * thresh * thresh;
|
||
|
|
fact = fact / (fact + btScalar(1.0));
|
||
|
|
swing1 *= fact;
|
||
|
|
}
|
||
|
|
|
||
|
|
if (m_swingSpan2 >= btScalar(0.05f))
|
||
|
|
{
|
||
|
|
b1Axis3 = getRigidBodyA().getCenterOfMassTransform().getBasis() * this->m_rbAFrame.getBasis().getColumn(2);
|
||
|
|
swx = b2Axis1.dot(b1Axis1);
|
||
|
|
swy = b2Axis1.dot(b1Axis3);
|
||
|
|
swing2 = btAtan2Fast(swy, swx);
|
||
|
|
fact = (swy*swy + swx*swx) * thresh * thresh;
|
||
|
|
fact = fact / (fact + btScalar(1.0));
|
||
|
|
swing2 *= fact;
|
||
|
|
}
|
||
|
|
|
||
|
|
btScalar RMaxAngle1Sq = 1.0f / (m_swingSpan1*m_swingSpan1);
|
||
|
|
btScalar RMaxAngle2Sq = 1.0f / (m_swingSpan2*m_swingSpan2);
|
||
|
|
btScalar EllipseAngle = btFabs(swing1*swing1)* RMaxAngle1Sq + btFabs(swing2*swing2) * RMaxAngle2Sq;
|
||
|
|
|
||
|
|
if (EllipseAngle > 1.0f)
|
||
|
|
{
|
||
|
|
m_swingCorrection = EllipseAngle-1.0f;
|
||
|
|
m_solveSwingLimit = true;
|
||
|
|
// Calculate necessary axis & factors
|
||
|
|
m_swingAxis = b2Axis1.cross(b1Axis2* b2Axis1.dot(b1Axis2) + b1Axis3* b2Axis1.dot(b1Axis3));
|
||
|
|
m_swingAxis.normalize();
|
||
|
|
btScalar swingAxisSign = (b2Axis1.dot(b1Axis1) >= 0.0f) ? 1.0f : -1.0f;
|
||
|
|
m_swingAxis *= swingAxisSign;
|
||
|
|
}
|
||
|
|
|
||
|
|
// Twist limits
|
||
|
|
if (m_twistSpan >= btScalar(0.))
|
||
|
|
{
|
||
|
|
btVector3 b2Axis2 = getRigidBodyB().getCenterOfMassTransform().getBasis() * this->m_rbBFrame.getBasis().getColumn(1);
|
||
|
|
btQuaternion rotationArc = shortestArcQuat(b2Axis1,b1Axis1);
|
||
|
|
btVector3 TwistRef = quatRotate(rotationArc,b2Axis2);
|
||
|
|
btScalar twist = btAtan2Fast( TwistRef.dot(b1Axis3), TwistRef.dot(b1Axis2) );
|
||
|
|
m_twistAngle = twist;
|
||
|
|
|
||
|
|
// btScalar lockedFreeFactor = (m_twistSpan > btScalar(0.05f)) ? m_limitSoftness : btScalar(0.);
|
||
|
|
btScalar lockedFreeFactor = (m_twistSpan > btScalar(0.05f)) ? btScalar(1.0f) : btScalar(0.);
|
||
|
|
if (twist <= -m_twistSpan*lockedFreeFactor)
|
||
|
|
{
|
||
|
|
m_twistCorrection = -(twist + m_twistSpan);
|
||
|
|
m_solveTwistLimit = true;
|
||
|
|
m_twistAxis = (b2Axis1 + b1Axis1) * 0.5f;
|
||
|
|
m_twistAxis.normalize();
|
||
|
|
m_twistAxis *= -1.0f;
|
||
|
|
}
|
||
|
|
else if (twist > m_twistSpan*lockedFreeFactor)
|
||
|
|
{
|
||
|
|
m_twistCorrection = (twist - m_twistSpan);
|
||
|
|
m_solveTwistLimit = true;
|
||
|
|
m_twistAxis = (b2Axis1 + b1Axis1) * 0.5f;
|
||
|
|
m_twistAxis.normalize();
|
||
|
|
}
|
||
|
|
}
|
||
|
|
}
|
||
|
|
#endif //__SPU__
|
||
|
|
|
||
|
|
static btVector3 vTwist(1,0,0); // twist axis in constraint's space
|
||
|
|
|
||
|
|
|
||
|
|
|
||
|
|
void btConeTwistConstraint::calcAngleInfo2(const btTransform& transA, const btTransform& transB, const btMatrix3x3& invInertiaWorldA,const btMatrix3x3& invInertiaWorldB)
|
||
|
|
{
|
||
|
|
m_swingCorrection = btScalar(0.);
|
||
|
|
m_twistLimitSign = btScalar(0.);
|
||
|
|
m_solveTwistLimit = false;
|
||
|
|
m_solveSwingLimit = false;
|
||
|
|
// compute rotation of A wrt B (in constraint space)
|
||
|
|
if (m_bMotorEnabled && (!m_useSolveConstraintObsolete))
|
||
|
|
{ // it is assumed that setMotorTarget() was alredy called
|
||
|
|
// and motor target m_qTarget is within constraint limits
|
||
|
|
// TODO : split rotation to pure swing and pure twist
|
||
|
|
// compute desired transforms in world
|
||
|
|
btTransform trPose(m_qTarget);
|
||
|
|
btTransform trA = transA * m_rbAFrame;
|
||
|
|
btTransform trB = transB * m_rbBFrame;
|
||
|
|
btTransform trDeltaAB = trB * trPose * trA.inverse();
|
||
|
|
btQuaternion qDeltaAB = trDeltaAB.getRotation();
|
||
|
|
btVector3 swingAxis = btVector3(qDeltaAB.x(), qDeltaAB.y(), qDeltaAB.z());
|
||
|
|
m_swingAxis = swingAxis;
|
||
|
|
m_swingAxis.normalize();
|
||
|
|
m_swingCorrection = qDeltaAB.getAngle();
|
||
|
|
if(!btFuzzyZero(m_swingCorrection))
|
||
|
|
{
|
||
|
|
m_solveSwingLimit = true;
|
||
|
|
}
|
||
|
|
return;
|
||
|
|
}
|
||
|
|
|
||
|
|
|
||
|
|
{
|
||
|
|
// compute rotation of A wrt B (in constraint space)
|
||
|
|
btQuaternion qA = transA.getRotation() * m_rbAFrame.getRotation();
|
||
|
|
btQuaternion qB = transB.getRotation() * m_rbBFrame.getRotation();
|
||
|
|
btQuaternion qAB = qB.inverse() * qA;
|
||
|
|
// split rotation into cone and twist
|
||
|
|
// (all this is done from B's perspective. Maybe I should be averaging axes...)
|
||
|
|
btVector3 vConeNoTwist = quatRotate(qAB, vTwist); vConeNoTwist.normalize();
|
||
|
|
btQuaternion qABCone = shortestArcQuat(vTwist, vConeNoTwist); qABCone.normalize();
|
||
|
|
btQuaternion qABTwist = qABCone.inverse() * qAB; qABTwist.normalize();
|
||
|
|
|
||
|
|
if (m_swingSpan1 >= m_fixThresh && m_swingSpan2 >= m_fixThresh)
|
||
|
|
{
|
||
|
|
btScalar swingAngle, swingLimit = 0; btVector3 swingAxis;
|
||
|
|
computeConeLimitInfo(qABCone, swingAngle, swingAxis, swingLimit);
|
||
|
|
|
||
|
|
if (swingAngle > swingLimit * m_limitSoftness)
|
||
|
|
{
|
||
|
|
m_solveSwingLimit = true;
|
||
|
|
|
||
|
|
// compute limit ratio: 0->1, where
|
||
|
|
// 0 == beginning of soft limit
|
||
|
|
// 1 == hard/real limit
|
||
|
|
m_swingLimitRatio = 1.f;
|
||
|
|
if (swingAngle < swingLimit && m_limitSoftness < 1.f - SIMD_EPSILON)
|
||
|
|
{
|
||
|
|
m_swingLimitRatio = (swingAngle - swingLimit * m_limitSoftness)/
|
||
|
|
(swingLimit - swingLimit * m_limitSoftness);
|
||
|
|
}
|
||
|
|
|
||
|
|
// swing correction tries to get back to soft limit
|
||
|
|
m_swingCorrection = swingAngle - (swingLimit * m_limitSoftness);
|
||
|
|
|
||
|
|
// adjustment of swing axis (based on ellipse normal)
|
||
|
|
adjustSwingAxisToUseEllipseNormal(swingAxis);
|
||
|
|
|
||
|
|
// Calculate necessary axis & factors
|
||
|
|
m_swingAxis = quatRotate(qB, -swingAxis);
|
||
|
|
|
||
|
|
m_twistAxisA.setValue(0,0,0);
|
||
|
|
|
||
|
|
m_kSwing = btScalar(1.) /
|
||
|
|
(computeAngularImpulseDenominator(m_swingAxis,invInertiaWorldA) +
|
||
|
|
computeAngularImpulseDenominator(m_swingAxis,invInertiaWorldB));
|
||
|
|
}
|
||
|
|
}
|
||
|
|
else
|
||
|
|
{
|
||
|
|
// you haven't set any limits;
|
||
|
|
// or you're trying to set at least one of the swing limits too small. (if so, do you really want a conetwist constraint?)
|
||
|
|
// anyway, we have either hinge or fixed joint
|
||
|
|
btVector3 ivA = transA.getBasis() * m_rbAFrame.getBasis().getColumn(0);
|
||
|
|
btVector3 jvA = transA.getBasis() * m_rbAFrame.getBasis().getColumn(1);
|
||
|
|
btVector3 kvA = transA.getBasis() * m_rbAFrame.getBasis().getColumn(2);
|
||
|
|
btVector3 ivB = transB.getBasis() * m_rbBFrame.getBasis().getColumn(0);
|
||
|
|
btVector3 target;
|
||
|
|
btScalar x = ivB.dot(ivA);
|
||
|
|
btScalar y = ivB.dot(jvA);
|
||
|
|
btScalar z = ivB.dot(kvA);
|
||
|
|
if((m_swingSpan1 < m_fixThresh) && (m_swingSpan2 < m_fixThresh))
|
||
|
|
{ // fixed. We'll need to add one more row to constraint
|
||
|
|
if((!btFuzzyZero(y)) || (!(btFuzzyZero(z))))
|
||
|
|
{
|
||
|
|
m_solveSwingLimit = true;
|
||
|
|
m_swingAxis = -ivB.cross(ivA);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
else
|
||
|
|
{
|
||
|
|
if(m_swingSpan1 < m_fixThresh)
|
||
|
|
{ // hinge around Y axis
|
||
|
|
if(!(btFuzzyZero(y)))
|
||
|
|
{
|
||
|
|
m_solveSwingLimit = true;
|
||
|
|
if(m_swingSpan2 >= m_fixThresh)
|
||
|
|
{
|
||
|
|
y = btScalar(0.f);
|
||
|
|
btScalar span2 = btAtan2(z, x);
|
||
|
|
if(span2 > m_swingSpan2)
|
||
|
|
{
|
||
|
|
x = btCos(m_swingSpan2);
|
||
|
|
z = btSin(m_swingSpan2);
|
||
|
|
}
|
||
|
|
else if(span2 < -m_swingSpan2)
|
||
|
|
{
|
||
|
|
x = btCos(m_swingSpan2);
|
||
|
|
z = -btSin(m_swingSpan2);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
}
|
||
|
|
}
|
||
|
|
else
|
||
|
|
{ // hinge around Z axis
|
||
|
|
if(!btFuzzyZero(z))
|
||
|
|
{
|
||
|
|
m_solveSwingLimit = true;
|
||
|
|
if(m_swingSpan1 >= m_fixThresh)
|
||
|
|
{
|
||
|
|
z = btScalar(0.f);
|
||
|
|
btScalar span1 = btAtan2(y, x);
|
||
|
|
if(span1 > m_swingSpan1)
|
||
|
|
{
|
||
|
|
x = btCos(m_swingSpan1);
|
||
|
|
y = btSin(m_swingSpan1);
|
||
|
|
}
|
||
|
|
else if(span1 < -m_swingSpan1)
|
||
|
|
{
|
||
|
|
x = btCos(m_swingSpan1);
|
||
|
|
y = -btSin(m_swingSpan1);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
}
|
||
|
|
}
|
||
|
|
target[0] = x * ivA[0] + y * jvA[0] + z * kvA[0];
|
||
|
|
target[1] = x * ivA[1] + y * jvA[1] + z * kvA[1];
|
||
|
|
target[2] = x * ivA[2] + y * jvA[2] + z * kvA[2];
|
||
|
|
target.normalize();
|
||
|
|
m_swingAxis = -ivB.cross(target);
|
||
|
|
m_swingCorrection = m_swingAxis.length();
|
||
|
|
m_swingAxis.normalize();
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
if (m_twistSpan >= btScalar(0.f))
|
||
|
|
{
|
||
|
|
btVector3 twistAxis;
|
||
|
|
computeTwistLimitInfo(qABTwist, m_twistAngle, twistAxis);
|
||
|
|
|
||
|
|
if (m_twistAngle > m_twistSpan*m_limitSoftness)
|
||
|
|
{
|
||
|
|
m_solveTwistLimit = true;
|
||
|
|
|
||
|
|
m_twistLimitRatio = 1.f;
|
||
|
|
if (m_twistAngle < m_twistSpan && m_limitSoftness < 1.f - SIMD_EPSILON)
|
||
|
|
{
|
||
|
|
m_twistLimitRatio = (m_twistAngle - m_twistSpan * m_limitSoftness)/
|
||
|
|
(m_twistSpan - m_twistSpan * m_limitSoftness);
|
||
|
|
}
|
||
|
|
|
||
|
|
// twist correction tries to get back to soft limit
|
||
|
|
m_twistCorrection = m_twistAngle - (m_twistSpan * m_limitSoftness);
|
||
|
|
|
||
|
|
m_twistAxis = quatRotate(qB, -twistAxis);
|
||
|
|
|
||
|
|
m_kTwist = btScalar(1.) /
|
||
|
|
(computeAngularImpulseDenominator(m_twistAxis,invInertiaWorldA) +
|
||
|
|
computeAngularImpulseDenominator(m_twistAxis,invInertiaWorldB));
|
||
|
|
}
|
||
|
|
|
||
|
|
if (m_solveSwingLimit)
|
||
|
|
m_twistAxisA = quatRotate(qA, -twistAxis);
|
||
|
|
}
|
||
|
|
else
|
||
|
|
{
|
||
|
|
m_twistAngle = btScalar(0.f);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
|
||
|
|
|
||
|
|
// given a cone rotation in constraint space, (pre: twist must already be removed)
|
||
|
|
// this method computes its corresponding swing angle and axis.
|
||
|
|
// more interestingly, it computes the cone/swing limit (angle) for this cone "pose".
|
||
|
|
void btConeTwistConstraint::computeConeLimitInfo(const btQuaternion& qCone,
|
||
|
|
btScalar& swingAngle, // out
|
||
|
|
btVector3& vSwingAxis, // out
|
||
|
|
btScalar& swingLimit) // out
|
||
|
|
{
|
||
|
|
swingAngle = qCone.getAngle();
|
||
|
|
if (swingAngle > SIMD_EPSILON)
|
||
|
|
{
|
||
|
|
vSwingAxis = btVector3(qCone.x(), qCone.y(), qCone.z());
|
||
|
|
vSwingAxis.normalize();
|
||
|
|
if (fabs(vSwingAxis.x()) > SIMD_EPSILON)
|
||
|
|
{
|
||
|
|
// non-zero twist?! this should never happen.
|
||
|
|
int wtf = 0; wtf = wtf;
|
||
|
|
}
|
||
|
|
|
||
|
|
// Compute limit for given swing. tricky:
|
||
|
|
// Given a swing axis, we're looking for the intersection with the bounding cone ellipse.
|
||
|
|
// (Since we're dealing with angles, this ellipse is embedded on the surface of a sphere.)
|
||
|
|
|
||
|
|
// For starters, compute the direction from center to surface of ellipse.
|
||
|
|
// This is just the perpendicular (ie. rotate 2D vector by PI/2) of the swing axis.
|
||
|
|
// (vSwingAxis is the cone rotation (in z,y); change vars and rotate to (x,y) coords.)
|
||
|
|
btScalar xEllipse = vSwingAxis.y();
|
||
|
|
btScalar yEllipse = -vSwingAxis.z();
|
||
|
|
|
||
|
|
// Now, we use the slope of the vector (using x/yEllipse) and find the length
|
||
|
|
// of the line that intersects the ellipse:
|
||
|
|
// x^2 y^2
|
||
|
|
// --- + --- = 1, where a and b are semi-major axes 2 and 1 respectively (ie. the limits)
|
||
|
|
// a^2 b^2
|
||
|
|
// Do the math and it should be clear.
|
||
|
|
|
||
|
|
swingLimit = m_swingSpan1; // if xEllipse == 0, we have a pure vSwingAxis.z rotation: just use swingspan1
|
||
|
|
if (fabs(xEllipse) > SIMD_EPSILON)
|
||
|
|
{
|
||
|
|
btScalar surfaceSlope2 = (yEllipse*yEllipse)/(xEllipse*xEllipse);
|
||
|
|
btScalar norm = 1 / (m_swingSpan2 * m_swingSpan2);
|
||
|
|
norm += surfaceSlope2 / (m_swingSpan1 * m_swingSpan1);
|
||
|
|
btScalar swingLimit2 = (1 + surfaceSlope2) / norm;
|
||
|
|
swingLimit = sqrt(swingLimit2);
|
||
|
|
}
|
||
|
|
|
||
|
|
// test!
|
||
|
|
/*swingLimit = m_swingSpan2;
|
||
|
|
if (fabs(vSwingAxis.z()) > SIMD_EPSILON)
|
||
|
|
{
|
||
|
|
btScalar mag_2 = m_swingSpan1*m_swingSpan1 + m_swingSpan2*m_swingSpan2;
|
||
|
|
btScalar sinphi = m_swingSpan2 / sqrt(mag_2);
|
||
|
|
btScalar phi = asin(sinphi);
|
||
|
|
btScalar theta = atan2(fabs(vSwingAxis.y()),fabs(vSwingAxis.z()));
|
||
|
|
btScalar alpha = 3.14159f - theta - phi;
|
||
|
|
btScalar sinalpha = sin(alpha);
|
||
|
|
swingLimit = m_swingSpan1 * sinphi/sinalpha;
|
||
|
|
}*/
|
||
|
|
}
|
||
|
|
else if (swingAngle < 0)
|
||
|
|
{
|
||
|
|
// this should never happen!
|
||
|
|
int wtf = 0; wtf = wtf;
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
btVector3 btConeTwistConstraint::GetPointForAngle(btScalar fAngleInRadians, btScalar fLength) const
|
||
|
|
{
|
||
|
|
// compute x/y in ellipse using cone angle (0 -> 2*PI along surface of cone)
|
||
|
|
btScalar xEllipse = btCos(fAngleInRadians);
|
||
|
|
btScalar yEllipse = btSin(fAngleInRadians);
|
||
|
|
|
||
|
|
// Use the slope of the vector (using x/yEllipse) and find the length
|
||
|
|
// of the line that intersects the ellipse:
|
||
|
|
// x^2 y^2
|
||
|
|
// --- + --- = 1, where a and b are semi-major axes 2 and 1 respectively (ie. the limits)
|
||
|
|
// a^2 b^2
|
||
|
|
// Do the math and it should be clear.
|
||
|
|
|
||
|
|
float swingLimit = m_swingSpan1; // if xEllipse == 0, just use axis b (1)
|
||
|
|
if (fabs(xEllipse) > SIMD_EPSILON)
|
||
|
|
{
|
||
|
|
btScalar surfaceSlope2 = (yEllipse*yEllipse)/(xEllipse*xEllipse);
|
||
|
|
btScalar norm = 1 / (m_swingSpan2 * m_swingSpan2);
|
||
|
|
norm += surfaceSlope2 / (m_swingSpan1 * m_swingSpan1);
|
||
|
|
btScalar swingLimit2 = (1 + surfaceSlope2) / norm;
|
||
|
|
swingLimit = sqrt(swingLimit2);
|
||
|
|
}
|
||
|
|
|
||
|
|
// convert into point in constraint space:
|
||
|
|
// note: twist is x-axis, swing 1 and 2 are along the z and y axes respectively
|
||
|
|
btVector3 vSwingAxis(0, xEllipse, -yEllipse);
|
||
|
|
btQuaternion qSwing(vSwingAxis, swingLimit);
|
||
|
|
btVector3 vPointInConstraintSpace(fLength,0,0);
|
||
|
|
return quatRotate(qSwing, vPointInConstraintSpace);
|
||
|
|
}
|
||
|
|
|
||
|
|
// given a twist rotation in constraint space, (pre: cone must already be removed)
|
||
|
|
// this method computes its corresponding angle and axis.
|
||
|
|
void btConeTwistConstraint::computeTwistLimitInfo(const btQuaternion& qTwist,
|
||
|
|
btScalar& twistAngle, // out
|
||
|
|
btVector3& vTwistAxis) // out
|
||
|
|
{
|
||
|
|
btQuaternion qMinTwist = qTwist;
|
||
|
|
twistAngle = qTwist.getAngle();
|
||
|
|
|
||
|
|
if (twistAngle > SIMD_PI) // long way around. flip quat and recalculate.
|
||
|
|
{
|
||
|
|
qMinTwist = operator-(qTwist);
|
||
|
|
twistAngle = qMinTwist.getAngle();
|
||
|
|
}
|
||
|
|
if (twistAngle < 0)
|
||
|
|
{
|
||
|
|
// this should never happen
|
||
|
|
int wtf = 0; wtf = wtf;
|
||
|
|
}
|
||
|
|
|
||
|
|
vTwistAxis = btVector3(qMinTwist.x(), qMinTwist.y(), qMinTwist.z());
|
||
|
|
if (twistAngle > SIMD_EPSILON)
|
||
|
|
vTwistAxis.normalize();
|
||
|
|
}
|
||
|
|
|
||
|
|
|
||
|
|
void btConeTwistConstraint::adjustSwingAxisToUseEllipseNormal(btVector3& vSwingAxis) const
|
||
|
|
{
|
||
|
|
// the swing axis is computed as the "twist-free" cone rotation,
|
||
|
|
// but the cone limit is not circular, but elliptical (if swingspan1 != swingspan2).
|
||
|
|
// so, if we're outside the limits, the closest way back inside the cone isn't
|
||
|
|
// along the vector back to the center. better (and more stable) to use the ellipse normal.
|
||
|
|
|
||
|
|
// convert swing axis to direction from center to surface of ellipse
|
||
|
|
// (ie. rotate 2D vector by PI/2)
|
||
|
|
btScalar y = -vSwingAxis.z();
|
||
|
|
btScalar z = vSwingAxis.y();
|
||
|
|
|
||
|
|
// do the math...
|
||
|
|
if (fabs(z) > SIMD_EPSILON) // avoid division by 0. and we don't need an update if z == 0.
|
||
|
|
{
|
||
|
|
// compute gradient/normal of ellipse surface at current "point"
|
||
|
|
btScalar grad = y/z;
|
||
|
|
grad *= m_swingSpan2 / m_swingSpan1;
|
||
|
|
|
||
|
|
// adjust y/z to represent normal at point (instead of vector to point)
|
||
|
|
if (y > 0)
|
||
|
|
y = fabs(grad * z);
|
||
|
|
else
|
||
|
|
y = -fabs(grad * z);
|
||
|
|
|
||
|
|
// convert ellipse direction back to swing axis
|
||
|
|
vSwingAxis.setZ(-y);
|
||
|
|
vSwingAxis.setY( z);
|
||
|
|
vSwingAxis.normalize();
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
|
||
|
|
|
||
|
|
void btConeTwistConstraint::setMotorTarget(const btQuaternion &q)
|
||
|
|
{
|
||
|
|
btTransform trACur = m_rbA.getCenterOfMassTransform();
|
||
|
|
btTransform trBCur = m_rbB.getCenterOfMassTransform();
|
||
|
|
btTransform trABCur = trBCur.inverse() * trACur;
|
||
|
|
btQuaternion qABCur = trABCur.getRotation();
|
||
|
|
btTransform trConstraintCur = (trBCur * m_rbBFrame).inverse() * (trACur * m_rbAFrame);
|
||
|
|
btQuaternion qConstraintCur = trConstraintCur.getRotation();
|
||
|
|
|
||
|
|
btQuaternion qConstraint = m_rbBFrame.getRotation().inverse() * q * m_rbAFrame.getRotation();
|
||
|
|
setMotorTargetInConstraintSpace(qConstraint);
|
||
|
|
}
|
||
|
|
|
||
|
|
|
||
|
|
void btConeTwistConstraint::setMotorTargetInConstraintSpace(const btQuaternion &q)
|
||
|
|
{
|
||
|
|
m_qTarget = q;
|
||
|
|
|
||
|
|
// clamp motor target to within limits
|
||
|
|
{
|
||
|
|
btScalar softness = 1.f;//m_limitSoftness;
|
||
|
|
|
||
|
|
// split into twist and cone
|
||
|
|
btVector3 vTwisted = quatRotate(m_qTarget, vTwist);
|
||
|
|
btQuaternion qTargetCone = shortestArcQuat(vTwist, vTwisted); qTargetCone.normalize();
|
||
|
|
btQuaternion qTargetTwist = qTargetCone.inverse() * m_qTarget; qTargetTwist.normalize();
|
||
|
|
|
||
|
|
// clamp cone
|
||
|
|
if (m_swingSpan1 >= btScalar(0.05f) && m_swingSpan2 >= btScalar(0.05f))
|
||
|
|
{
|
||
|
|
btScalar swingAngle, swingLimit; btVector3 swingAxis;
|
||
|
|
computeConeLimitInfo(qTargetCone, swingAngle, swingAxis, swingLimit);
|
||
|
|
|
||
|
|
if (fabs(swingAngle) > SIMD_EPSILON)
|
||
|
|
{
|
||
|
|
if (swingAngle > swingLimit*softness)
|
||
|
|
swingAngle = swingLimit*softness;
|
||
|
|
else if (swingAngle < -swingLimit*softness)
|
||
|
|
swingAngle = -swingLimit*softness;
|
||
|
|
qTargetCone = btQuaternion(swingAxis, swingAngle);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
// clamp twist
|
||
|
|
if (m_twistSpan >= btScalar(0.05f))
|
||
|
|
{
|
||
|
|
btScalar twistAngle; btVector3 twistAxis;
|
||
|
|
computeTwistLimitInfo(qTargetTwist, twistAngle, twistAxis);
|
||
|
|
|
||
|
|
if (fabs(twistAngle) > SIMD_EPSILON)
|
||
|
|
{
|
||
|
|
// eddy todo: limitSoftness used here???
|
||
|
|
if (twistAngle > m_twistSpan*softness)
|
||
|
|
twistAngle = m_twistSpan*softness;
|
||
|
|
else if (twistAngle < -m_twistSpan*softness)
|
||
|
|
twistAngle = -m_twistSpan*softness;
|
||
|
|
qTargetTwist = btQuaternion(twistAxis, twistAngle);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
m_qTarget = qTargetCone * qTargetTwist;
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
///override the default global value of a parameter (such as ERP or CFM), optionally provide the axis (0..5).
|
||
|
|
///If no axis is provided, it uses the default axis for this constraint.
|
||
|
|
void btConeTwistConstraint::setParam(int num, btScalar value, int axis)
|
||
|
|
{
|
||
|
|
switch(num)
|
||
|
|
{
|
||
|
|
case BT_CONSTRAINT_ERP :
|
||
|
|
case BT_CONSTRAINT_STOP_ERP :
|
||
|
|
if((axis >= 0) && (axis < 3))
|
||
|
|
{
|
||
|
|
m_linERP = value;
|
||
|
|
m_flags |= BT_CONETWIST_FLAGS_LIN_ERP;
|
||
|
|
}
|
||
|
|
else
|
||
|
|
{
|
||
|
|
m_biasFactor = value;
|
||
|
|
}
|
||
|
|
break;
|
||
|
|
case BT_CONSTRAINT_CFM :
|
||
|
|
case BT_CONSTRAINT_STOP_CFM :
|
||
|
|
if((axis >= 0) && (axis < 3))
|
||
|
|
{
|
||
|
|
m_linCFM = value;
|
||
|
|
m_flags |= BT_CONETWIST_FLAGS_LIN_CFM;
|
||
|
|
}
|
||
|
|
else
|
||
|
|
{
|
||
|
|
m_angCFM = value;
|
||
|
|
m_flags |= BT_CONETWIST_FLAGS_ANG_CFM;
|
||
|
|
}
|
||
|
|
break;
|
||
|
|
default:
|
||
|
|
btAssertConstrParams(0);
|
||
|
|
break;
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
///return the local value of parameter
|
||
|
|
btScalar btConeTwistConstraint::getParam(int num, int axis) const
|
||
|
|
{
|
||
|
|
btScalar retVal = 0;
|
||
|
|
switch(num)
|
||
|
|
{
|
||
|
|
case BT_CONSTRAINT_ERP :
|
||
|
|
case BT_CONSTRAINT_STOP_ERP :
|
||
|
|
if((axis >= 0) && (axis < 3))
|
||
|
|
{
|
||
|
|
btAssertConstrParams(m_flags & BT_CONETWIST_FLAGS_LIN_ERP);
|
||
|
|
retVal = m_linERP;
|
||
|
|
}
|
||
|
|
else if((axis >= 3) && (axis < 6))
|
||
|
|
{
|
||
|
|
retVal = m_biasFactor;
|
||
|
|
}
|
||
|
|
else
|
||
|
|
{
|
||
|
|
btAssertConstrParams(0);
|
||
|
|
}
|
||
|
|
break;
|
||
|
|
case BT_CONSTRAINT_CFM :
|
||
|
|
case BT_CONSTRAINT_STOP_CFM :
|
||
|
|
if((axis >= 0) && (axis < 3))
|
||
|
|
{
|
||
|
|
btAssertConstrParams(m_flags & BT_CONETWIST_FLAGS_LIN_CFM);
|
||
|
|
retVal = m_linCFM;
|
||
|
|
}
|
||
|
|
else if((axis >= 3) && (axis < 6))
|
||
|
|
{
|
||
|
|
btAssertConstrParams(m_flags & BT_CONETWIST_FLAGS_ANG_CFM);
|
||
|
|
retVal = m_angCFM;
|
||
|
|
}
|
||
|
|
else
|
||
|
|
{
|
||
|
|
btAssertConstrParams(0);
|
||
|
|
}
|
||
|
|
break;
|
||
|
|
default :
|
||
|
|
btAssertConstrParams(0);
|
||
|
|
}
|
||
|
|
return retVal;
|
||
|
|
}
|
||
|
|
|
||
|
|
|
||
|
|
|