Commit fdab27a6 authored by Sebastian Eibl's avatar Sebastian Eibl
Browse files

[API] made setMass protected

static calcMass and calcInertia functions are now available for all geometries
parent 7d572372
......@@ -1693,6 +1693,25 @@ inline const Matrix3<Type> fabs( const Matrix3<Type>& m )
//**********************************************************************************************************************
//**********************************************************************************************************************
/*!\fn bool math::isinf( const Matrix3<Type>& m )
// \brief Checks the given matrix for infinite elements.
//
// \param m The matrix to be checked for infinite elements.
// \return \a true if at least one element of the matrix is infinite, \a false otherwise.
*/
template< typename Type >
inline bool isinf( const Matrix3<Type>& m )
{
if( math::isinf( m[0] ) || math::isinf( m[1] ) || math::isinf( m[2] ) ||
math::isinf( m[3] ) || math::isinf( m[4] ) || math::isinf( m[5] ) ||
math::isinf( m[6] ) || math::isinf( m[7] ) || math::isinf( m[8] ) )
return true;
else return false;
}
//**********************************************************************************************************************
//**********************************************************************************************************************
// \brief Calculates the tensor product of two vectors.
//**********************************************************************************************************************
......
......@@ -90,14 +90,15 @@ void ConvexPolyhedron::init( const Vec3& gpos, const Vec3& rpos, const Quat& q,
q_ = q; // Setting the orientation
R_ = q_.toRotationMatrix(); // Setting the rotation matrix
// Calculating the ConvexPolyhedron mass
mass_ = getVolume() * Material::getDensity( getMaterial() );
// Calculating the moment of inertia
calcInertia();
setGlobal( global );
setMass( infiniteMass ); // sets inverse mass and interatio tensor
if (infiniteMass)
{
setMassAndInertiaToInfinity();
} else
{
// sets inverse mass and interatio tensor
setMassAndInertia( getVolume() * Material::getDensity( getMaterial() ), mesh::computeIntertiaTensor( mesh_ ) );
}
setCommunicating( communicating );
setFinite( true );
......@@ -159,17 +160,6 @@ void ConvexPolyhedron::calcBoundingBox()
}
//*************************************************************************************************
//*************************************************************************************************
/*!\brief Calculation of the moment of inertia in reference to the body frame of the ConvexPolyhedron.
*
* \return void
*/
void ConvexPolyhedron::calcInertia()
{
I_ = mesh::computeIntertiaTensor( mesh_ );
}
//*************************************************************************************************
//*************************************************************************************************
/*!\brief Calculation of the volume of the ConvexPolyhedron.
......
......@@ -127,7 +127,6 @@ protected:
/*!\name Utility functions */
//@{
virtual void calcBoundingBox(); // Calculation of the axis-aligned bounding box
void calcInertia(); // Calculation of the moment of inertia
virtual TriangleMesh::VertexHandle supportVertex( const TriangleMesh::Normal & d, const TriangleMesh::VertexHandle startVertex ) const;
//@}
//**********************************************************************************************
......
......@@ -288,7 +288,7 @@ inline void HardContactSemiImplicitTimesteppingSolvers::timestep( const real_t d
for( auto body = globalBodyStorage_->begin(); body != globalBodyStorage_->end(); ++body, ++j ) {
body->wake(); // BUGFIX: Force awaking of all bodies!
body->index_ = j;
WALBERLA_ASSERT( body->hasInfiniteMass(), "fatal" );
WALBERLA_CHECK( body->hasInfiniteMass(), "Global bodies need to have infinite mass as they are not communicated!" );
initializeVelocityCorrections( *body, bodyCache.dv_[j], bodyCache.dw_[j], dt ); // use applied external forces to calculate starting velocity
......@@ -347,7 +347,7 @@ inline void HardContactSemiImplicitTimesteppingSolvers::timestep( const real_t d
for( auto body = globalBodyStorage_->begin(); body != globalBodyStorage_->end(); ++body, ++j ) {
body->wake(); // BUGFIX: Force awaking of all bodies!
body->index_ = j;
WALBERLA_ASSERT( body->hasInfiniteMass(), "fatal" );
WALBERLA_CHECK( body->hasInfiniteMass(), "Global bodies need to have infinite mass as they are not communicated!" );
initializeVelocityCorrections( *body, bodyCache.dv_[j], bodyCache.dw_[j], dt ); // use applied external forces to calculate starting velocity
......
......@@ -75,15 +75,15 @@ Box::Box( id_t sid, id_t uid, const Vec3& gpos, const Vec3& rpos, const Quat& q,
q_ = q; // Setting the orientation
R_ = q_.toRotationMatrix(); // Setting the rotation matrix
// Calculating the sphere mass
mass_ = calcMass( lengths, Material::getDensity( material ) );
invMass_ = real_c(1) / mass_;
// Calculating the moment of inertia
calcInertia();
setGlobal( global );
setMass( infiniteMass );
if (infiniteMass)
{
setMassAndInertiaToInfinity();
} else
{
auto mass = calcMass( lengths, Material::getDensity( material ) );
setMassAndInertia( mass, calcInertia( lengths, mass ) );
}
setCommunicating( communicating );
setFinite( true );
......@@ -130,8 +130,8 @@ Box::~Box()
bool Box::containsRelPointImpl( real_t px, real_t py, real_t pz ) const
{
return std::fabs(px) <= real_t(0.5)*lengths_[0] &&
std::fabs(py) <= real_t(0.5)*lengths_[1] &&
std::fabs(pz) <= real_t(0.5)*lengths_[2];
std::fabs(py) <= real_t(0.5)*lengths_[1] &&
std::fabs(pz) <= real_t(0.5)*lengths_[2];
}
//*************************************************************************************************
......@@ -307,7 +307,7 @@ void Box::calcBoundingBox()
const real_t ylength( real_t(0.5) * ( fabs(R_[3]*lengths_[0]) + fabs(R_[4]*lengths_[1]) + fabs(R_[5]*lengths_[2]) ) + contactThreshold );
const real_t zlength( real_t(0.5) * ( fabs(R_[6]*lengths_[0]) + fabs(R_[7]*lengths_[1]) + fabs(R_[8]*lengths_[2]) ) + contactThreshold );
aabb_ = math::AABB(
gpos_[0] - xlength,
gpos_[0] - xlength,
gpos_[1] - ylength,
gpos_[2] - zlength,
gpos_[0] + xlength,
......@@ -326,12 +326,12 @@ void Box::calcBoundingBox()
*
* \return void
*/
void Box::calcInertia()
Mat3 Box::calcInertia(const Vec3& length, const real_t mass)
{
I_[0] = mass_/static_cast<real_t>( 12 ) * ( lengths_[1]*lengths_[1] + lengths_[2]*lengths_[2] );
I_[4] = mass_/static_cast<real_t>( 12 ) * ( lengths_[0]*lengths_[0] + lengths_[2]*lengths_[2] );
I_[8] = mass_/static_cast<real_t>( 12 ) * ( lengths_[0]*lengths_[0] + lengths_[1]*lengths_[1] );
Iinv_ = I_.getInverse();
return Mat3::makeDiagonalMatrix(
mass/static_cast<real_t>( 12 ) * ( length[1]*length[1] + length[2]*length[2] ),
mass/static_cast<real_t>( 12 ) * ( length[0]*length[0] + length[2]*length[2] ),
mass/static_cast<real_t>( 12 ) * ( length[0]*length[0] + length[1]*length[1] ));
}
//*************************************************************************************************
......
......@@ -96,6 +96,7 @@ public:
static inline real_t calcMass( const Vec3& l, real_t density );
static inline real_t calcDensity( real_t x, real_t y, real_t z, real_t mass );
static inline real_t calcDensity( const Vec3& l, real_t mass );
static inline Mat3 calcInertia( const Vec3& length, const real_t mass); // Calculation of the moment of inertia
//@}
//**********************************************************************************************
......@@ -143,7 +144,6 @@ protected:
/*!\name Utility functions */
//@{
virtual void calcBoundingBox(); // Calculation of the axis-aligned bounding box
void calcInertia(); // Calculation of the moment of inertia
//@}
//**********************************************************************************************
......
......@@ -75,15 +75,15 @@ Capsule::Capsule( id_t sid, id_t uid, const Vec3& gpos, const Vec3& rpos, const
q_ = q; // Setting the orientation
R_ = q_.toRotationMatrix(); // Setting the rotation matrix
// Calculating the capsule mass
mass_ = radius*radius * math::M_PI * ( real_t(4)/real_t(3) * radius + length ) * Material::getDensity( material );
invMass_ = real_t(1) / mass_;
// Calculating the moment of inertia
calcInertia();
setGlobal( global );
setMass( infiniteMass );
if (infiniteMass)
{
setMassAndInertiaToInfinity();
} else
{
auto mass = calcMass( radius, length, Material::getDensity( material ) );
setMassAndInertia( mass, calcInertia( radius, length, Material::getDensity( material ) ) );
}
setCommunicating( communicating );
setFinite( true );
......@@ -192,29 +192,25 @@ void Capsule::calcBoundingBox()
*
* \return void
*/
void Capsule::calcInertia()
Mat3 Capsule::calcInertia( const real_t radius, const real_t length, const real_t density)
{
const real_t density( calcDensity( radius_, length_, mass_ ) );
const real_t sphereMass( real_t (4)/real_t (3) * math::M_PI * radius_*radius_*radius_ * density );
const real_t cylinderMass( math::M_PI * radius_*radius_ * length_ * density );
const real_t sphereMass( real_t (4)/real_t (3) * math::M_PI * radius*radius*radius * density );
const real_t cylinderMass( math::M_PI * radius*radius * length * density );
// 'Ia' represent the moment of inertia along the x-axis. 'Ia' contains the following two parts:
// - cylinder : I = (1/2)*mass*radius^2
// - sphere : I = (2/5)*mass*radius^2
const real_t Ia( radius_*radius_ * ( real_t (0.5)*cylinderMass + real_t (0.4)*sphereMass ) );
const real_t Ia( radius*radius * ( real_t (0.5)*cylinderMass + real_t (0.4)*sphereMass ) );
// 'Ib' represent the moment of inertia along the y- and z-axis. 'Ib' contains the following two parts,
// where full_length is the full length of the cylinder part and half_length is (1/2)*full_length:
// - cylinder : I = mass*( (1/4)*radius^2 + (1/12)*full_length^2 )
// - sphere : I = mass*( (2/5)*radius^2 + half_length^2 + (3/4)*half_length*radius )
const real_t Ib( cylinderMass*( real_t (0.25)*radius_*radius_ + real_t (1)/real_t (12)*length_*length_ ) +
sphereMass*( real_t (0.4)*radius_*radius_ + real_t (0.375)*radius_*length_ + real_t (0.25)*length_*length_ ) );
const real_t Ib( cylinderMass*( real_t (0.25)*radius*radius + real_t (1)/real_t (12)*length*length ) +
sphereMass*( real_t (0.4)*radius*radius + real_t (0.375)*radius*length + real_t (0.25)*length*length ) );
// Setting the moment of inertia (capsule is aligned along the x-axis)
I_[0] = Ia;
I_[4] = Ib;
I_[8] = Ib;
Iinv_ = I_.getInverse();
return Mat3::makeDiagonalMatrix(Ia, Ib, Ib);
}
//*************************************************************************************************
......
......@@ -106,6 +106,7 @@ public:
static inline real_t calcVolume( real_t radius, real_t length );
static inline real_t calcMass( real_t radius, real_t length, real_t density );
static inline real_t calcDensity( real_t radius, real_t length, real_t mass );
static Mat3 calcInertia( const real_t radius, const real_t length, const real_t density);
//@}
//**********************************************************************************************
......@@ -139,7 +140,6 @@ protected:
/*!\name Utility functions */
//@{
virtual void calcBoundingBox(); // Calculation of the axis-aligned bounding box
void calcInertia(); // Calculation of the moment of inertia
//@}
//**********************************************************************************************
......
......@@ -61,7 +61,7 @@ CylindricalBoundary::CylindricalBoundary( id_t sid, id_t uid, const Vec3& gpos,
//boundaries are always considered locally and have infinite mass
setGlobal( true );
setCommunicating( false );
setMass( true );
setMassAndInertiaToInfinity();
setFinite( false );
// Checking the radius
......
......@@ -84,15 +84,15 @@ Ellipsoid::Ellipsoid( id_t const typeId, id_t sid, id_t uid, const Vec3& gpos, c
q_ = q; // Setting the orientation
R_ = q_.toRotationMatrix(); // Setting the rotation matrix
// Calculating the Ellipsoid mass
mass_ = calcMass(semiAxes_, Material::getDensity( material ));
invMass_ = real_c(1) / mass_;
// Calculating the moment of inertia
calcInertia();
setGlobal( global );
setMass( infiniteMass );
if (infiniteMass)
{
setMassAndInertiaToInfinity();
} else
{
auto mass = calcMass( semiAxes, Material::getDensity( material ) );
setMassAndInertia( mass, calcInertia( mass, semiAxes ) );
}
setCommunicating( communicating );
setFinite( true );
......
......@@ -55,8 +55,8 @@ namespace pe {
* \brief Base class for the Ellipsoid geometry.
*
* The Ellipsoid class represents the base class for the Ellipsoid geometry. It provides
* the basic functionality of a Ellipsoid. An Ellipsoid is obtained from a sphere by deforming it by means
* of directional scalings. Its three semi-axes corresponding to the x, y, z direction are labeled
* the basic functionality of a Ellipsoid. An Ellipsoid is obtained from a sphere by deforming it by means
* of directional scalings. Its three semi-axes corresponding to the x, y, z direction are labeled
* a, b, c.
*/
class Ellipsoid : public GeomPrimitive
......@@ -66,11 +66,11 @@ public:
/*!\name Constructors */
//@{
explicit Ellipsoid( id_t sid, id_t uid, const Vec3& gpos, const Vec3& rpos, const Quat& q,
const Vec3& semiAxes, MaterialID material,
const bool global, const bool communicating, const bool infiniteMass );
const Vec3& semiAxes, MaterialID material,
const bool global, const bool communicating, const bool infiniteMass );
explicit Ellipsoid( id_t const typeID, id_t sid, id_t uid, const Vec3& gpos, const Vec3& rpos, const Quat& q,
const Vec3& semiAxes, MaterialID material,
const bool global, const bool communicating, const bool infiniteMass );
const Vec3& semiAxes, MaterialID material,
const bool global, const bool communicating, const bool infiniteMass );
//@}
//**********************************************************************************************
......@@ -116,6 +116,7 @@ public:
static inline real_t calcVolume( const Vec3& semiAxes );
static inline real_t calcMass( const Vec3& semiAxes, real_t density );
static inline real_t calcDensity( const Vec3& semiAxes, real_t mass );
static inline Mat3 calcInertia( const real_t mass, const Vec3& semiAxes );
//@}
//**********************************************************************************************
......@@ -132,7 +133,6 @@ protected:
/*!\name Utility functions */
//@{
inline virtual void calcBoundingBox(); // Calculation of the axis-aligned bounding box
inline void calcInertia(); // Calculation of the moment of inertia
//@}
//**********************************************************************************************
......@@ -140,7 +140,7 @@ protected:
/*!\name Member variables */
//@{
Vec3 semiAxes_; //!< Semi-axes of the Ellipsoid.
/*!< The radius is constrained to values larger than 0.0. */
/*!< The radius is constrained to values larger than 0.0. */
//@}
//**********************************************************************************************
private:
......@@ -264,14 +264,14 @@ inline void Ellipsoid::calcBoundingBox()
const real_t ylength( fabs(R_[3]*semiAxes_[0]) + fabs(R_[4]*semiAxes_[1]) + fabs(R_[5]*semiAxes_[2]) + contactThreshold );
const real_t zlength( fabs(R_[6]*semiAxes_[0]) + fabs(R_[7]*semiAxes_[1]) + fabs(R_[8]*semiAxes_[2]) + contactThreshold );
aabb_ = math::AABB(
gpos_[0] - xlength,
gpos_[0] - xlength,
gpos_[1] - ylength,
gpos_[2] - zlength,
gpos_[0] + xlength,
gpos_[1] + ylength,
gpos_[2] + zlength
);
// WALBERLA_ASSERT( aabb_.isValid() , "Invalid bounding box detected" );
// WALBERLA_ASSERT( aabb_.isValid() , "Invalid bounding box detected" );
WALBERLA_ASSERT( aabb_.contains( gpos_ ), "Invalid bounding box detected("<< getSystemID() <<")\n" << "pos: " << gpos_ << "\nlengths: " << xlength << "," << ylength << ", " << zlength<< "\nvel: " << getLinearVel() << "\nbox: " << aabb_ );
}
//*************************************************************************************************
......@@ -282,12 +282,12 @@ inline void Ellipsoid::calcBoundingBox()
*
* \return void
*/
inline void Ellipsoid::calcInertia()
inline Mat3 Ellipsoid::calcInertia( const real_t mass, const Vec3& semiAxes )
{
I_[0] = real_c(0.2) * mass_ * (semiAxes_[1] * semiAxes_[1] + semiAxes_[2] * semiAxes_[2]);
I_[4] = real_c(0.2) * mass_ * (semiAxes_[2] * semiAxes_[2] + semiAxes_[0] * semiAxes_[0]);
I_[8] = real_c(0.2) * mass_ * (semiAxes_[0] * semiAxes_[0] + semiAxes_[1] * semiAxes_[1]);
Iinv_ = I_.getInverse();
return Mat3::makeDiagonalMatrix(
real_c(0.2) * mass * (semiAxes[1] * semiAxes[1] + semiAxes[2] * semiAxes[2]),
real_c(0.2) * mass * (semiAxes[2] * semiAxes[2] + semiAxes[0] * semiAxes[0]),
real_c(0.2) * mass * (semiAxes[0] * semiAxes[0] + semiAxes[1] * semiAxes[1]));
}
//*************************************************************************************************
......@@ -303,11 +303,11 @@ inline Vec3 Ellipsoid::support( const Vec3& d ) const
auto len = d.sqrLength();
if (!math::equal(len, real_t(0)))
{
Vec3 d_loc = vectorFromWFtoBF(d);
Vec3 norm_vec(d_loc[0] * semiAxes_[0], d_loc[1] * semiAxes_[1], d_loc[2] * semiAxes_[2]);
real_t norm_length = norm_vec.length();
Vec3 local_support = (real_t(1.0) / norm_length) * Vec3(semiAxes_[0] * semiAxes_[0] * d_loc[0],
semiAxes_[1] * semiAxes_[1] * d_loc[1], semiAxes_[2] * semiAxes_[2] * d_loc[2]);
Vec3 d_loc = vectorFromWFtoBF(d);
Vec3 norm_vec(d_loc[0] * semiAxes_[0], d_loc[1] * semiAxes_[1], d_loc[2] * semiAxes_[2]);
real_t norm_length = norm_vec.length();
Vec3 local_support = (real_t(1.0) / norm_length) * Vec3(semiAxes_[0] * semiAxes_[0] * d_loc[0],
semiAxes_[1] * semiAxes_[1] * d_loc[1], semiAxes_[2] * semiAxes_[2] * d_loc[2]);
return pointFromBFtoWF(local_support);
} else
{
......
......@@ -71,14 +71,14 @@ Plane::Plane( id_t sid, id_t uid,
//planes are always considered locally and have infinite mass
setGlobal( true );
setCommunicating( false );
setMass( true );
setMassAndInertiaToInfinity();
setFinite( false );
// Checking the mass properties of the plane
WALBERLA_ASSERT_FLOAT_EQUAL( mass_ , real_c(0), "Mass of plane is not 0" );
WALBERLA_ASSERT ( std::isinf(mass_), "Mass of plane is not infinity" );
WALBERLA_ASSERT_FLOAT_EQUAL( invMass_ , real_c(0), "Inverse mass of plane is not 0" );
WALBERLA_ASSERT( I_ == Mat3(real_c(0)), "Moment of inertia of plane is not 0\n" << I_ );
WALBERLA_ASSERT( Iinv_ == Mat3(real_c(0)), "Inverse moment of inertia of plane is not 0" );
WALBERLA_ASSERT ( isinf(I_), "Moment of inertia of plane is not infinity\n" << I_ );
WALBERLA_ASSERT_FLOAT_EQUAL( Iinv_ , Mat3(real_c(0)), "Inverse moment of inertia of plane is not 0" );
// Checking the plane normal
// Since the plane constructor is never directly called but only used in a small number
......
......@@ -46,9 +46,6 @@ PlaneID createPlane( BodyStorage& globalStorage, id_t uid, Vec3 normal, const Ve
PlaneID plane = new Plane( sid, uid, gpos, normal, normal*gpos, material );
plane->setCommunicating( false );
plane->setMass( true );
globalStorage.add(plane);
// Logging the successful creation of the plane
......
......@@ -166,12 +166,12 @@ public:
inline void resetSB() ;
inline void setFinite ( const bool finite );
inline void setMass ( bool infinite );
inline void setVisible ( bool visible );
inline void setPosition ( real_t px, real_t py, real_t pz );
inline void setPosition ( const Vec3& gpos );
inline void setOrientation( real_t r, real_t i, real_t j, real_t k );
inline void setOrientation( const Quat& q );
inline void setMassAndInertiaToInfinity();
inline void setRelLinearVel ( real_t vx, real_t vy, real_t vz );
inline void setRelLinearVel ( const Vec3& lvel );
......@@ -357,8 +357,10 @@ protected:
virtual void rotateImpl ( const Quat& dq );
virtual void rotateAroundOriginImpl( const Quat& dq );
virtual void rotateAroundPointImpl ( const Vec3& point, const Quat& dq );
virtual bool containsRelPointImpl ( real_t px, real_t py, real_t pz ) const;
virtual bool isSurfaceRelPointImpl ( real_t px, real_t py, real_t pz ) const;
virtual bool containsRelPointImpl ( real_t px, real_t py, real_t pz ) const;
virtual bool isSurfaceRelPointImpl ( real_t px, real_t py, real_t pz ) const;
inline void setMassAndInertia ( const real_t mass, const Mat3& inertia );
inline void calcRelPosition();
//@}
......@@ -706,7 +708,7 @@ inline bool RigidBody::isFixed() const
*/
inline bool RigidBody::hasInfiniteMass() const
{
return (isIdentical(invMass_, real_t(0))) && ( Iinv_.isZero() );
return std::isinf(getMass());
}
//*************************************************************************************************
......@@ -2313,39 +2315,30 @@ inline void RigidBody::setFinite( const bool finite )
//*************************************************************************************************
/*!\brief Setting the global position (the center of mass) of the rigid body fixed.
*
* \return void
*
* This function either fixes the global position (the center of mass) of a finite
* rigid body. If the body is contained in a superordinate body, fixing the contained
* body will also fix the global position of the superordinate body.
* Fixation is permanent and a fixed body cannot be unfixed anymore.
/*!\brief Sets mass and inertia of a rigid body. Also calculates inverse values.
*
* In case of a <b>MPI parallel simulation</b>, changing the settings of a (local) rigid body
* on one process may invalidate the settings of the rigid body on another process. In order to
* synchronize all rigid bodies after local changes, the simulation has to be synchronized
* by the user. Note that any changes on remote rigid
* bodies are neglected and overwritten by the settings of the rigid body on its local process!
* \param mass mass to be set (may be infinity)
* \param inertia inertia to be set (if mass is infinity this one will be ignored)
*/
inline void RigidBody::setMass( bool infinite )
inline void RigidBody::setMassAndInertia( const real_t mass, const Mat3& inertia )
{
if (infinite)
if ( std::isinf( mass ) )
{
// Adjusting the inverse mass and inverse moment of inertia
mass_ = std::numeric_limits<real_t>::infinity();
I_ = Mat3::makeDiagonalMatrix(std::numeric_limits<real_t>::infinity());
invMass_ = real_c(0);
Iinv_ = Mat3( real_c(0) );
// Setting the linear and angular velocity to zero
v_ = Vec3(0);
w_ = Vec3(0);
} else
{
// Adjusting the inverse mass and inverse moment of inertia
mass_ = mass;
I_ = inertia;
invMass_ = real_c(1) / mass_;
Iinv_ = I_.getInverse();
}
if (hasSuperBody()) WALBERLA_LOG_WARNING("Changing the mass of a body contained in a Union is currently untested!!!");
// Signaling the fixation change to the superordinate body
signalFixation();
}
......@@ -2555,6 +2548,19 @@ inline void RigidBody::setOrientation( const Quat& q )
}
//*************************************************************************************************
//*************************************************************************************************
/*!\fn void RigidBody::setMassAndInertiaToInfinity( )
* \brief Setting the mass to infinity. This will also make the inertia tensor infinite.
*
* \attention This cannot be undone!
*/
//*************************************************************************************************
inline void RigidBody::setMassAndInertiaToInfinity()
{
setMassAndInertia( std::numeric_limits<real_t>::infinity(), Mat3(std::numeric_limits<real_t>::infinity()) );
}
//*************************************************************************************************
//=================================================================================================
......@@ -3261,7 +3267,7 @@ inline Vec3 RigidBody::supportContactThreshold(const Vec3& /*d*/) const
inline void RigidBody::fix()
{
setCommunicating( false );
setMass( true );
setMassAndInertia( getMass(), getInertia() );
// Signaling the fixation change to the superordinate body
signalFixation();
......
......@@ -83,15 +83,15 @@ Sphere::Sphere( id_t const typeId, id_t sid, id_t uid, const Vec3& gpos, const V
q_ = q; // Setting the orientation
R_ = q_.toRotationMatrix(); // Setting the rotation matrix
// Calculating the sphere mass
mass_ = calcMass( radius, Material::getDensity( material ) );
invMass_ = real_c(1) / mass_;
// Calculating the moment of inertia
calcInertia();
setGlobal( global );
setMass( infiniteMass );
if (infiniteMass)
{
setMassAndInertia( std::numeric_limits<real_t>::infinity(), Mat3(real_t(0)) );
} else
{
auto mass = calcMass( radius, Material::getDensity( material ) );
setMassAndInertia( mass, calcInertia( radius, mass ) );
}
setCommunicating( communicating );
setFinite( true );
......@@ -194,20 +194,20 @@ void Sphere::print( std::ostream& os, const char* tab ) const
<< tab << " Linear velocity = " << getLinearVel() << "\n"
<< tab << " Angular velocity = " << getAngularVel() << "\n";
// if( verboseMode )
// {
os << tab << " Bounding box = " << getAABB() << "\n"
<< tab << " Quaternion = " << getQuaternion() << "\n"
<< tab << " Rotation matrix = ( " << setw(9) << R_[0] << " , " << setw(9) << R_[1] << " , " << setw(9) << R_[2] << " )\n"
<< tab << " ( " << setw(9) << R_[3] << " , " << setw(9) << R_[4] << " , " << setw(9) << R_[5] << " )\n"
<< tab << " ( " << setw(9) << R_[6] << " , " << setw(9) << R_[7] << " , " << setw(9) << R_[8] << " )\n";
os << std::setiosflags(std::ios::right)
<< tab << " Moment of inertia = ( " << setw(9) << I_[0] << " , " << setw(9) << I_[1] << " , " << setw(9) << I_[2] << " )\n"
<< tab << " ( " << setw(9) << I_[3] << " , " << setw(9) << I_[4] << " , " << setw(9) << I_[5] << " )\n"
<< tab << " ( " << setw(9) << I_[6] << " , " << setw(9) << I_[7] << " , " << setw(9) << I_[8] << " )\n"
<< std::resetiosflags(std::ios::right);
// }
// if( verboseMode )
// {
os << tab << " Bounding box = " << getAABB() << "\n"
<< tab << " Quaternion = " << getQuaternion() << "\n"
<< tab << " Rotation matrix = ( " << setw(9) << R_[0] << " , " << setw(9) << R_[1] << " , " << setw(9) << R_[2] << " )\n"
<< tab << " ( " << setw(9) << R_[3] << " , " << setw(9) << R_[4] << " , " << setw(9) <<