Commit f0a69f7d authored by Tobias Leemann's avatar Tobias Leemann Committed by Sebastian Eibl

[API] Union redesign

[API] removed rpos from all body constructors

Changes: Bodies that are part of a union now store all attributes relative to its respective coordinate frame. This affects gpos and q (Rotation), which are implicitly relative once a body becomes part of the union. These properties became private and the global Position can be accessed via getGlobalPos() which performs a recursive calcuation of the position of this body in the world frame. This avoids error-prone update procedures, if a Union is tranlated or rotated.
parent b1cecb31
......@@ -74,7 +74,7 @@ inline mesh::pe::ConvexPolyhedronPtr instantiate( mpi::RecvBuffer& buffer, const
mesh::pe::ConvexPolyhedronParameters subobjparam;
unmarshal( buffer, subobjparam );
correctBodyPosition(domain, block.center(), subobjparam.gpos_);
auto cp = std::make_unique<mesh::pe::ConvexPolyhedron>( subobjparam.sid_, subobjparam.uid_, subobjparam.gpos_, subobjparam.rpos_, subobjparam.q_, subobjparam.mesh_, subobjparam.material_, false, subobjparam.communicating_, subobjparam.infiniteMass_ );
auto cp = std::make_unique<mesh::pe::ConvexPolyhedron>( subobjparam.sid_, subobjparam.uid_, subobjparam.gpos_, subobjparam.q_, subobjparam.mesh_, subobjparam.material_, false, subobjparam.communicating_, subobjparam.infiniteMass_ );
cp->setLinearVel( subobjparam.v_ );
cp->setAngularVel( subobjparam.w_ );
cp->MPITrait.setOwner( subobjparam.mpiTrait_.owner_ );
......
......@@ -62,18 +62,18 @@ namespace pe {
* \param communicating specifies if the ConvexPolyhedron should take part in synchronization (syncNextNeighbour, syncShadowOwner)
* \param infiniteMass specifies if the ConvexPolyhedron has infinite mass and will be treated as an obstacle
*/
ConvexPolyhedron::ConvexPolyhedron( id_t sid, id_t uid, const Vec3& gpos, const Vec3& rpos, const Quat& q,
ConvexPolyhedron::ConvexPolyhedron( id_t sid, id_t uid, const Vec3& gpos, const Quat& q,
const TriangleMesh & mesh, MaterialID material,
const bool global, const bool communicating, const bool infiniteMass )
: GeomPrimitive( getStaticTypeID(), sid, uid, material ), // Initialization of the parent class
mesh_( mesh )
{
init( gpos, rpos, q, global, communicating, infiniteMass );
init( gpos, q, global, communicating, infiniteMass );
}
//*************************************************************************************************
void ConvexPolyhedron::init( const Vec3& gpos, const Vec3& rpos, const Quat& q,
void ConvexPolyhedron::init( const Vec3& gpos, const Quat& q,
const bool global, const bool communicating, const bool infiniteMass )
{
WALBERLA_ASSERT_FLOAT_EQUAL( (toWalberla( computeCentroid( mesh_ ) ) - Vec3() ).length(), real_t(0) );
......@@ -82,14 +82,21 @@ void ConvexPolyhedron::init( const Vec3& gpos, const Vec3& rpos, const Quat& q,
mesh_.request_face_normals();
mesh_.update_face_normals();
// Calculate the bounding sphere radius first, as setPosition will trigger a call to calcBoundingBox(), which needs it
real_t maxSqRadius(0);
for(auto vh : mesh_.vertices())
{
real_t sqRadius = mesh_.point( vh ).sqrnorm();
if( sqRadius > maxSqRadius )
maxSqRadius = sqRadius;
}
boundingSphereRadius_ = std::sqrt( maxSqRadius );
// Setting the center of the ConvexPolyhedron
gpos_ = gpos;
setPosition(gpos);
// Initializing the instantiated ConvexPolyhedron
rpos_ = rpos; // Setting the relative position
q_ = q; // Setting the orientation
R_ = q_.toRotationMatrix(); // Setting the rotation matrix
setOrientation(q);
setGlobal( global );
if (infiniteMass)
{
......@@ -103,14 +110,7 @@ void ConvexPolyhedron::init( const Vec3& gpos, const Vec3& rpos, const Quat& q,
setFinite( true );
// Setting the axis-aligned bounding box
real_t maxSqRadius(0);
for(auto vh : mesh_.vertices())
{
real_t sqRadius = mesh_.point( vh ).sqrnorm();
if( sqRadius > maxSqRadius )
maxSqRadius = sqRadius;
}
boundingSphereRadius_ = std::sqrt( maxSqRadius );
ConvexPolyhedron::calcBoundingBox();
octandVertices_[0] = supportVertex( TriangleMesh::Normal( real_t( 1), real_t( 1), real_t( 1) ), *mesh_.vertices_begin() );
......@@ -365,11 +365,12 @@ void ConvexPolyhedron::print( std::ostream& os, const char* tab ) const
// if( verboseMode )
// {
Mat3 R = getRotation();
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";
<< 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"
......
......@@ -65,7 +65,7 @@ public:
//**Constructors********************************************************************************
/*!\name Constructors */
//@{
ConvexPolyhedron( id_t sid, id_t uid, const Vec3& gpos, const Vec3& rpos, const Quat& q,
ConvexPolyhedron( id_t sid, id_t uid, const Vec3& gpos, const Quat& q,
const TriangleMesh & mesh, MaterialID material,
const bool global, const bool communicating, const bool infiniteMass );
//@}
......@@ -80,7 +80,7 @@ public:
//**********************************************************************************************
public:
void init( const Vec3& gpos, const Vec3& rpos, const Quat& q,
void init( const Vec3& gpos, const Quat& q,
const bool global, const bool communicating, const bool infiniteMass );
//**Get functions*******************************************************************************
......
......@@ -78,7 +78,7 @@ ConvexPolyhedronID createConvexPolyhedron( BodyStorage& globalStorage, BlockStor
WALBERLA_CHECK_EQUAL(communicating, false, "Global bodies can not be communicating!" );
WALBERLA_CHECK_EQUAL(infiniteMass, true, "Global bodies must have infinite mass!" );
auto cp = std::make_unique<ConvexPolyhedron>(sid, uid, gpos, Vec3(0,0,0), Quat(), mesh, material, global, false, true);
auto cp = std::make_unique<ConvexPolyhedron>(sid, uid, gpos, Quat(), mesh, material, global, false, true);
poly = static_cast<ConvexPolyhedronID>(&globalStorage.add(std::move(cp)));
} else
{
......@@ -88,7 +88,7 @@ ConvexPolyhedronID createConvexPolyhedron( BodyStorage& globalStorage, BlockStor
const id_t sid( UniqueID<RigidBody>::create() );
BodyStorage& bs = (*block.getData<Storage>(storageID))[0];
auto cp = std::make_unique<ConvexPolyhedron>(sid, uid, gpos, Vec3(0,0,0), Quat(), mesh, material, global, communicating, infiniteMass);
auto cp = std::make_unique<ConvexPolyhedron>(sid, uid, gpos, Quat(), mesh, material, global, communicating, infiniteMass);
cp->MPITrait.setOwner(Owner(MPIManager::instance()->rank(), block.getId().getID()));
poly = static_cast<ConvexPolyhedronID>(&bs.add( std::move(cp) ) );
}
......
......@@ -66,10 +66,6 @@ void marshal( mpi::SendBuffer& buffer, const RigidBody& obj ) {
buffer << obj.hasInfiniteMass();
buffer << obj.getPosition();
buffer << obj.hasSuperBody();
if( obj.hasSuperBody() )
{
buffer << obj.getRelPosition();
}
buffer << obj.getQuaternion();
if( !obj.hasSuperBody() )
{
......@@ -96,12 +92,6 @@ void unmarshal( mpi::RecvBuffer& buffer, RigidBodyParameters& objparam ) {
buffer >> objparam.infiniteMass_;
buffer >> objparam.gpos_;
buffer >> objparam.hasSuperBody_;
if( objparam.hasSuperBody_ )
{
buffer >> objparam.rpos_;
}
buffer >> objparam.q_;
if( !objparam.hasSuperBody_ )
......
......@@ -75,7 +75,7 @@ struct RigidBodyParameters {
MPIRigidBodyTraitParameter mpiTrait_;
bool communicating_, infiniteMass_;
id_t sid_, uid_;
Vec3 gpos_, rpos_, v_, w_;
Vec3 gpos_, v_, w_;
bool hasSuperBody_;
Quat q_;
};
......
......@@ -64,7 +64,7 @@ inline BoxPtr instantiate( mpi::RecvBuffer& buffer, const math::AABB& domain, co
BoxParameters subobjparam;
unmarshal( buffer, subobjparam );
correctBodyPosition(domain, block.center(), subobjparam.gpos_);
auto bx = std::make_unique<Box>( subobjparam.sid_, subobjparam.uid_, subobjparam.gpos_, subobjparam.rpos_, subobjparam.q_, subobjparam.lengths_, subobjparam.material_, false, subobjparam.communicating_, subobjparam.infiniteMass_ );
auto bx = std::make_unique<Box>( subobjparam.sid_, subobjparam.uid_, subobjparam.gpos_, subobjparam.q_, subobjparam.lengths_, subobjparam.material_, false, subobjparam.communicating_, subobjparam.infiniteMass_ );
bx->setLinearVel( subobjparam.v_ );
bx->setAngularVel( subobjparam.w_ );
bx->MPITrait.setOwner( subobjparam.mpiTrait_.owner_ );
......
......@@ -66,7 +66,7 @@ inline CapsulePtr instantiate( mpi::RecvBuffer& buffer, const math::AABB& domain
CapsuleParameters subobjparam;
unmarshal( buffer, subobjparam );
correctBodyPosition(domain, block.center(), subobjparam.gpos_);
auto cp = std::make_unique<Capsule>( subobjparam.sid_, subobjparam.uid_, subobjparam.gpos_, subobjparam.rpos_, subobjparam.q_, subobjparam.radius_, subobjparam.length_, subobjparam.material_, false, subobjparam.communicating_, subobjparam.infiniteMass_ );
auto cp = std::make_unique<Capsule>( subobjparam.sid_, subobjparam.uid_, subobjparam.gpos_, subobjparam.q_, subobjparam.radius_, subobjparam.length_, subobjparam.material_, false, subobjparam.communicating_, subobjparam.infiniteMass_ );
cp->setLinearVel( subobjparam.v_ );
cp->setAngularVel( subobjparam.w_ );
cp->MPITrait.setOwner( subobjparam.mpiTrait_.owner_ );
......
......@@ -66,7 +66,7 @@ inline EllipsoidPtr instantiate( mpi::RecvBuffer& buffer, const math::AABB& doma
EllipsoidParameters subobjparam;
unmarshal( buffer, subobjparam );
correctBodyPosition(domain, block.center(), subobjparam.gpos_);
auto el = std::make_unique<Ellipsoid>( subobjparam.sid_, subobjparam.uid_, subobjparam.gpos_, subobjparam.rpos_, subobjparam.q_, subobjparam.semiAxes_, subobjparam.material_, false, subobjparam.communicating_, subobjparam.infiniteMass_ );
auto el = std::make_unique<Ellipsoid>( subobjparam.sid_, subobjparam.uid_, subobjparam.gpos_, subobjparam.q_, subobjparam.semiAxes_, subobjparam.material_, false, subobjparam.communicating_, subobjparam.infiniteMass_ );
el->setLinearVel( subobjparam.v_ );
el->setAngularVel( subobjparam.w_ );
el->MPITrait.setOwner( subobjparam.mpiTrait_.owner_ );
......
......@@ -66,7 +66,7 @@ inline SpherePtr instantiate( mpi::RecvBuffer& buffer, const math::AABB& domain,
SphereParameters subobjparam;
unmarshal( buffer, subobjparam );
correctBodyPosition(domain, block.center(), subobjparam.gpos_);
auto sp = std::make_unique<Sphere>( subobjparam.sid_, subobjparam.uid_, subobjparam.gpos_, subobjparam.rpos_, subobjparam.q_, subobjparam.radius_, subobjparam.material_, false, subobjparam.communicating_, subobjparam.infiniteMass_ );
auto sp = std::make_unique<Sphere>( subobjparam.sid_, subobjparam.uid_, subobjparam.gpos_, subobjparam.q_, subobjparam.radius_, subobjparam.material_, false, subobjparam.communicating_, subobjparam.infiniteMass_ );
sp->setLinearVel( subobjparam.v_ );
sp->setAngularVel( subobjparam.w_ );
sp->MPITrait.setOwner( subobjparam.mpiTrait_.owner_ );
......
......@@ -66,7 +66,7 @@ inline SquirmerPtr instantiate( mpi::RecvBuffer& buffer, const math::AABB& domai
SquirmerParameters subobjparam;
unmarshal( buffer, subobjparam );
correctBodyPosition(domain, block.center(), subobjparam.gpos_);
auto sq = std::make_unique<Squirmer>( subobjparam.sid_, subobjparam.uid_, subobjparam.gpos_, subobjparam.rpos_, subobjparam.q_, subobjparam.radius_, subobjparam.squirmerVelocity_, subobjparam.squirmerBeta_, subobjparam.material_, false, subobjparam.communicating_, subobjparam.infiniteMass_ );
auto sq = std::make_unique<Squirmer>( subobjparam.sid_, subobjparam.uid_, subobjparam.gpos_, subobjparam.q_, subobjparam.radius_, subobjparam.squirmerVelocity_, subobjparam.squirmerBeta_, subobjparam.material_, false, subobjparam.communicating_, subobjparam.infiniteMass_ );
sq->setLinearVel( subobjparam.v_ );
sq->setAngularVel( subobjparam.w_ );
sq->MPITrait.setOwner( subobjparam.mpiTrait_.owner_ );
......
......@@ -119,8 +119,7 @@ inline std::unique_ptr<Union<BodyTypes...>> instantiate( mpi::RecvBuffer& buffer
correctBodyPosition(domain, block.center(), subobjparam.gpos_);
auto un = std::make_unique<Union<BodyTypes...>>( subobjparam.sid_,
subobjparam.uid_,
subobjparam.gpos_,
subobjparam.rpos_,
Vec3(),
subobjparam.q_,
false,
subobjparam.communicating_,
......@@ -140,6 +139,10 @@ inline std::unique_ptr<Union<BodyTypes...>> instantiate( mpi::RecvBuffer& buffer
un->setLinearVel( subobjparam.v_ );
un->setAngularVel( subobjparam.w_ );
newBody = un.get();
// Checks with global data of the union
WALBERLA_ASSERT_FLOAT_EQUAL(subobjparam.m_, un->getMass());
WALBERLA_ASSERT_FLOAT_EQUAL(subobjparam.I_, un->getBodyInertia());
WALBERLA_ASSERT_FLOAT_EQUAL(subobjparam.gpos_, un->getPosition());
return un;
}
......
......@@ -55,7 +55,7 @@ namespace pe {
* \param visible Specifies if the box is visible in a visualization.
* \param fixed \a true to fix the box, \a false to unfix it.
*/
Box::Box( id_t sid, id_t uid, const Vec3& gpos, const Vec3& rpos, const Quat& q,
Box::Box( id_t sid, id_t uid, const Vec3& gpos, const Quat& q,
const Vec3& lengths, MaterialID material,
const bool global, const bool communicating, const bool infiniteMass )
: GeomPrimitive( getStaticTypeID(), sid, uid, material ) // Initialization of the parent class
......@@ -70,10 +70,8 @@ Box::Box( id_t sid, id_t uid, const Vec3& gpos, const Vec3& rpos, const Quat& q,
WALBERLA_ASSERT_GREATER( lengths[2], real_t(0), "Invalid side length in z-dimension" );
// Initializing the instantiated box
gpos_ = gpos;
rpos_ = rpos; // Setting the relative position
q_ = q; // Setting the orientation
R_ = q_.toRotationMatrix(); // Setting the rotation matrix
setPosition(gpos);
setOrientation(q); // Setting the orientation
setGlobal( global );
if (infiniteMass)
......@@ -267,11 +265,12 @@ void Box::print( std::ostream& os, const char* tab ) const
//if( verboseMode )
{
Mat3 R = getRotation();
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";
<< 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"
......@@ -302,22 +301,23 @@ void Box::print( std::ostream& os, const char* tab ) const
void Box::calcBoundingBox()
{
using std::fabs;
const real_t xlength( real_t(0.5) * ( fabs(R_[0]*lengths_[0]) + fabs(R_[1]*lengths_[1]) + fabs(R_[2]*lengths_[2]) ) + contactThreshold );
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 );
Mat3 R = getRotation();
const real_t xlength( real_t(0.5) * ( fabs(R[0]*lengths_[0]) + fabs(R[1]*lengths_[1]) + fabs(R[2]*lengths_[2]) ) + contactThreshold );
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_[1] - ylength,
gpos_[2] - zlength,
gpos_[0] + xlength,
gpos_[1] + ylength,
gpos_[2] + zlength
getPosition()[0] - xlength,
getPosition()[1] - ylength,
getPosition()[2] - zlength,
getPosition()[0] + xlength,
getPosition()[1] + ylength,
getPosition()[2] + zlength
);
//WALBERLA_ASSERT( aabb_.isValid() , "Invalid bounding box detected" );
WALBERLA_ASSERT( aabb_.contains( gpos_ ), "Invalid bounding box detected" );
WALBERLA_ASSERT( aabb_.contains( getPosition() ), "Invalid bounding box detected" );
}
//*************************************************************************************************
......
......@@ -65,7 +65,7 @@ public:
//**Constructors********************************************************************************
/*!\name Constructors */
//@{
explicit Box( id_t sid, id_t uid, const Vec3& gpos, const Vec3& rpos, const Quat& q,
explicit Box( id_t sid, id_t uid, const Vec3& gpos, const Quat& q,
const Vec3& lengths, MaterialID material,
const bool global, const bool communicating, const bool infiniteMass );
//@}
......@@ -440,7 +440,7 @@ inline Vec3 Box::support( const Vec3& d ) const
math::sign(bfD[1])*lengths_[1]*real_t(0.5),
math::sign(bfD[2])*lengths_[2]*real_t(0.5) );
return gpos_ + vectorFromBFtoWF(relativSupport);
return getPosition() + vectorFromBFtoWF(relativSupport);
}
//*************************************************************************************************
......
......@@ -50,7 +50,7 @@ BoxID createBox( BodyStorage& globalStorage, BlockStorage& blocks, BlockDa
const id_t sid = UniqueID<RigidBody>::createGlobal();
WALBERLA_ASSERT_EQUAL(communicating, false);
WALBERLA_ASSERT_EQUAL(infiniteMass, true);
BoxPtr bx = std::make_unique<Box>(sid, uid, gpos, Vec3(0,0,0), Quat(), lengths, material, global, false, true);
BoxPtr bx = std::make_unique<Box>(sid, uid, gpos, Quat(), lengths, material, global, false, true);
box = static_cast<BoxID>(&globalStorage.add(std::move(bx)));
} else
{
......@@ -60,7 +60,7 @@ BoxID createBox( BodyStorage& globalStorage, BlockStorage& blocks, BlockDa
const id_t sid( UniqueID<RigidBody>::create() );
BodyStorage& bs = (*block.getData<Storage>(storageID))[0];
BoxPtr bx = std::make_unique<Box>(sid, uid, gpos, Vec3(0,0,0), Quat(), lengths, material, global, communicating, infiniteMass);
BoxPtr bx = std::make_unique<Box>(sid, uid, gpos, Quat(), lengths, material, global, communicating, infiniteMass);
bx->MPITrait.setOwner(Owner(MPIManager::instance()->rank(), block.getId().getID()));
box = static_cast<BoxID>(&bs.add(std::move(bx)));
}
......
......@@ -55,7 +55,7 @@ namespace pe {
*
* The capsule is created lying along the x-axis.
*/
Capsule::Capsule( id_t sid, id_t uid, const Vec3& gpos, const Vec3& rpos, const Quat& q,
Capsule::Capsule( id_t sid, id_t uid, const Vec3& gpos, const Quat& q,
real_t radius, real_t length, MaterialID material,
const bool global, const bool communicating, const bool infiniteMass )
: GeomPrimitive( getStaticTypeID(), sid, uid, material ) // Initializing the base object
......@@ -70,10 +70,8 @@ Capsule::Capsule( id_t sid, id_t uid, const Vec3& gpos, const Vec3& rpos, const
WALBERLA_ASSERT_GREATER( length, real_t(0), "Invalid capsule length" );
// Initializing the instantiated capsule
gpos_ = gpos;
rpos_ = rpos; // Setting the relative position
q_ = q; // Setting the orientation
R_ = q_.toRotationMatrix(); // Setting the rotation matrix
setPosition(gpos);
setOrientation(q);
setGlobal( global );
if (infiniteMass)
......@@ -171,18 +169,20 @@ bool Capsule::isSurfaceRelPointImpl( real_t px, real_t py, real_t pz ) const
*/
void Capsule::calcBoundingBox()
{
const real_t xlength( std::fabs( R_[0]*length_ )*real_t (0.5) + radius_ + contactThreshold );
const real_t ylength( std::fabs( R_[3]*length_ )*real_t (0.5) + radius_ + contactThreshold );
const real_t zlength( std::fabs( R_[6]*length_ )*real_t (0.5) + radius_ + contactThreshold );
Mat3 R = getRotation();
Vec3 gpos = getPosition();
const real_t xlength( std::fabs( R[0]*length_ )*real_t (0.5) + radius_ + contactThreshold );
const real_t ylength( std::fabs( R[3]*length_ )*real_t (0.5) + radius_ + contactThreshold );
const real_t zlength( std::fabs( R[6]*length_ )*real_t (0.5) + radius_ + contactThreshold );
aabb_ = math::AABB(
gpos_[0] - xlength,
gpos_[1] - ylength,
gpos_[2] - zlength,
gpos_[0] + xlength,
gpos_[1] + ylength,
gpos_[2] + zlength);
WALBERLA_ASSERT( aabb_.contains( gpos_ ), "Invalid bounding box detected" );
gpos[0] - xlength,
gpos[1] - ylength,
gpos[2] - zlength,
gpos[0] + xlength,
gpos[1] + ylength,
gpos[2] + zlength);
WALBERLA_ASSERT( aabb_.contains( getPosition() ), "Invalid bounding box detected" );
}
//*************************************************************************************************
......@@ -243,11 +243,12 @@ void Capsule::print( std::ostream& os, const char* tab ) const
<< tab << " Linear velocity = " << getLinearVel() << "\n"
<< tab << " Angular velocity = " << getAngularVel() << "\n";
Mat3 R = getRotation();
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";
<< 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"
......
......@@ -70,7 +70,7 @@ public:
//**Constructors********************************************************************************
/*!\name Constructors */
//@{
explicit Capsule( id_t sid, id_t uid, const Vec3& gpos, const Vec3& rpos, const Quat& q,
explicit Capsule( id_t sid, id_t uid, const Vec3& gpos, const Quat& q,
real_t radius, real_t length, MaterialID material,
const bool global, const bool communicating, const bool infiniteMass );
//@}
......@@ -97,7 +97,6 @@ public:
/*!\name Utility functions */
//@{
inline virtual Vec3 support( const Vec3& d ) const;
inline virtual Vec3 supportContactThreshold( const Vec3& d ) const;
//@}
//**********************************************************************************************
......@@ -278,30 +277,11 @@ inline Vec3 Capsule::support( const Vec3& d ) const
const Vec3 supportSegment = Vec3( math::sign(bfD[0])*length_*real_t (0.5), real_t (0.0), real_t (0.0));
const Vec3 supportSphere = radius_ * dnorm;
return gpos_ + vectorFromBFtoWF(supportSegment) + supportSphere;
return getPosition() + vectorFromBFtoWF(supportSegment) + supportSphere;
}
//*************************************************************************************************
//*************************************************************************************************
/*!\brief Estimates the point which is farthest in direction \a d.
*
* \param d The normalized search direction in world-frame coordinates
* \return The support point in world-frame coordinates in direction a\ d extended by a vector in
* direction \a d of length \a pe::contactThreshold.
*/
inline Vec3 Capsule::supportContactThreshold( const Vec3& d ) const
{
auto len = d.sqrLength();
if (math::equal(len, real_t(0)))
return Vec3(0,0,0);
return support(d) + d*contactThreshold;
}
//*************************************************************************************************
//=================================================================================================
//
......
......@@ -50,7 +50,7 @@ CapsuleID createCapsule( BodyStorage& globalStorage, BlockStorage& blocks, Blo
const id_t sid = UniqueID<RigidBody>::createGlobal();
WALBERLA_ASSERT_EQUAL(communicating, false);
WALBERLA_ASSERT_EQUAL(infiniteMass, true);
CapsulePtr cp = std::make_unique<Capsule>(sid, uid, gpos, Vec3(0,0,0), Quat(), radius, length, material, global, false, true);
CapsulePtr cp = std::make_unique<Capsule>(sid, uid, gpos, Quat(), radius, length, material, global, false, true);
capsule = static_cast<CapsuleID>(&globalStorage.add(std::move(cp)));
} else
{
......@@ -60,7 +60,7 @@ CapsuleID createCapsule( BodyStorage& globalStorage, BlockStorage& blocks, Blo
const id_t sid( UniqueID<RigidBody>::create() );
BodyStorage& bs = (*block.getData<Storage>(storageID))[0];
CapsulePtr cp = std::make_unique<Capsule>(sid, uid, gpos, Vec3(0,0,0), Quat(), radius, length, material, global, communicating, infiniteMass);
CapsulePtr cp = std::make_unique<Capsule>(sid, uid, gpos, Quat(), radius, length, material, global, communicating, infiniteMass);
cp->MPITrait.setOwner(Owner(MPIManager::instance()->rank(), block.getId().getID()));
capsule = static_cast<CapsuleID>(&bs.add( std::move(cp) ));
}
......
......@@ -71,9 +71,8 @@ CylindricalBoundary::CylindricalBoundary( id_t sid, id_t uid, const Vec3& gpos,
WALBERLA_ASSERT_GREATER( radius, real_t(0), "Invalid cylinder radius" );
// Initializing the instantiated cylinder
gpos_ = gpos;
q_ = Quat(); // Setting the orientation
R_ = q_.toRotationMatrix(); // Setting the rotation matrix
setPosition(gpos);
setOrientation(Quat());
// Setting the axis-aligned bounding box
CylindricalBoundary::calcBoundingBox();
......@@ -147,7 +146,7 @@ void CylindricalBoundary::calcBoundingBox()
+math::Limits<real_t>::inf(),
+math::Limits<real_t>::inf());
WALBERLA_ASSERT( aabb_.contains( gpos_ ), "Invalid bounding box detected" );
WALBERLA_ASSERT( aabb_.contains( getPosition() ), "Invalid bounding box detected" );
}
//*************************************************************************************************
......@@ -179,12 +178,12 @@ void CylindricalBoundary::print( std::ostream& os, const char* tab ) const
<< tab << " Global position = " << getPosition() << "\n"
<< tab << " Linear velocity = " << getLinearVel() << "\n"
<< tab << " Angular velocity = " << getAngularVel() << "\n";
Mat3 R = getRotation();
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";
<< 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";
}
//*************************************************************************************************
......
......@@ -58,12 +58,12 @@ namespace pe {
* \param communicating specifies if the Ellipsoid should take part in synchronization (syncNextNeighbour, syncShadowOwner)
* \param infiniteMass specifies if the Ellipsoid has infinite mass and will be treated as an obstacle
*/
Ellipsoid::Ellipsoid( id_t sid, id_t uid, const Vec3& gpos, const Vec3& rpos, const Quat& q,
Ellipsoid::Ellipsoid( id_t sid, id_t uid, const Vec3& gpos, const Quat& q,
const Vec3& semiAxes, MaterialID material,
const bool global, const bool communicating, const bool infiniteMass )
: Ellipsoid::Ellipsoid( getStaticTypeID(), sid, uid, gpos, rpos, q, semiAxes, material, global, communicating, infiniteMass )
: Ellipsoid::Ellipsoid( getStaticTypeID(), sid, uid, gpos, q, semiAxes, material, global, communicating, infiniteMass )
{}
Ellipsoid::Ellipsoid( id_t const typeId, id_t sid, id_t uid, const Vec3& gpos, const Vec3& rpos, const Quat& q,
Ellipsoid::Ellipsoid( id_t const typeId, id_t sid, id_t uid, const Vec3& gpos, const Quat& q,
const Vec3& semiAxes, MaterialID material,
const bool global, const bool communicating, const bool infiniteMass )
: GeomPrimitive( typeId, sid, uid, material ) // Initialization of the parent class
......@@ -76,13 +76,10 @@ Ellipsoid::Ellipsoid( id_t const typeId, id_t sid, id_t uid, const Vec3& gpos, c
WALBERLA_ASSERT( semiAxes_[0] > real_c(0), "Invalid Ellipsoid radius" );
WALBERLA_ASSERT( semiAxes_[1] > real_c(0), "Invalid Ellipsoid radius" );
WALBERLA_ASSERT( semiAxes_[2] > real_c(0), "Invalid Ellipsoid radius" );
// Setting the center of the Ellipsoid
gpos_ = gpos;
// Initializing the instantiated Ellipsoid
rpos_ = rpos; // Setting the relative position
q_ = q; // Setting the orientation
R_ = q_.toRotationMatrix(); // Setting the rotation matrix
// Setting the center of the Ellipsoid
setPosition(gpos);
setOrientation(q);
setGlobal( global );
if (infiniteMass)
......@@ -182,6 +179,8 @@ void Ellipsoid::print( std::ostream& os, const char* tab ) const
{
using std::setw;
Mat3 R = getRotation();
os << tab << " Ellipsoid " << uid_ << " with semi-axis " << semiAxes_ << "\n";
os << tab << " Fixed: " << isFixed() << " , sleeping: " << !isAwake() << "\n";
......@@ -199,9 +198,9 @@ void Ellipsoid::print( std::ostream& os, const char* tab ) const
// {
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";
<< tab << " Rotation matrix = ( " << setw(9) << R[1] << " , " << 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"
......
......@@ -65,10 +65,10 @@ public:
//**Constructors********************************************************************************
/*!\name Constructors */
//@{
explicit Ellipsoid( id_t sid, id_t uid, const Vec3& gpos, const Vec3& rpos, const Quat& q,
explicit Ellipsoid( id_t sid, id_t uid, const Vec3& gpos, const Quat& q,
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,
explicit Ellipsoid( id_t const typeID, id_t sid, id_t uid, const Vec3& rpos, const Quat& q,
const Vec3& semiAxes, MaterialID material,
const bool global, const bool communicating, const bool infiniteMass );
//@}
......@@ -259,20 +259,21 @@ inline real_t Ellipsoid::calcDensity(const Vec3& semiAxes, real_t mass )
inline void Ellipsoid::calcBoundingBox()
{
using std::fabs;
const real_t xlength( fabs(R_[0]*semiAxes_[0]) + fabs(R_[1]*semiAxes_[1]) + fabs(R_[2]*semiAxes_[2]) + contactThreshold );
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 );
Mat3 R(getRotation());
Vec3 gpos = getPosition();
const real_t xlength( fabs(R[0]*semiAxes_[0]) + fabs(R[1]*semiAxes_[1]) + fabs(R[2]*semiAxes_[2]) + contactThreshold );
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_[1] - ylength,
gpos_[2] - zlength,
gpos_[0] + xlength,
gpos_[1] + ylength,
gpos_[2] + zlength
gpos[0] - xlength,
gpos[1] - ylength,
gpos[2] - zlength,
gpos[0] + xlength,
gpos[1] + ylength,
gpos[2] + zlength
);