| 
 | 1 | +#include "TimeStepRigidBodyModel.h"  | 
 | 2 | +#include "Demos/Utils/TimeManager.h"  | 
 | 3 | +#include "PositionBasedDynamics/PositionBasedDynamics.h"  | 
 | 4 | +#include "PositionBasedDynamics/TimeIntegration.h"  | 
 | 5 | + | 
 | 6 | +using namespace PBD;  | 
 | 7 | +using namespace std;  | 
 | 8 | + | 
 | 9 | +TimeStepRigidBodyModel::TimeStepRigidBodyModel()  | 
 | 10 | +{  | 
 | 11 | +	m_velocityUpdateMethod = 0;  | 
 | 12 | +}  | 
 | 13 | + | 
 | 14 | +TimeStepRigidBodyModel::~TimeStepRigidBodyModel(void)  | 
 | 15 | +{  | 
 | 16 | +}  | 
 | 17 | + | 
 | 18 | +void TimeStepRigidBodyModel::step(RigidBodyModel &model)  | 
 | 19 | +{  | 
 | 20 | + 	TimeManager *tm = TimeManager::getCurrent ();  | 
 | 21 | + 	const float h = tm->getTimeStepSize();  | 
 | 22 | +   | 
 | 23 | + 	clearAccelerations(model);  | 
 | 24 | +	RigidBodyModel::RigidBodyVector &rb = model.getRigidBodies();  | 
 | 25 | +	for (size_t i = 0; i < rb.size(); i++)  | 
 | 26 | + 	{   | 
 | 27 | +		rb[i].getLastPosition() = rb[i].getOldPosition();  | 
 | 28 | +		rb[i].getOldPosition() = rb[i].getPosition();  | 
 | 29 | +		TimeIntegration::semiImplicitEuler(h, rb[i].getMass(), rb[i].getPosition(), rb[i].getVelocity(), rb[i].getAcceleration());  | 
 | 30 | +		rb[i].getLastRotation() = rb[i].getOldRotation();  | 
 | 31 | +		rb[i].getOldRotation() = rb[i].getRotation();  | 
 | 32 | +		TimeIntegration::semiImplicitEulerRotation(h, rb[i].getMass(), rb[i].getInertiaTensorInverseW(), rb[i].getRotation(), rb[i].getAngularVelocity(), rb[i].getTorque());  | 
 | 33 | +		rb[i].rotationUpdated();  | 
 | 34 | + 	}  | 
 | 35 | +   | 
 | 36 | + 	constraintProjection(model);  | 
 | 37 | +   | 
 | 38 | + 	// Update velocities	  | 
 | 39 | +	for (size_t i = 0; i < rb.size(); i++)  | 
 | 40 | + 	{  | 
 | 41 | +		if (m_velocityUpdateMethod == 0)  | 
 | 42 | +		{  | 
 | 43 | +			TimeIntegration::velocityUpdateFirstOrder(h, rb[i].getMass(), rb[i].getPosition(), rb[i].getOldPosition(), rb[i].getVelocity());  | 
 | 44 | +			TimeIntegration::angularVelocityUpdateFirstOrder(h, rb[i].getMass(), rb[i].getRotation(), rb[i].getOldRotation(), rb[i].getAngularVelocity());  | 
 | 45 | +		}  | 
 | 46 | +		else  | 
 | 47 | +		{  | 
 | 48 | +			TimeIntegration::velocityUpdateSecondOrder(h, rb[i].getMass(), rb[i].getPosition(), rb[i].getOldPosition(), rb[i].getLastPosition(), rb[i].getVelocity());  | 
 | 49 | +			TimeIntegration::angularVelocityUpdateSecondOrder(h, rb[i].getMass(), rb[i].getRotation(), rb[i].getOldRotation(), rb[i].getLastRotation(), rb[i].getAngularVelocity());  | 
 | 50 | +		}  | 
 | 51 | + 	}  | 
 | 52 | + | 
 | 53 | +	// compute new time	  | 
 | 54 | +	tm->setTime (tm->getTime () + h);  | 
 | 55 | +}  | 
 | 56 | + | 
 | 57 | +/** Clear accelerations and add gravitation.  | 
 | 58 | + */  | 
 | 59 | +void TimeStepRigidBodyModel::clearAccelerations(RigidBodyModel &model)  | 
 | 60 | +{  | 
 | 61 | +	RigidBodyModel::RigidBodyVector &rb = model.getRigidBodies();  | 
 | 62 | + 	const Eigen::Vector3f grav(0.0f, -9.81f, 0.0f);  | 
 | 63 | + 	for (size_t i=0; i < rb.size(); i++)  | 
 | 64 | + 	{  | 
 | 65 | + 		// Clear accelerations of dynamic particles  | 
 | 66 | + 		if (rb[i].getMass() != 0.0)  | 
 | 67 | + 		{  | 
 | 68 | +			Eigen::Vector3f &a = rb[i].getAcceleration();  | 
 | 69 | + 			a = grav;  | 
 | 70 | + 		}  | 
 | 71 | + 	}  | 
 | 72 | +}  | 
 | 73 | + | 
 | 74 | +void TimeStepRigidBodyModel::reset(RigidBodyModel &model)  | 
 | 75 | +{  | 
 | 76 | + | 
 | 77 | +}  | 
 | 78 | + | 
 | 79 | +void TimeStepRigidBodyModel::constraintProjection(RigidBodyModel &model)  | 
 | 80 | +{  | 
 | 81 | +	const unsigned int maxIter = 5;  | 
 | 82 | +	unsigned int iter = 0;  | 
 | 83 | + | 
 | 84 | +	RigidBodyModel::RigidBodyVector &rb = model.getRigidBodies();  | 
 | 85 | +	RigidBodyModel::BallJointVector &bj = model.getBallJoints();  | 
 | 86 | + | 
 | 87 | + 	while (iter < maxIter)  | 
 | 88 | + 	{  | 
 | 89 | +		for (unsigned int i = 0; i < bj.size(); i++)  | 
 | 90 | +		{  | 
 | 91 | +			model.updateBallJoint(i);  | 
 | 92 | + | 
 | 93 | +			RigidBody &rb1 = rb[bj[i].m_index[0]];  | 
 | 94 | +			RigidBody &rb2 = rb[bj[i].m_index[1]];  | 
 | 95 | + | 
 | 96 | +			Eigen::Vector3f corr_x1, corr_x2;  | 
 | 97 | +			Eigen::Quaternionf corr_q1, corr_q2;  | 
 | 98 | +			const bool res = PositionBasedDynamics::solveRigidBodyBallJoint(  | 
 | 99 | +				bj[i].m_points[0],  | 
 | 100 | +				rb1.getMass(),  | 
 | 101 | +				rb1.getPosition(),  | 
 | 102 | +				rb1.getInertiaTensorInverseW(),  | 
 | 103 | +				rb1.getRotation(),  | 
 | 104 | +				bj[i].m_points[1],  | 
 | 105 | +				rb2.getMass(),  | 
 | 106 | +				rb2.getPosition(),  | 
 | 107 | +				rb2.getInertiaTensorInverseW(),  | 
 | 108 | +				rb2.getRotation(),  | 
 | 109 | +				corr_x1,  | 
 | 110 | +				corr_q1,  | 
 | 111 | +				corr_x2,  | 
 | 112 | +				corr_q2);  | 
 | 113 | + | 
 | 114 | +			if (res)  | 
 | 115 | +			{  | 
 | 116 | +				if (rb1.getMass() != 0.0f)  | 
 | 117 | +				{  | 
 | 118 | +					rb1.getPosition() += corr_x1;  | 
 | 119 | +					rb1.getRotation().coeffs() += corr_q1.coeffs();  | 
 | 120 | +					rb1.getRotation().normalize();  | 
 | 121 | +					rb1.rotationUpdated();  | 
 | 122 | +				}  | 
 | 123 | +				if (rb2.getMass() != 0.0f)  | 
 | 124 | +				{  | 
 | 125 | +					rb2.getPosition() += corr_x2;  | 
 | 126 | +					rb2.getRotation().coeffs() += corr_q2.coeffs();  | 
 | 127 | +					rb2.getRotation().normalize();  | 
 | 128 | +					rb2.rotationUpdated();  | 
 | 129 | +				}  | 
 | 130 | +			}  | 
 | 131 | +		}  | 
 | 132 | + 		iter++;  | 
 | 133 | + 	}  | 
 | 134 | +}  | 
0 commit comments