From f0a69f7da7784ffd1e7c74136073c181de28f8e0 Mon Sep 17 00:00:00 2001
From: Tobias Leemann <tobias.leemann@fau.de>
Date: Wed, 10 Jul 2019 13:35:21 +0200
Subject: [PATCH] [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.
---
src/mesh/pe/communication/ConvexPolyhedron.h | 2 +-
src/mesh/pe/rigid_body/ConvexPolyhedron.cpp | 39 +-
src/mesh/pe/rigid_body/ConvexPolyhedron.h | 4 +-
.../pe/rigid_body/ConvexPolyhedronFactory.cpp | 4 +-
src/pe/communication/Marshalling.cpp | 10 -
src/pe/communication/Marshalling.h | 2 +-
src/pe/communication/rigidbody/Box.h | 2 +-
src/pe/communication/rigidbody/Capsule.h | 2 +-
src/pe/communication/rigidbody/Ellipsoid.h | 2 +-
src/pe/communication/rigidbody/Sphere.h | 2 +-
src/pe/communication/rigidbody/Squirmer.h | 2 +-
src/pe/communication/rigidbody/Union.h | 7 +-
src/pe/rigidbody/Box.cpp | 38 +-
src/pe/rigidbody/Box.h | 4 +-
src/pe/rigidbody/BoxFactory.cpp | 4 +-
src/pe/rigidbody/Capsule.cpp | 39 +-
src/pe/rigidbody/Capsule.h | 24 +-
src/pe/rigidbody/CapsuleFactory.cpp | 4 +-
src/pe/rigidbody/CylindricalBoundary.cpp | 15 +-
src/pe/rigidbody/Ellipsoid.cpp | 23 +-
src/pe/rigidbody/Ellipsoid.h | 27 +-
src/pe/rigidbody/EllipsoidFactory.cpp | 4 +-
src/pe/rigidbody/Plane.cpp | 169 ++------
src/pe/rigidbody/Plane.h | 8 -
src/pe/rigidbody/RigidBody.cpp | 9 +-
src/pe/rigidbody/RigidBody.h | 376 ++++++++----------
src/pe/rigidbody/Sphere.cpp | 23 +-
src/pe/rigidbody/Sphere.h | 48 +--
src/pe/rigidbody/SphereFactory.cpp | 4 +-
src/pe/rigidbody/Squirmer.cpp | 8 +-
src/pe/rigidbody/Squirmer.h | 2 +-
src/pe/rigidbody/SquirmerFactory.cpp | 4 +-
src/pe/rigidbody/Union.h | 309 ++++----------
src/pe/rigidbody/UnionFactory.h | 10 +-
tests/mesh/MeshMarshalling.cpp | 42 +-
tests/mesh/MeshPeRaytracing.cpp | 4 +-
tests/pe/BodyFlags.cpp | 4 +-
tests/pe/BodyStorage.cpp | 4 +-
tests/pe/CMakeLists.txt | 5 +
tests/pe/CheckVitalParameters.cpp | 2 +-
tests/pe/Collision.cpp | 32 +-
tests/pe/CollisionTobiasGJK.cpp | 50 +--
tests/pe/Marshalling.cpp | 14 +-
tests/pe/Raytracing.cpp | 20 +-
tests/pe/RigidBody.cpp | 12 +-
tests/pe/SetBodyTypeIDs.cpp | 4 +-
tests/pe/SimpleCCD.cpp | 6 +-
tests/pe/Synchronization.cpp | 2 +-
tests/pe/SynchronizationLargeBody.cpp | 2 +-
tests/pe/Union.cpp | 185 ++++++++-
tests/pe/UnionBehavior.cpp | 259 ++++++++++++
tests/pe/VolumeInertia.cpp | 8 +-
.../geometry/PeIntersectionRatioTest.cpp | 9 +-
53 files changed, 1004 insertions(+), 890 deletions(-)
create mode 100644 tests/pe/UnionBehavior.cpp
diff --git a/src/mesh/pe/communication/ConvexPolyhedron.h b/src/mesh/pe/communication/ConvexPolyhedron.h
index df14bb176..1527530b4 100644
--- a/src/mesh/pe/communication/ConvexPolyhedron.h
+++ b/src/mesh/pe/communication/ConvexPolyhedron.h
@@ -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_ );
diff --git a/src/mesh/pe/rigid_body/ConvexPolyhedron.cpp b/src/mesh/pe/rigid_body/ConvexPolyhedron.cpp
index eb3c3c3ee..62ce1d1c5 100644
--- a/src/mesh/pe/rigid_body/ConvexPolyhedron.cpp
+++ b/src/mesh/pe/rigid_body/ConvexPolyhedron.cpp
@@ -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"
diff --git a/src/mesh/pe/rigid_body/ConvexPolyhedron.h b/src/mesh/pe/rigid_body/ConvexPolyhedron.h
index 48f7c8ff1..10e6f62a7 100644
--- a/src/mesh/pe/rigid_body/ConvexPolyhedron.h
+++ b/src/mesh/pe/rigid_body/ConvexPolyhedron.h
@@ -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*******************************************************************************
diff --git a/src/mesh/pe/rigid_body/ConvexPolyhedronFactory.cpp b/src/mesh/pe/rigid_body/ConvexPolyhedronFactory.cpp
index 6f29a33cf..4229a58f1 100644
--- a/src/mesh/pe/rigid_body/ConvexPolyhedronFactory.cpp
+++ b/src/mesh/pe/rigid_body/ConvexPolyhedronFactory.cpp
@@ -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) ) );
}
diff --git a/src/pe/communication/Marshalling.cpp b/src/pe/communication/Marshalling.cpp
index bdeb9fd9d..9a8e930f6 100644
--- a/src/pe/communication/Marshalling.cpp
+++ b/src/pe/communication/Marshalling.cpp
@@ -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_ )
diff --git a/src/pe/communication/Marshalling.h b/src/pe/communication/Marshalling.h
index 225ffaa4f..413636274 100644
--- a/src/pe/communication/Marshalling.h
+++ b/src/pe/communication/Marshalling.h
@@ -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_;
};
diff --git a/src/pe/communication/rigidbody/Box.h b/src/pe/communication/rigidbody/Box.h
index 77f543cbc..78d698f48 100644
--- a/src/pe/communication/rigidbody/Box.h
+++ b/src/pe/communication/rigidbody/Box.h
@@ -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_ );
diff --git a/src/pe/communication/rigidbody/Capsule.h b/src/pe/communication/rigidbody/Capsule.h
index 3ef02f07d..e319d579a 100644
--- a/src/pe/communication/rigidbody/Capsule.h
+++ b/src/pe/communication/rigidbody/Capsule.h
@@ -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_ );
diff --git a/src/pe/communication/rigidbody/Ellipsoid.h b/src/pe/communication/rigidbody/Ellipsoid.h
index a292c5b1d..de8787b45 100644
--- a/src/pe/communication/rigidbody/Ellipsoid.h
+++ b/src/pe/communication/rigidbody/Ellipsoid.h
@@ -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_ );
diff --git a/src/pe/communication/rigidbody/Sphere.h b/src/pe/communication/rigidbody/Sphere.h
index 3e7da5601..846a621e5 100644
--- a/src/pe/communication/rigidbody/Sphere.h
+++ b/src/pe/communication/rigidbody/Sphere.h
@@ -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_ );
diff --git a/src/pe/communication/rigidbody/Squirmer.h b/src/pe/communication/rigidbody/Squirmer.h
index 3c9b3fbe0..171b9c742 100644
--- a/src/pe/communication/rigidbody/Squirmer.h
+++ b/src/pe/communication/rigidbody/Squirmer.h
@@ -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_ );
diff --git a/src/pe/communication/rigidbody/Union.h b/src/pe/communication/rigidbody/Union.h
index bfaf37881..644f16a8f 100644
--- a/src/pe/communication/rigidbody/Union.h
+++ b/src/pe/communication/rigidbody/Union.h
@@ -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;
}
diff --git a/src/pe/rigidbody/Box.cpp b/src/pe/rigidbody/Box.cpp
index ffd27fc8e..0705a4ca1 100644
--- a/src/pe/rigidbody/Box.cpp
+++ b/src/pe/rigidbody/Box.cpp
@@ -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" );
}
+
//*************************************************************************************************
diff --git a/src/pe/rigidbody/Box.h b/src/pe/rigidbody/Box.h
index 5da39e784..266aead5d 100644
--- a/src/pe/rigidbody/Box.h
+++ b/src/pe/rigidbody/Box.h
@@ -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);
}
//*************************************************************************************************
diff --git a/src/pe/rigidbody/BoxFactory.cpp b/src/pe/rigidbody/BoxFactory.cpp
index 6f471ba55..2a7036004 100644
--- a/src/pe/rigidbody/BoxFactory.cpp
+++ b/src/pe/rigidbody/BoxFactory.cpp
@@ -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)));
}
diff --git a/src/pe/rigidbody/Capsule.cpp b/src/pe/rigidbody/Capsule.cpp
index 7b58b0ce6..f1d9ad9ee 100644
--- a/src/pe/rigidbody/Capsule.cpp
+++ b/src/pe/rigidbody/Capsule.cpp
@@ -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"
diff --git a/src/pe/rigidbody/Capsule.h b/src/pe/rigidbody/Capsule.h
index c2d4c6b46..8718f4c90 100644
--- a/src/pe/rigidbody/Capsule.h
+++ b/src/pe/rigidbody/Capsule.h
@@ -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;
-}
-//*************************************************************************************************
-
-
-
//=================================================================================================
//
diff --git a/src/pe/rigidbody/CapsuleFactory.cpp b/src/pe/rigidbody/CapsuleFactory.cpp
index 325af9701..349e30b2f 100644
--- a/src/pe/rigidbody/CapsuleFactory.cpp
+++ b/src/pe/rigidbody/CapsuleFactory.cpp
@@ -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) ));
}
diff --git a/src/pe/rigidbody/CylindricalBoundary.cpp b/src/pe/rigidbody/CylindricalBoundary.cpp
index 66f639e49..58f4f97ea 100644
--- a/src/pe/rigidbody/CylindricalBoundary.cpp
+++ b/src/pe/rigidbody/CylindricalBoundary.cpp
@@ -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";
}
//*************************************************************************************************
diff --git a/src/pe/rigidbody/Ellipsoid.cpp b/src/pe/rigidbody/Ellipsoid.cpp
index 3c68d5875..d47b400d5 100644
--- a/src/pe/rigidbody/Ellipsoid.cpp
+++ b/src/pe/rigidbody/Ellipsoid.cpp
@@ -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"
diff --git a/src/pe/rigidbody/Ellipsoid.h b/src/pe/rigidbody/Ellipsoid.h
index 573d77348..1cda47381 100644
--- a/src/pe/rigidbody/Ellipsoid.h
+++ b/src/pe/rigidbody/Ellipsoid.h
@@ -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
);
// 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_ );
+ WALBERLA_ASSERT( aabb_.contains( gpos ), "Invalid bounding box detected("<< getSystemID() <<")\n" << "pos: " << getPosition() << "\nlengths: " << xlength << "," << ylength << ", " << zlength<< "\nvel: " << getLinearVel() << "\nbox: " << aabb_ );
}
//*************************************************************************************************
diff --git a/src/pe/rigidbody/EllipsoidFactory.cpp b/src/pe/rigidbody/EllipsoidFactory.cpp
index a5e46b350..2136a8414 100644
--- a/src/pe/rigidbody/EllipsoidFactory.cpp
+++ b/src/pe/rigidbody/EllipsoidFactory.cpp
@@ -48,7 +48,7 @@ EllipsoidID createEllipsoid( BodyStorage& globalStorage, BlockStorage& blocks, B
const id_t sid = UniqueID<RigidBody>::createGlobal();
WALBERLA_ASSERT_EQUAL(communicating, false);
WALBERLA_ASSERT_EQUAL(infiniteMass, true);
- EllipsoidPtr el = std::make_unique<Ellipsoid>(sid, uid, gpos, Vec3(0,0,0), Quat(), semiAxes, material, global, false, true);
+ EllipsoidPtr el = std::make_unique<Ellipsoid>(sid, uid, gpos, Quat(), semiAxes, material, global, false, true);
ellipsoid = static_cast<EllipsoidID>(&globalStorage.add(std::move(el)));
} else
{
@@ -58,7 +58,7 @@ EllipsoidID createEllipsoid( BodyStorage& globalStorage, BlockStorage& blocks, B
const id_t sid( UniqueID<RigidBody>::create() );
BodyStorage& bs = (*block.getData<Storage>(storageID))[0];
- EllipsoidPtr el = std::make_unique<Ellipsoid>(sid, uid, gpos, Vec3(0,0,0), Quat(), semiAxes, material, global, communicating, infiniteMass);
+ EllipsoidPtr el = std::make_unique<Ellipsoid>(sid, uid, gpos, Quat(), semiAxes, material, global, communicating, infiniteMass);
el->MPITrait.setOwner(Owner(MPIManager::instance()->rank(), block.getId().getID()));
ellipsoid = static_cast<EllipsoidID>(&bs.add(std::move(el)));
}
diff --git a/src/pe/rigidbody/Plane.cpp b/src/pe/rigidbody/Plane.cpp
index 872edbe3e..0d501c60d 100644
--- a/src/pe/rigidbody/Plane.cpp
+++ b/src/pe/rigidbody/Plane.cpp
@@ -87,17 +87,15 @@ Plane::Plane( id_t sid, id_t uid,
WALBERLA_ASSERT_FLOAT_EQUAL( normal.sqrLength(), real_c(1) , "Invalid plane normal" );
// Setting the global position (anchor point) of the plane
- gpos_ = gpos;
+ setPosition(gpos);
// Calculating the orientation and rotation
// The default normal of a plane is <0,0,1>. The rotation of the plane is calculated
// as the rotation of this default normal to the specified normal.
if( normal[0]*normal[0] + normal[1]*normal[1] < math::Limits<real_t>::accuracy() )
- q_ = Quat( Vec3( 1, 0, 0 ), std::acos( normal[2] ) );
+ setOrientation( Quat(Vec3( 1, 0, 0 ), std::acos( normal[2] ) ) );
else
- q_ = Quat( Vec3( -normal[1], normal[0], real_c(0) ), std::acos( normal[2] ) );
-
- R_ = q_.toRotationMatrix();
+ setOrientation( Quat( Vec3( -normal[1], normal[0], real_c(0) ), std::acos( normal[2] ) ) );
// Setting the axis-aligned bounding box
Plane::calcBoundingBox();
@@ -175,11 +173,8 @@ void Plane::calcBoundingBox()
*/
void Plane::setPositionImpl( real_t px, real_t py, real_t pz )
{
- gpos_[0] = px;
- gpos_[1] = py;
- gpos_[2] = pz;
- d_ = normal_ * gpos_;
-
+ RigidBody::setPositionImpl(px,py,pz);
+ d_ = normal_ * getPosition();
Plane::calcBoundingBox(); // Updating the axis-aligned bounding box of the plane
#if MOBILE_INFINITE
wake(); // Waking the box from sleep mode
@@ -210,16 +205,15 @@ void Plane::setPositionImpl( real_t px, real_t py, real_t pz )
*/
void Plane::setOrientationImpl( real_t r, real_t i, real_t j, real_t k )
{
- q_.set(r, i, j, k); // Updating the orientation of the plane
- R_ = q_.toRotationMatrix(); // Updating the rotation of the plane
-
+ RigidBody::setOrientationImpl(r,i,j,k);
+ Mat3 R = getRotation();
// Updating the normal of the plane ( R * <0,0,1> )
- normal_[0] = R_[2];
- normal_[1] = R_[5];
- normal_[2] = R_[8];
+ normal_[0] = R[2];
+ normal_[1] = R[5];
+ normal_[2] = R[8];
// Updating the displacement from the origin
- d_ = normal_ * gpos_;
+ d_ = normal_ * getPosition();
Plane::calcBoundingBox(); // Updating the axis-aligned bounding box of the plane
#if MOBILE_INFINITE
@@ -256,10 +250,9 @@ void Plane::setOrientationImpl( real_t r, real_t i, real_t j, real_t k )
*/
void Plane::translateImpl( real_t dx, real_t dy, real_t dz )
{
- gpos_[0] += dx;
- gpos_[1] += dy;
- gpos_[2] += dz;
- d_ = normal_ * gpos_;
+
+ setPosition(getPosition() + Vec3(dx,dy,dz));
+ d_ = normal_ * getPosition();
Plane::calcBoundingBox(); // Updating the axis-aligned bounding box
#if MOBILE_INFINITE
@@ -301,16 +294,16 @@ void Plane::rotateImpl( const Quat& dq )
// if( ExclusiveSection::isActive() )
// throw std::logic_error( "Invalid rotation of a plane inside an exclusive section" );
- q_ = dq * q_; // Updating the orientation of the plane
- R_ = q_.toRotationMatrix(); // Updating the rotation of the plane
+ setOrientation(dq * getQuaternion());
// Updating the normal of the plane ( R * <0,0,1> )
- normal_[0] = R_[2];
- normal_[1] = R_[5];
- normal_[2] = R_[8];
+ Mat3 R = getRotation();
+ normal_[0] = R[2];
+ normal_[1] = R[5];
+ normal_[2] = R[8];
// Updating the displacement from the origin
- d_ = normal_ * gpos_;
+ d_ = normal_ * getPosition();
Plane::calcBoundingBox(); // Updating the axis-aligned bounding box of the plane
#if MOBILE_INFINITE
@@ -342,14 +335,16 @@ void Plane::rotateImpl( const Quat& dq )
*/
void Plane::rotateAroundOriginImpl( const Quat& dq )
{
- gpos_ = dq.rotate( gpos_ ); // Updating the global position of the plane
- q_ = dq * q_; // Updating the orientation of the plane
- R_ = q_.toRotationMatrix(); // Updating the rotation of the plane
+
+ setPosition(dq.rotate( getPosition() )); // Updating the global position of the plane
+ setOrientation(dq * getQuaternion()); // Updating the orientation of the plane
+
// Updating the normal of the plane ( R * <0,0,1> )
- normal_[0] = R_[2];
- normal_[1] = R_[5];
- normal_[2] = R_[8];
+ Mat3 R = getRotation();
+ normal_[0] = R[2];
+ normal_[1] = R[5];
+ normal_[2] = R[8];
Plane::calcBoundingBox(); // Updating the axis-aligned bounding box of the plane
signalRotation(); // Signaling the change of orientation to the superordinate body
@@ -382,19 +377,19 @@ void Plane::rotateAroundOriginImpl( const Quat& dq )
*/
void Plane::rotateAroundPointImpl( const Vec3& point, const Quat& dq )
{
- const Vec3 dp( gpos_ - point );
+ const Vec3 dp( getPosition() - point );
- gpos_ = point + dq.rotate( dp ); // Updating the global position of the box
- q_ = dq * q_; // Updating the orientation of the box
- R_ = q_.toRotationMatrix(); // Updating the rotation of the box
+ setPosition(point + dq.rotate( dp ) ); // Updating the global position of the box
+ setOrientation(dq * getQuaternion()); // Updating the orientation of the box
// Updating the normal of the plane ( R * <0,0,1> )
- normal_[0] = R_[2];
- normal_[1] = R_[5];
- normal_[2] = R_[8];
+ Mat3 R = getRotation();
+ normal_[0] = R[2];
+ normal_[1] = R[5];
+ normal_[2] = R[8];
// Updating the displacement from the origin
- d_ = normal_ * gpos_;
+ d_ = normal_ * getPosition();
Plane::calcBoundingBox(); // Updating the axis-aligned bounding box of the plane
#if MOBILE_INFINITE
@@ -447,85 +442,6 @@ bool Plane::isSurfaceRelPointImpl( real_t /*px*/, real_t /*py*/, real_t pz ) con
//*************************************************************************************************
-
-
-
-//=================================================================================================
-//
-// SIMULATION FUNCTIONS
-//
-//=================================================================================================
-
-//*************************************************************************************************
-/*!\brief Translation update of a subordinate plane.
- *
- * \param dp Change in the global position of the superordinate body.
- * \return void
- *
- * This update function is triggered by the superordinate body in case of a translational
- * movement. This movement involves a change of the global position, the displacement from
- * the origin and the axis-aligned bounding box.
- */
-void Plane::update( const Vec3& dp )
-{
- // Checking the state of the plane
- WALBERLA_ASSERT( checkInvariants(), "Invalid plane state detected" );
- WALBERLA_ASSERT( hasSuperBody(), "Invalid superordinate body detected" );
-
- // Updating the global position and the displacement
- gpos_ += dp;
- d_ = normal_ * gpos_;
-
- // Setting the axis-aligned bounding box
- Plane::calcBoundingBox();
-
- // Checking the state of the plane
- WALBERLA_ASSERT( checkInvariants(), "Invalid plane state detected" );
-}
-//*************************************************************************************************
-
-
-//*************************************************************************************************
-/*!\brief Rotation update of a subordinate plane.
- *
- * \param dq Change in the orientation of the superordinate body.
- * \return void
- *
- * This update function is triggered by the superordinate body in case of a rotational
- * movement. This movement involves a change in the global position, the displacement from
- * the origin, the orientation/rotation, the normal and the axis-aligned bounding box of
- * the plane.
- */
-void Plane::update( const Quat& dq )
-{
- // Checking the state of the plane
- WALBERLA_ASSERT( checkInvariants(), "Invalid plane state detected" );
- WALBERLA_ASSERT( hasSuperBody(), "Invalid superordinate body detected" );
-
- // Calculating the new orientation and rotation
- q_ = dq * q_;
- R_ = q_.toRotationMatrix();
-
- // Updating the normal of the plane ( R * <0,0,1> )
- normal_[0] = R_[2];
- normal_[1] = R_[5];
- normal_[2] = R_[8];
-
- // Updating the global position and the displacement
- gpos_ = sb_->getPosition() + ( sb_->getRotation() * rpos_ );
- d_ = normal_ * gpos_;
-
- // Setting the axis-aligned bounding box
- Plane::calcBoundingBox();
-
- // Checking the state of the plane
- WALBERLA_ASSERT( checkInvariants(), "Invalid plane state detected" );
-}
-//*************************************************************************************************
-
-
-
-
//=================================================================================================
//
// OUTPUT FUNCTIONS
@@ -546,14 +462,15 @@ void Plane::print( std::ostream& os, const char* tab ) const
os << tab << " Plane " << uid_ << " with normal " << normal_ << " and displacement " << d_ << "\n"
<< tab << " System ID = " << sid_ << "\n"
<< tab << " Material = " << Material::getName( material_ ) << "\n"
- << tab << " Global position = " << gpos_ << "\n";
+ << tab << " Global position = " << getPosition() << "\n";
- os << tab << " Relative position = " << rpos_ << "\n"
+ os << tab << " Relative position = " << getRelPosition() << "\n"
<< tab << " Bounding box = " << aabb_ << "\n"
- << tab << " Quaternion = " << q_ << "\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 << " Quaternion = " << getQuaternion()<< "\n";
+ Mat3 R = getRotation();
+ os << 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";
}
//*************************************************************************************************
diff --git a/src/pe/rigidbody/Plane.h b/src/pe/rigidbody/Plane.h
index 589d11038..94368f79f 100644
--- a/src/pe/rigidbody/Plane.h
+++ b/src/pe/rigidbody/Plane.h
@@ -138,14 +138,6 @@ protected:
//@}
//**********************************************************************************************
- //**Simulation functions************************************************************************
- /*!\name Simulation functions */
- //@{
- virtual void update( const Vec3& dp ); // Translation update of a subordinate plane
- virtual void update( const Quat& dq ); // Rotation update of a subordinate plane
- //@}
- //**********************************************************************************************
-
//**Member variables****************************************************************************
/*!\name Member variables */
//@{
diff --git a/src/pe/rigidbody/RigidBody.cpp b/src/pe/rigidbody/RigidBody.cpp
index f30f35ce6..78743f796 100644
--- a/src/pe/rigidbody/RigidBody.cpp
+++ b/src/pe/rigidbody/RigidBody.cpp
@@ -37,17 +37,13 @@ RigidBody::RigidBody( id_t const typeID, id_t sid, id_t uid )
, mass_( 0 ) // Total mass of the rigid body
, invMass_( 0 ) // Inverse total mass of the rigid body
, motion_(sleepThreshold) // The current motion of the rigid body
- , gpos_() // Global position of the center of mass
- , rpos_() // Relative position within the body frame of the superordinate body
, v_() // Linear velocity
, w_() // Angular velocity
, force_() // Total force
, torque_() // Total torque
, I_( real_c(0) ) // Moment of inertia
, Iinv_( real_c(0) ) // Inverse moment of inertia
- , q_() // Orientation of the body frame
- , R_() // Rigid body rotation
- , manager_(nullptr) // The rigid body manager responsible for the rigid body
+ , manager_(nullptr) // The rigid body manager responsible for the rigid body
, finite_ (true) // Finiteness flag
, visible_(true) // Visibility flag
, remote_ (false) // Remote flag
@@ -56,6 +52,9 @@ RigidBody::RigidBody( id_t const typeID, id_t sid, id_t uid )
, toBeDeleted_(false) // deletion flag
, sid_ (sid) // System-specific body index
, uid_ (uid) // User-specific body ID
+ , gpos_() // Global position of the center of mass
+ , q_() // Orientation of the body frame
+ , R_() // Rigid body rotation
, typeID_(typeID) // geometry type
{
sb_ = this; // The superordinate rigid body
diff --git a/src/pe/rigidbody/RigidBody.h b/src/pe/rigidbody/RigidBody.h
index d987ca4e7..1ac988c75 100644
--- a/src/pe/rigidbody/RigidBody.h
+++ b/src/pe/rigidbody/RigidBody.h
@@ -98,14 +98,15 @@ public:
inline bool isMarkedForDeletion() const { return toBeDeleted_; }
inline id_t getSystemID() const;
inline id_t getID() const;
- inline const Vec3& getRelPosition() const;
- inline const Vec3& getPosition() const;
- inline const Vec3 getRelLinearVel() const;
- inline const Vec3& getLinearVel() const;
- inline const Vec3 getRelAngularVel() const;
- inline const Vec3& getAngularVel() const;
- inline const Quat& getQuaternion() const;
- inline const Mat3& getRotation() const;
+ inline const Vec3 getRelPosition() const;
+ inline const Vec3 getPosition() const;
+ inline const Vec3 getBodyLinearVel() const;
+ inline const Vec3 getLinearVel() const;
+ inline const Vec3 getBodyAngularVel() const;
+ inline const Vec3 & getAngularVel() const;
+ inline const Quat getRelQuaternion() const;
+ inline const Quat getQuaternion() const;
+ inline const Mat3 getRotation() const;
inline real_t getMass() const;
inline real_t getInvMass() const;
virtual inline real_t getVolume() const;
@@ -149,8 +150,10 @@ public:
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 setRelPosition( 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 setRelOrientation( const Quat& q );
inline void setMassAndInertiaToInfinity();
inline void setRelLinearVel ( real_t vx, real_t vy, real_t vz );
@@ -309,8 +312,6 @@ protected:
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();
//@}
//**********************************************************************************************
@@ -321,14 +322,6 @@ protected:
//@}
//**********************************************************************************************
- //**Simulation functions************************************************************************
- /*!\name Simulation functions */
- //@{
- virtual void update( const Vec3& dp ); // Translation update of a subordinate rigid body
- virtual void update( const Quat& dq ); // Rotation update of a subordinate rigid body
- //@}
- //**********************************************************************************************
-
//**Functions for internal changes in compound geometries***************************************
/*!\name Functions for internal changes in compound geometries */
//@{
@@ -357,12 +350,7 @@ protected:
real_t motion_; //!< The current motion of the rigid body.
/*!< If this value drops under the specified sleep threshold, the
rigid body will be put to sleep. */
- Vec3 gpos_; //!< The global position of the center of mass of this rigid body.
- Vec3 rpos_; //!< The relative position within the body frame of the superordinate body.
- /*!< If the body is contained within a superordinate Union the relative
- position gives the position of the body's center of mass within the
- body frame of the superordinate body. If the body is not contained
- in a Union, the relative position is 0. */
+
mutable Vec3 v_; //!< The linear velocity of this rigid body.
mutable Vec3 w_; //!< Angular velocity of this rigid body.
Vec3 force_; //!< Total force (external+contact) acting in the body's center of mass.
@@ -379,8 +367,7 @@ protected:
I_{zx} & I_{zy} & I_{zz} \\
\end{array}\right)\f] */
Mat3 Iinv_; //!< The inverse moment of inertia within the body frame.
- Quat q_; //!< The orientation of the body frame in the global world frame.
- Mat3 R_; //!< The rotation in reference to the global frame of reference.
+
//@}
//**********************************************************************************************
@@ -390,6 +377,7 @@ protected:
BodyID sb_; //!< The superordinate rigid body.
/*!< This data member is the connection to the superordinate body,
which is either the enclosing Union or the rigid body itself. */
+
bool finite_; //!< Finiteness flag.
/*!< The flag value indicates if the rigid body is finite (\a true)
or infinite (\a false). */
@@ -417,7 +405,15 @@ protected:
AABB aabb_; //!< Axis-aligned bounding box for the rigid body.
//@}
//**********************************************************************************************
+
private:
+
+ Vec3 gpos_; /*!< The position of the center of mass of this rigid body. If the body is contained in union
+ this and the other properties are relative to its center / orientation. Use the respective function
+ getPosition() / getRotation() to access the actual global Position in the world frame */
+ Quat q_; //!< The orientation of the body frame in the global world frame.
+ Mat3 R_; //!< The rotation in reference to the global frame of reference.
+
id_t typeID_; //< identify geometry type
};
@@ -694,29 +690,38 @@ inline id_t RigidBody::getID() const
//*************************************************************************************************
+
+
+
//*************************************************************************************************
-/*!\brief Returns the relative position of the rigid body within the superordinate body.
- *
- * \return The relative position of the rigid body.
+/*!\brief Returns the global position of the center of mass of the rigid body.
*
- * If the rigid body is not contained in a superordinate body, the returned relative position will
- * be \f$ \left(\begin{array}{*{3}{c}} 0 & 0 & 0 \end{array}\right) \f$.
+ * \return The global position of the center of mass.
*/
-inline const Vec3& RigidBody::getRelPosition() const
+inline const Vec3 RigidBody::getPosition() const
{
- return rpos_;
+ if( !hasSuperBody() ){
+ return gpos_;
+ }
+ return sb_->getPosition() + sb_->getRotation() * getRelPosition();
+
}
//*************************************************************************************************
-
//*************************************************************************************************
-/*!\brief Returns the global position of the center of mass of the rigid body.
+/*!\brief Returns the relative position of the rigid body within the superordinate body.
*
- * \return The global position of the center of mass.
+ * \return The relative position of the rigid body.
+ *
+ * If this body is not containted in a superbody the zero-vector is returned.
*/
-inline const Vec3& RigidBody::getPosition() const
+inline const Vec3 RigidBody::getRelPosition() const
{
- return gpos_;
+ if( hasSuperBody() ){
+ return gpos_;
+ }
+ return Vec3(0, 0, 0);
+
}
//*************************************************************************************************
@@ -729,11 +734,9 @@ inline const Vec3& RigidBody::getPosition() const
* This function returns the linear velocity of the center of mass of the rigid body in reference
* to the body's own frame of reference.
*/
-inline const Vec3 RigidBody::getRelLinearVel() const
+inline const Vec3 RigidBody::getBodyLinearVel() const
{
- if( hasSuperBody() )
- v_ = sb_->velFromWF( gpos_ );
- return R_.getTranspose() * v_;
+ return getRotation().getTranspose() * getLinearVel();
}
//*************************************************************************************************
@@ -746,11 +749,13 @@ inline const Vec3 RigidBody::getRelLinearVel() const
* This function returns the linear velocity of the center of mass of the rigid body in reference
* to the global world frame.
*/
-inline const Vec3& RigidBody::getLinearVel() const
+inline const Vec3 RigidBody::getLinearVel() const
{
- if( hasSuperBody() )
- v_ = sb_->velFromWF( gpos_ );
+ if( hasSuperBody() ){
+ return sb_->getLinearVel() + (sb_->getAngularVel() % (getPosition() - sb_->getPosition()));
+ }
return v_;
+
}
//*************************************************************************************************
@@ -760,14 +765,12 @@ inline const Vec3& RigidBody::getLinearVel() const
*
* \return The relative angular velocity.
*
- * This function returns the angluar velocity of the center of mass in reference to the body's
+ * This function returns the angular velocity of the center of mass in reference to the body's
* own frame of reference.
*/
-inline const Vec3 RigidBody::getRelAngularVel() const
+inline const Vec3 RigidBody::getBodyAngularVel() const
{
- if( hasSuperBody() )
- w_ = sb_->getAngularVel();
- return R_.getTranspose() * w_;
+ return getRotation().getTranspose() * getAngularVel();
}
//*************************************************************************************************
@@ -780,10 +783,11 @@ inline const Vec3 RigidBody::getRelAngularVel() const
* This function returns the angluar velocity of the center of mass in reference to the global
* world frame.
*/
-inline const Vec3& RigidBody::getAngularVel() const
+inline const Vec3 & RigidBody::getAngularVel() const
{
- if( hasSuperBody() )
- w_ = sb_->getAngularVel();
+ if( hasSuperBody() ) {
+ return sb_->getAngularVel();
+ }
return w_;
}
//*************************************************************************************************
@@ -797,13 +801,34 @@ inline const Vec3& RigidBody::getAngularVel() const
* This function returns the quaternion of the rigid body, which represents the orientation of
* the body in reference to the global world frame.
*/
-inline const Quat& RigidBody::getQuaternion() const
+inline const Quat RigidBody::getQuaternion() const
{
+ if( hasSuperBody() ){
+ return getRelQuaternion() * sb_->getQuaternion();
+ }
return q_;
}
//*************************************************************************************************
+//*************************************************************************************************
+/*!\brief Returns the relative orientation of the rigid body to its superbody.
+ *
+ * \return The relative orientation of the rigid body.
+ *
+ * This function returns the quaternion, which represents the orientation of
+ * the body in reference to coordinate system of the superbody. If this body has no superbody
+ * the union-quaternion is returned.
+ */
+inline const Quat RigidBody::getRelQuaternion() const
+{
+ if( hasSuperBody() ){
+ return q_;
+ }
+ return Quat();
+}
+//*************************************************************************************************
+
//*************************************************************************************************
/*!\brief Returns the rotation of the rigid body.
*
@@ -812,8 +837,11 @@ inline const Quat& RigidBody::getQuaternion() const
* This function returns the rotation matrix of the rigid body, which represents the rotation of
* the body in reference to the global world frame.
*/
-inline const Mat3& RigidBody::getRotation() const
+inline const Mat3 RigidBody::getRotation() const
{
+ if( hasSuperBody() ) {
+ return getQuaternion().toRotationMatrix();
+ }
return R_;
}
//*************************************************************************************************
@@ -874,7 +902,7 @@ inline const Mat3& RigidBody::getBodyInertia() const
*/
inline const Mat3 RigidBody::getInertia() const
{
- return math::transformMatrixRART(R_, I_);
+ return math::transformMatrixRART(getRotation(), I_);
}
//*************************************************************************************************
@@ -898,7 +926,7 @@ inline const Mat3& RigidBody::getInvBodyInertia() const
*/
inline const Mat3 RigidBody::getInvInertia() const
{
- return math::transformMatrixRART(R_, Iinv_);
+ return math::transformMatrixRART(getRotation(), Iinv_);
}
//*************************************************************************************************
@@ -941,6 +969,7 @@ inline real_t RigidBody::getKineticEnergy() const
return real_c(0.5) * getMass() * getLinearVel() * getLinearVel();
}
//*************************************************************************************************
+
//*************************************************************************************************
/*!\brief Returns the rotational energy of the rigid body
*
@@ -951,6 +980,7 @@ inline real_t RigidBody::getRotationalEnergy() const
return real_c(0.5) * getAngularVel() * (getInertia() * getAngularVel());
}
//*************************************************************************************************
+
//*************************************************************************************************
/*!\brief Returns the energy of the rigid body
*
@@ -976,7 +1006,7 @@ inline real_t RigidBody::getEnergy() const
*/
inline const Vec3 RigidBody::vectorFromBFtoWF( real_t vx, real_t vy, real_t vz ) const
{
- return R_ * Vec3( vx, vy, vz );
+ return getRotation() * Vec3( vx, vy, vz );
}
//*************************************************************************************************
@@ -992,7 +1022,7 @@ inline const Vec3 RigidBody::vectorFromBFtoWF( real_t vx, real_t vy, real_t vz )
*/
inline const Vec3 RigidBody::vectorFromBFtoWF( const Vec3& v ) const
{
- return R_ * v;
+ return getRotation() * v;
}
//*************************************************************************************************
@@ -1010,7 +1040,7 @@ inline const Vec3 RigidBody::vectorFromBFtoWF( const Vec3& v ) const
*/
inline const Vec3 RigidBody::pointFromBFtoWF( real_t px, real_t py, real_t pz ) const
{
- return gpos_ + ( R_ * Vec3( px, py, pz ) );
+ return getPosition() + ( getRotation() * Vec3( px, py, pz ) );
}
//*************************************************************************************************
@@ -1026,7 +1056,7 @@ inline const Vec3 RigidBody::pointFromBFtoWF( real_t px, real_t py, real_t pz )
*/
inline const Vec3 RigidBody::pointFromBFtoWF( const Vec3& rpos ) const
{
- return gpos_ + ( R_ * rpos );
+ return getPosition() + ( getRotation() * rpos );
}
//*************************************************************************************************
@@ -1044,7 +1074,7 @@ inline const Vec3 RigidBody::pointFromBFtoWF( const Vec3& rpos ) const
*/
inline const Vec3 RigidBody::vectorFromWFtoBF( real_t vx, real_t vy, real_t vz ) const
{
- return R_.getTranspose() * Vec3( vx, vy, vz );
+ return getRotation() .getTranspose() * Vec3( vx, vy, vz );
}
//*************************************************************************************************
@@ -1060,7 +1090,7 @@ inline const Vec3 RigidBody::vectorFromWFtoBF( real_t vx, real_t vy, real_t vz )
*/
inline const Vec3 RigidBody::vectorFromWFtoBF( const Vec3& v ) const
{
- return R_.getTranspose() * v;
+ return getRotation() .getTranspose() * v;
}
//*************************************************************************************************
@@ -1078,7 +1108,7 @@ inline const Vec3 RigidBody::vectorFromWFtoBF( const Vec3& v ) const
*/
inline const Vec3 RigidBody::pointFromWFtoBF( real_t px, real_t py, real_t pz ) const
{
- return R_.getTranspose() * ( Vec3( px, py, pz ) - gpos_ );
+ return getRotation().getTranspose() * ( Vec3( px, py, pz ) - getPosition() );
}
//*************************************************************************************************
@@ -1094,7 +1124,7 @@ inline const Vec3 RigidBody::pointFromWFtoBF( real_t px, real_t py, real_t pz )
*/
inline const Vec3 RigidBody::pointFromWFtoBF( const Vec3& gpos ) const
{
- return R_.getTranspose() * ( gpos - gpos_ );
+ return getRotation().getTranspose() * ( gpos - getPosition() );
}
//*************************************************************************************************
@@ -1279,23 +1309,6 @@ inline void RigidBody::setAngularVel( const Vec3& avel )
//*************************************************************************************************
-//*************************************************************************************************
-/*!\brief Calculation of the relative position within a superordinate body.
- *
- * \return void
- *
- * The function calculates the relative position depending on its current global position, the
- * current global position of the superordinate body and the rotation of the superordinate body.
- */
-inline void RigidBody::calcRelPosition()
-{
- rpos_ = sb_->R_.getTranspose()*( gpos_ - sb_->gpos_ );
-}
-//*************************************************************************************************
-
-
-
-
//=================================================================================================
//
// FORCE FUNCTIONS
@@ -1374,8 +1387,12 @@ inline void RigidBody::resetSB()
inline void RigidBody::setForce( const Vec3& f )
{
// Increasing the force on this rigid body
+ if(hasSuperBody()){
+ sb_->setForce(f);
+ sb_->setTorque((getPosition() - sb_->getPosition()) % f);
+ return;
+ }
force_ = f;
-
wake();
}
//*************************************************************************************************
@@ -1389,9 +1406,13 @@ inline void RigidBody::setForce( const Vec3& f )
*/
inline void RigidBody::setTorque( const Vec3& tau )
{
- // Increasing the force on this rigid body
+ // Increasing the torque on this rigid body
+ // Increasing the force on this rigid body
+ if(hasSuperBody()){
+ sb_->setTorque(tau);
+ return;
+ }
torque_ = tau;
-
wake();
}
//*************************************************************************************************
@@ -1479,7 +1500,7 @@ inline void RigidBody::addForce( const Vec3& f )
wake();
} else
{
- sb_->addForceAtPos( f, gpos_ );
+ sb_->addForceAtPos( f, getPosition() );
}
}
//*************************************************************************************************
@@ -2099,6 +2120,7 @@ inline void RigidBody::setVisible( bool visible )
//*************************************************************************************************
inline void RigidBody::setPositionImpl( real_t px, real_t py, real_t pz )
{
+ if ( hasSuperBody() ) return;
if (isFixed())
WALBERLA_ABORT("Trying to move a fixed body: " << *this);
@@ -2112,6 +2134,36 @@ inline void RigidBody::setPositionImpl( real_t px, real_t py, real_t pz )
}
//*************************************************************************************************
+//*************************************************************************************************
+/*!\fn void RigidBody::setRelPosition( const Vec3 &gpos )
+ * \brief Setting the global position of the rigid body.
+ *
+ * \param gpos The relative position.
+ * \return void
+ *
+ * This function sets the relatvie position of the rigid body w. r. t. the superbodies position.
+ *
+ * \b Note:
+ * - Setting the position of a rigid body without superbody will have no effect.
+ * - 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 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!
+ */
+//*************************************************************************************************
+inline void RigidBody::setRelPosition(const Vec3 &gpos) {
+
+ if ( !hasSuperBody() ) return;
+ if (isFixed())
+ WALBERLA_ABORT("Trying to move a fixed body: " << *this);
+
+ gpos_ = gpos;
+ calcBoundingBox(); // Updating the axis-aligned bounding box of the box
+ wake(); // Waking the box from sleep mode
+}
+//*************************************************************************************************
+
//*************************************************************************************************
/*!\fn void RigidBody::setPosition( real_t px, real_t py, real_t pz )
* \brief Setting the global position of the rigid body.
@@ -2194,6 +2246,7 @@ inline void RigidBody::setPosition( const Vec3& gpos )
//*************************************************************************************************
inline void RigidBody::setOrientationImpl( real_t r, real_t i, real_t j, real_t k )
{
+ if ( hasSuperBody() ) return;
if (isFixed())
WALBERLA_ABORT("Trying to rotate a fixed body: " << *this);
@@ -2263,6 +2316,38 @@ inline void RigidBody::setOrientation( const Quat& q )
}
//*************************************************************************************************
+//*************************************************************************************************
+/*!\fn void RigidBody::setRelOrientation( const Quat &q )
+ * \brief Setting the relative position of the rigid body.
+ *
+ * \param q The relative orientation.
+ * \return void
+ *
+ * This function sets the relative orientation of the rigid body w. r. t. the superbodies position.
+ *
+ * \b Note:
+ * - Setting the relative orientation of a rigid body without superbody will have no effect.
+ * - 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 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!
+ */
+//*************************************************************************************************
+inline void RigidBody::setRelOrientation(const Quat &q) {
+
+ if ( !hasSuperBody() ) return;
+
+ if (isFixed())
+ WALBERLA_ABORT("Trying to move a fixed body: " << *this);
+
+ q_ = q;
+
+ calcBoundingBox(); // Updating the axis-aligned bounding box of the box
+ wake(); // Waking the box from sleep mode
+}
+//*************************************************************************************************
+
//*************************************************************************************************
/*!\fn void RigidBody::setMassAndInertiaToInfinity( )
* \brief Setting the mass to infinity. This will also make the inertia tensor infinite.
@@ -2990,125 +3075,6 @@ inline void RigidBody::fix()
//*************************************************************************************************
-//=================================================================================================
-//
-// SIMULATION FUNCTIONS
-//
-//=================================================================================================
-
-//*************************************************************************************************
-/*!\brief Translation update of a subordinate rigid body.
- *
- * \param dp Change in the global position of the superordinate rigid body.
- * \return void
- *
- * This update function is triggered by the superordinate body in case of a translational
- * movement. This movement involves a change in the global position and the axis-aligned
- * bounding box.
- */
-inline void RigidBody::update( const Vec3& dp )
-{
- // Checking the state of the sphere
- WALBERLA_ASSERT( checkInvariants(), "Invalid sphere state detected" );
- WALBERLA_ASSERT( hasSuperBody(), "Invalid superordinate body detected" );
-
- // Updating the global position
- gpos_ += dp;
-
- // Setting the axis-aligned bounding box
- calcBoundingBox();
-
- // Checking the state of the sphere
- WALBERLA_ASSERT( checkInvariants(), "Invalid sphere state detected" );
-}
-//*************************************************************************************************
-
-
-//*************************************************************************************************
-/*!\brief Rotation update of a subordinate rigid body.
- *
- * \param dq Change in the orientation of the superordinate rigid body.
- * \return void
- *
- * This update function is triggered by the superordinate body in case of a rotational movement.
- * This movement involves a change in the global position, the orientation/rotation and the
- * axis-aligned bounding box of the sphere.
- */
-inline void RigidBody::update( const Quat& dq )
-{
- // Checking the state of the sphere
- WALBERLA_ASSERT( checkInvariants(), "Invalid sphere state detected" );
- WALBERLA_ASSERT( hasSuperBody(), "Invalid superordinate body detected" );
-
- // Calculating the new global position
- gpos_ = sb_->getPosition() + ( sb_->getRotation() * rpos_ );
-
- // Calculating the new orientation and rotation
- q_ = dq * q_;
- R_ = q_.toRotationMatrix();
-
- // Setting the axis-aligned bounding box
- calcBoundingBox();
-
- // Checking the state of the sphere
- WALBERLA_ASSERT( checkInvariants(), "Invalid sphere state detected" );
-}
-//*************************************************************************************************
-
-
-
-
-//=================================================================================================
-//
-// SIMULATION FUNCTIONS
-//
-//=================================================================================================
-
-//*************************************************************************************************
-/*!\fn void RigidBody::update( const Vec3& dp )
- * \brief Translation update of a subordinate rigid body.
- *
- * \param dp Change in the global position of the superordinate rigid body.
- * \return void
- *
- * This function calculates the necessary updates for a subordinate rigid body contained
- * in a superordinate rigid body that has performed a translational movement. This function
- * is triggered automatically by the superordinate body in case of a translational movement.
- * All classes deriving from the RigidBody class have to implement this function to update
- * the properties of the rigid body. The following properties of the body might change due
- * to this translation. All derived classes have to make sure these properties are updated
- * correctly:
- *
- * - the global position
- * - the axis-aligned bounding box
- */
-//*************************************************************************************************
-
-
-//*************************************************************************************************
-/*!\fn void RigidBody::update( const Quat& dq )
- * \brief Rotation update of a subordinate rigid body.
- *
- * \param dq Change in the orientation of the superordinate rigid body.
- * \return void
- *
- * This function calculates the necessary updates for a subordinate rigid body contained
- * in a superordinate rigid body that has performed a rotational movement. The function
- * is triggered automatically by the superordinate body in case of a rotational movement.
- * All classes deriving from the RigidBody class have to implement this function to update
- * the properties of the rigid body. The following properties of the body might change due
- * to this rotation. All derived classes have to make sure these properties are updated
- * correctly:
- *
- * - the global position
- * - the orientation/rotation (i.e. the quaterion and the rotation matrix)
- * - the axis-aligned bounding box
- */
-//*************************************************************************************************
-
-
-
-
//=================================================================================================
//
// FUNCTIONS FOR INTERNAL CHANGES IN COMPOUND GEOMETRIES
diff --git a/src/pe/rigidbody/Sphere.cpp b/src/pe/rigidbody/Sphere.cpp
index 2d45389cb..c3a76c6d8 100644
--- a/src/pe/rigidbody/Sphere.cpp
+++ b/src/pe/rigidbody/Sphere.cpp
@@ -50,7 +50,6 @@ namespace pe {
* \param sid Unique system-specific ID for the sphere.
* \param uid User-specific ID for the sphere.
* \param gpos Global geometric center of the sphere.
- * \param rpos The relative position within the body frame of a superordinate body.
* \param q The orientation of the sphere's body frame in the global world frame.
* \param radius The radius of the sphere \f$ (0..\infty) \f$.
* \param material The material of the sphere.
@@ -58,12 +57,13 @@ namespace pe {
* \param communicating specifies if the sphere should take part in synchronization (syncNextNeighbour, syncShadowOwner)
* \param infiniteMass specifies if the sphere has infinite mass and will be treated as an obstacle
*/
-Sphere::Sphere( id_t sid, id_t uid, const Vec3& gpos, const Vec3& rpos, const Quat& q,
+Sphere::Sphere( id_t sid, id_t uid, const Vec3& gpos, const Quat& q,
real_t radius, MaterialID material,
const bool global, const bool communicating, const bool infiniteMass )
- : Sphere::Sphere( getStaticTypeID(), sid, uid, gpos, rpos, q, radius, material, global, communicating, infiniteMass )
+ : Sphere::Sphere( getStaticTypeID(), sid, uid, gpos, q, radius, material, global, communicating, infiniteMass )
{}
-Sphere::Sphere( id_t const typeId, id_t sid, id_t uid, const Vec3& gpos, const Vec3& rpos, const Quat& q,
+
+Sphere::Sphere( id_t const typeId, id_t sid, id_t uid, const Vec3& gpos, const Quat& q,
real_t radius, MaterialID material,
const bool global, const bool communicating, const bool infiniteMass )
: GeomPrimitive( typeId, sid, uid, material ) // Initialization of the parent class
@@ -76,12 +76,8 @@ Sphere::Sphere( id_t const typeId, id_t sid, id_t uid, const Vec3& gpos, const V
WALBERLA_ASSERT( radius > real_c(0), "Invalid sphere radius" );
// Setting the center of the sphere
- gpos_ = gpos;
-
- // Initializing the instantiated sphere
- 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)
@@ -196,11 +192,12 @@ void Sphere::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"
diff --git a/src/pe/rigidbody/Sphere.h b/src/pe/rigidbody/Sphere.h
index 8efc9bcae..c4a1c95fc 100644
--- a/src/pe/rigidbody/Sphere.h
+++ b/src/pe/rigidbody/Sphere.h
@@ -63,10 +63,10 @@ public:
//**Constructors********************************************************************************
/*!\name Constructors */
//@{
- explicit Sphere( id_t sid, id_t uid, const Vec3& gpos, const Vec3& rpos, const Quat& q,
+ explicit Sphere( id_t sid, id_t uid, const Vec3& gpos, const Quat& q,
real_t radius, MaterialID material,
const bool global, const bool communicating, const bool infiniteMass );
- explicit Sphere( id_t const typeID, id_t sid, id_t uid, const Vec3& gpos, const Vec3& rpos, const Quat& q,
+ explicit Sphere( id_t const typeID, id_t sid, id_t uid, const Vec3& gpos, const Quat& q,
real_t radius, MaterialID material,
const bool global, const bool communicating, const bool infiniteMass );
//@}
@@ -115,7 +115,6 @@ public:
/*!\name Utility functions */
//@{
inline virtual Vec3 support( const Vec3& d ) const;
- inline virtual Vec3 supportContactThreshold( const Vec3& d ) const;
//@}
//**********************************************************************************************
@@ -266,18 +265,18 @@ inline real_t Sphere::calcDensity( real_t radius, real_t mass )
inline void Sphere::calcBoundingBox()
{
const real_t length( radius_ + contactThreshold );
-
+ Vec3 gpos = getPosition();
aabb_ = AABB(
- gpos_[0] - length,
- gpos_[1] - length,
- gpos_[2] - length,
- gpos_[0] + length,
- gpos_[1] + length,
- gpos_[2] + length
+ gpos[0] - length,
+ gpos[1] - length,
+ gpos[2] - length,
+ gpos[0] + length,
+ gpos[1] + length,
+ gpos[2] + length
);
// WALBERLA_ASSERT( aabb_.isValid() , "Invalid bounding box detected" );
- WALBERLA_ASSERT( aabb_.contains( gpos_ ), "Invalid bounding box detected("<< getSystemID() <<")\n" << "pos: " << gpos_ << "\nlength: " << length << "\nvel: " << getLinearVel() << "\nbox: " << aabb_ );
+ WALBERLA_ASSERT( aabb_.contains( gpos ), "Invalid bounding box detected("<< getSystemID() <<")\n" << "pos: " << gpos << "\nlength: " << length << "\nvel: " << getLinearVel() << "\nbox: " << aabb_ );
}
//*************************************************************************************************
@@ -317,29 +316,6 @@ inline Vec3 Sphere::support( const Vec3& d ) const
//*************************************************************************************************
-//*************************************************************************************************
-/*!\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 Sphere::supportContactThreshold( const Vec3& d ) const
-{
- auto len = d.sqrLength();
- if (!math::equal(len, real_t(0)))
- {
- //WALBERLA_ASSERT_FLOAT_EQUAL( len, real_t(1), "search direction not normalized!");
- const Vec3 s = getPosition() + (getRadius() / sqrt(len) + contactThreshold) * d;
- //std::cout << "Support in direction " << d << " with center " << getPosition() << " (r=" << getRadius() << ") is " << s << std::endl;
- return s;
- } else
- {
- return Vec3(0,0,0);
- }
-}
-//*************************************************************************************************
-
//=================================================================================================
//
// UTILITY FUNCTIONS
@@ -409,7 +385,7 @@ inline real_t Sphere::getDepth( real_t px, real_t py, real_t pz ) const
*/
inline real_t Sphere::getDepth( const Vec3& gpos ) const
{
- return ( radius_ - ( gpos - gpos_ ).length() );
+ return ( radius_ - ( gpos - getPosition() ).length() );
}
//*************************************************************************************************
@@ -477,7 +453,7 @@ inline real_t Sphere::getDistance( real_t px, real_t py, real_t pz ) const
*/
inline real_t Sphere::getDistance( const Vec3& gpos ) const
{
- return ( ( gpos - gpos_ ).length() - radius_ );
+ return ( ( gpos - getPosition() ).length() - radius_ );
}
//*************************************************************************************************
diff --git a/src/pe/rigidbody/SphereFactory.cpp b/src/pe/rigidbody/SphereFactory.cpp
index d88774c45..dfe500daa 100644
--- a/src/pe/rigidbody/SphereFactory.cpp
+++ b/src/pe/rigidbody/SphereFactory.cpp
@@ -49,7 +49,7 @@ SphereID createSphere( BodyStorage& globalStorage, BlockStorage& blocks, BlockDa
const id_t sid = UniqueID<RigidBody>::createGlobal();
WALBERLA_ASSERT_EQUAL(communicating, false);
WALBERLA_ASSERT_EQUAL(infiniteMass, true);
- SpherePtr sp = std::make_unique<Sphere>(sid, uid, gpos, Vec3(0,0,0), Quat(), radius, material, global, false, true);
+ SpherePtr sp = std::make_unique<Sphere>(sid, uid, gpos, Quat(), radius, material, global, false, true);
sphere = static_cast<SphereID>(&globalStorage.add( std::move(sp) ));
} else
{
@@ -59,7 +59,7 @@ SphereID createSphere( BodyStorage& globalStorage, BlockStorage& blocks, BlockDa
const id_t sid( UniqueID<RigidBody>::create() );
BodyStorage& bs = (*block.getData<Storage>(storageID))[0];
- SpherePtr sp = std::make_unique<Sphere>(sid, uid, gpos, Vec3(0,0,0), Quat(), radius, material, global, communicating, infiniteMass);
+ SpherePtr sp = std::make_unique<Sphere>(sid, uid, gpos, Quat(), radius, material, global, communicating, infiniteMass);
sp->MPITrait.setOwner(Owner(MPIManager::instance()->rank(), block.getId().getID()));
sphere = static_cast<SphereID>(&bs.add( std::move(sp) ));
}
diff --git a/src/pe/rigidbody/Squirmer.cpp b/src/pe/rigidbody/Squirmer.cpp
index 99f650b00..e9d386c59 100644
--- a/src/pe/rigidbody/Squirmer.cpp
+++ b/src/pe/rigidbody/Squirmer.cpp
@@ -48,10 +48,10 @@ namespace pe {
* \param communicating specifies if the sphere should take part in synchronization (syncNextNeighbour, syncShadowOwner)
* \param infiniteMass specifies if the sphere has infinite mass and will be treated as an obstacle
*/
-Squirmer::Squirmer( id_t sid, id_t uid, const Vec3& gpos, const Vec3& rpos, const Quat& q,
+Squirmer::Squirmer( id_t sid, id_t uid, const Vec3& gpos, const Quat& q,
real_t radius, real_t squirmerVelocity, real_t squirmerBeta, MaterialID material,
const bool global, const bool communicating, const bool infiniteMass )
- : Sphere( getStaticTypeID(), sid, uid, gpos, rpos, q, radius, material, global, communicating, infiniteMass ) // Initialization of the parent class
+ : Sphere( getStaticTypeID(), sid, uid, gpos, q, radius, material, global, communicating, infiniteMass ) // Initialization of the parent class
, squirmerVelocity_(squirmerVelocity), squirmerBeta_(squirmerBeta)
{
}
@@ -106,7 +106,7 @@ const Vec3 Squirmer::velFromBF( real_t px, real_t py, real_t pz ) const
*/
const Vec3 Squirmer::velFromBF( const Vec3& rpos ) const
{
- return Sphere::velFromBF( rpos ) + getSquirmerVelocity( R_ * rpos );
+ return Sphere::velFromBF( rpos ) + getSquirmerVelocity( getRotation() * rpos );
}
//*************************************************************************************************
@@ -154,7 +154,7 @@ const Vec3 Squirmer::getSquirmerVelocity( const Vec3& rpos ) const
*/
const Vec3 Squirmer::velFromWF( const Vec3& gpos ) const
{
- return Sphere::velFromWF( gpos ) + getSquirmerVelocity( gpos - gpos_ );
+ return Sphere::velFromWF( gpos ) + getSquirmerVelocity( gpos - getPosition() );
}
//*************************************************************************************************
diff --git a/src/pe/rigidbody/Squirmer.h b/src/pe/rigidbody/Squirmer.h
index 5dcf205c4..0d9514324 100644
--- a/src/pe/rigidbody/Squirmer.h
+++ b/src/pe/rigidbody/Squirmer.h
@@ -32,7 +32,7 @@ public:
//**Constructors********************************************************************************
/*!\name Constructors */
//@{
- explicit Squirmer( id_t sid, id_t uid, const Vec3& gpos, const Vec3& rpos, const Quat& q,
+ explicit Squirmer( id_t sid, id_t uid, const Vec3& gpos, const Quat& q,
real_t radius, real_t squirmerVelocity, real_t squirmerBeta, MaterialID material,
const bool global, const bool communicating, const bool infiniteMass );
//@}
diff --git a/src/pe/rigidbody/SquirmerFactory.cpp b/src/pe/rigidbody/SquirmerFactory.cpp
index 7160ff2d0..caac4c9a0 100644
--- a/src/pe/rigidbody/SquirmerFactory.cpp
+++ b/src/pe/rigidbody/SquirmerFactory.cpp
@@ -48,7 +48,7 @@ SquirmerID createSquirmer( BodyStorage& globalStorage, BlockStorage& blocks, Blo
const id_t sid = UniqueID<RigidBody>::createGlobal();
WALBERLA_ASSERT_EQUAL(communicating, false);
WALBERLA_ASSERT_EQUAL(infiniteMass, true);
- SquirmerPtr sq = std::make_unique<Squirmer>(sid, uid, gpos, Vec3(0,0,0), Quat(), radius, squirmerVelocity, squirmerBeta, material, global, false, true);
+ SquirmerPtr sq = std::make_unique<Squirmer>(sid, uid, gpos, Quat(), radius, squirmerVelocity, squirmerBeta, material, global, false, true);
squirmer = static_cast<SquirmerID>(&globalStorage.add(std::move(sq)));
} else
{
@@ -58,7 +58,7 @@ SquirmerID createSquirmer( BodyStorage& globalStorage, BlockStorage& blocks, Blo
const id_t sid( UniqueID<RigidBody>::create() );
BodyStorage& bs = (*block.getData<Storage>(storageID))[0];
- SquirmerPtr sq = std::make_unique<Squirmer>(sid, uid, gpos, Vec3(0,0,0), Quat(), radius, squirmerVelocity, squirmerBeta, material, global, communicating, infiniteMass);
+ SquirmerPtr sq = std::make_unique<Squirmer>(sid, uid, gpos, Quat(), radius, squirmerVelocity, squirmerBeta, material, global, communicating, infiniteMass);
sq->MPITrait.setOwner(Owner(MPIManager::instance()->rank(), block.getId().getID()));
squirmer = static_cast<SquirmerID>(&bs.add( std::move(sq) ));
}
diff --git a/src/pe/rigidbody/Union.h b/src/pe/rigidbody/Union.h
index 3311851f5..2f462877f 100644
--- a/src/pe/rigidbody/Union.h
+++ b/src/pe/rigidbody/Union.h
@@ -55,11 +55,11 @@ namespace pe {
//*************************************************************************************************
/**
* \ingroup pe
- * \brief Base class for the sphere geometry.
+ * \brief Base class for the union geometry, a rigid assebly of bodies.
*
- * The Sphere class represents the base class for the sphere geometry. It provides
- * the basic functionality of a sphere. For a full description of the sphere geometry,
- * see the Sphere class description.
+ * The Union class represents the base class for the sphere geometry.
+ * For a full description of the union geometry,
+ * see the Union class description.
*/
template <typename... BodyTypes>
class Union : public RigidBody
@@ -79,7 +79,7 @@ public:
//**Constructors********************************************************************************
/*!\name Constructors */
//@{
- explicit Union( id_t sid, id_t uid, const Vec3& gpos, const Vec3& rpos, const Quat& q,
+ explicit Union( id_t sid, id_t uid, const Vec3& gpos, const Quat& q,
const bool global, const bool communicating, const bool infiniteMass );
//@}
//**********************************************************************************************
@@ -132,14 +132,6 @@ public:
virtual inline bool hasSubBodies() const WALBERLA_OVERRIDE { return true; }
- //**Simulation functions************************************************************************
- /*!\name Simulation functions */
- //@{
- virtual void update( const Vec3& dp ) WALBERLA_OVERRIDE; // Translation update of a subordinate rigid body
- virtual void update( const Quat& dq ) WALBERLA_OVERRIDE; // Rotation update of a subordinate rigid body
- //@}
- //**********************************************************************************************
-
//**Signal functions***************************************************************************
/*!\name Signal functions */
//@{
@@ -189,7 +181,7 @@ protected:
/*!\name Utility functions */
//@{
inline virtual void calcBoundingBox() WALBERLA_OVERRIDE; // Calculation of the axis-aligned bounding box
- void calcCenterOfMass(); ///< updates the center of mass (gpos)
+ inline void calcCenterOfMass(); // Compute mass and center of gravity
inline void calcInertia(); // Calculation of the moment of inertia
//@}
//**********************************************************************************************
@@ -228,15 +220,13 @@ private:
* \param infiniteMass specifies if the sphere has infinite mass and will be treated as an obstacle
*/
template <typename... BodyTypes>
-Union<BodyTypes...>::Union( id_t sid, id_t uid, const Vec3& gpos, const Vec3& rpos, const Quat& q,
+Union<BodyTypes...>::Union( id_t sid, id_t uid, const Vec3& gpos, const Quat& q,
const bool global, const bool communicating, const bool /*infiniteMass*/ )
: RigidBody( getStaticTypeID(), sid, uid ) // Initialization of the parent class
{
// Initializing the instantiated union
- gpos_ = gpos; // Setting the global center of mass
- rpos_ = rpos; // Setting the relative position
- q_ = q; // Setting the orientation
- R_ = q_.toRotationMatrix(); // Setting the rotation matrix
+ setPosition(gpos); // Setting the global center of mass
+ setOrientation(q); // Setting the orientation
calcCenterOfMass();
calcInertia();
@@ -258,7 +248,7 @@ Union<BodyTypes...>::Union( id_t sid, id_t uid, const Vec3& gpos, const Vec3& rp
//=================================================================================================
//*************************************************************************************************
-/*!\brief Destructor for the Sphere class.
+/*!\brief Destructor for the Union class.
*/
template <typename... BodyTypes>
Union<BodyTypes...>::~Union()
@@ -266,7 +256,7 @@ Union<BodyTypes...>::~Union()
// Clearing the bodies
bodies_.clear();
- // Logging the destruction of the sphere
+ // Logging the destruction of the union
WALBERLA_LOG_DETAIL( "Destroyed union " << sid_ );
}
//*************************************************************************************************
@@ -321,27 +311,31 @@ template <typename... BodyTypes>
void Union<BodyTypes...>::calcBoundingBox()
{
// Setting the bounding box of an empty union
+ Vec3 gpos = getPosition();
if( bodies_.isEmpty() ) {
aabb_ = math::AABB(
- gpos_[0] - real_t(0.01),
- gpos_[1] - real_t(0.01),
- gpos_[2] - real_t(0.01),
- gpos_[0] + real_t(0.01),
- gpos_[1] + real_t(0.01),
- gpos_[2] + real_t(0.01)
+ gpos[0] - real_t(0.01),
+ gpos[1] - real_t(0.01),
+ gpos[2] - real_t(0.01),
+ gpos[0] + real_t(0.01),
+ gpos[1] + real_t(0.01),
+ gpos[2] + real_t(0.01)
);
}
// Using the bounding box of the first contained bodies as initial bounding box
// and merging it with the bounding boxes of all other bodies
else {
+ bodies_.begin()->calcBoundingBox();
aabb_ = bodies_.begin()->getAABB();
- for( auto b=bodies_.begin()+1; b!=bodies_.end(); ++b )
- aabb_.merge( b->getAABB() );
+ for( auto &b : bodies_ ){
+ b.calcBoundingBox();
+ aabb_.merge( b.getAABB() );
+ }
}
WALBERLA_ASSERT( aabb_.checkInvariant() , "Invalid bounding box detected" );
- WALBERLA_ASSERT( aabb_.contains( gpos_ ), "Invalid bounding box detected" );
+ WALBERLA_ASSERT( aabb_.contains( getPosition() ), "Invalid bounding box detected" );
}
//*************************************************************************************************
@@ -364,64 +358,58 @@ void Union<BodyTypes...>::calcCenterOfMass()
// Don't calculate the center of mass of an empty union
if( bodies_.isEmpty() ) return;
- // Calculating the center of mass of a single body
- if( bodies_.size() == 1 )
- {
- const BodyID body( bodies_.begin().getBodyID() );
- gpos_ = body->getPosition();
- mass_ = body->getMass();
- if( !isFixed() && mass_ > real_t(0) )
- invMass_ = real_t(1) / mass_;
- }
-
- // Calculating the center of mass of a union containing several bodies
- else
+ // Resetting the center of mass
+ Vec3 gpos = Vec3(0,0,0);
+ real_t mass = real_t(0.0);
+ // Calculating the center of mass of a finite union
+ if( finite_ )
{
- // Resetting the center of mass
- gpos_.reset();
-
- // Calculating the center of mass of a finite union
- if( finite_ )
- {
- // Accumulating the mass of all contained rigid bodies
- for( auto b=bodies_.begin(); b!=bodies_.end(); ++b ) {
- WALBERLA_ASSERT( b->isFinite(), "Invalid infinite body in finite union detected" );
- mass_ += b->getMass();
- gpos_ += b->getPosition() * b->getMass();
- }
-
- // Calculating the center of mass for unions with non-zero mass
- if( mass_ > real_t(0) ) {
- if( !isFixed() ) invMass_ = real_t(1) / mass_;
- gpos_ /= mass_;
- }
-
- // Calculating the center of mass for unions with a mass of zero
- else {
- size_t counter( 0 );
- gpos_.reset();
+ // Accumulating the mass of all contained rigid bodies
+ for( auto &b : bodies_ ) {
+ WALBERLA_ASSERT( b.isFinite(), "Invalid infinite body in finite union detected" );
+ mass += b.getMass();
+ gpos += b.getPosition() * b.getMass();
+ }
- for( auto b=bodies_.begin(); b!=bodies_.end(); ++b ) {
- gpos_ += b->getPosition();
- ++counter;
- }
- gpos_ /= counter;
- }
+ mass_ = mass;
+ // Calculating the center of mass for unions with non-zero mass
+ if( mass > real_t(0) ) {
+ if( !isFixed() ) invMass_ = real_t(1) / mass;
+ gpos /= mass;
}
- // Calculating the center of mass of an infinite union
+ // Calculating the center of mass for unions with a mass of zero
else {
size_t counter( 0 );
- for( auto b=bodies_.begin(); b!=bodies_.end(); ++b ) {
- if( b->isFinite() ) continue;
- gpos_ += b->getPosition();
+ for( auto &b : bodies_ ) {
+ gpos += b.getPosition();
++counter;
}
- gpos_ /= counter;
+ gpos /= counter;
+ }
+ }
+
+ // Calculating the center of mass of an infinite union
+ else {
+ size_t counter( 0 );
+
+ for( auto &b : bodies_ ) {
+ if( b.isFinite() ) continue;
+ gpos += b.getPosition();
+ ++counter;
}
+
+ gpos /= counter;
+ }
+
+ // Set new center of mass and adapt all relative positions
+ Vec3 shift = getQuaternion().getInverse().rotate(gpos - getPosition());
+ gpos_ = gpos;
+ for( auto &b : bodies_) {
+ b.setRelPosition(b.getRelPosition() - shift);
}
// Checking the state of the union
@@ -430,8 +418,12 @@ void Union<BodyTypes...>::calcCenterOfMass()
//*************************************************************************************************
+
//*************************************************************************************************
/*!\brief Calculation of the moment of inertia in reference to the body frame of the union.
+ *
+ * Use this function only if the center of gravity is set correctly (e.g. after calling
+ * calcCenterOfMass)
*
* \return void
*/
@@ -455,7 +447,7 @@ void Union<BodyTypes...>::calcInertia()
for( auto b=bodies_.begin(); b!=bodies_.end(); ++b )
{
mass = b->getMass();
- pos = b->getPosition() - gpos_;
+ pos = getRotation() * b->getRelPosition();
I_ += b->getInertia();
I_[0] += mass * ( pos[1]*pos[1] + pos[2]*pos[2] );
@@ -468,9 +460,8 @@ void Union<BodyTypes...>::calcInertia()
I_[7] -= mass * pos[1] * pos[2];
I_[8] += mass * ( pos[0]*pos[0] + pos[1]*pos[1] );
}
-
// Rotating the moment of inertia from the global frame of reference to the body frame of reference
- I_ = R_ * I_ * R_;
+ I_ = getRotation().getTranspose() * I_ * getRotation();
// Calculating the inverse of the body moment of inertia
if( !isFixed() ) Iinv_ = I_.getInverse();
@@ -504,17 +495,7 @@ void Union<BodyTypes...>::calcInertia()
template <typename... BodyTypes>
void Union<BodyTypes...>::setPositionImpl( real_t px, real_t py, real_t pz )
{
- Vec3 gpos = Vec3( px, py, pz );
-
- // Calculating the position change
- const Vec3 dp( gpos - gpos_ );
-
- // Setting the global position
- gpos_ = gpos;
-
- // Updating the contained bodies
- for( auto& b : bodies_ )
- b.update( dp );
+ gpos_ = Vec3( px, py, pz );
Union<BodyTypes...>::calcBoundingBox(); // Setting the axis-aligned bounding box
wake(); // Waking the union from sleep mode
@@ -542,16 +523,12 @@ void Union<BodyTypes...>::setPositionImpl( real_t px, real_t py, real_t pz )
template <typename... BodyTypes>
void Union<BodyTypes...>::setOrientationImpl( real_t r, real_t i, real_t j, real_t k )
{
- const Quat q ( r, i, j, k );
- const Quat dq( q * q_.getInverse() );
+ if ( hasSuperBody() ) return;
+ const Quat q ( r, i, j, k );
q_ = q;
R_ = q_.toRotationMatrix();
- // Updating the contained bodies
- for( auto& b : bodies_ )
- b.update( dq );
-
Union<BodyTypes...>::calcBoundingBox(); // Setting the axis-aligned bounding box
wake(); // Waking the union from sleep mode
signalRotation(); // Signaling the change of orientation to the superordinate body
@@ -560,84 +537,6 @@ void Union<BodyTypes...>::setOrientationImpl( real_t r, real_t i, real_t j, real
-//=================================================================================================
-//
-// SIMULATION FUNCTIONS
-//
-//=================================================================================================
-
-//*************************************************************************************************
-/*!\brief Translation update of a subordinate union.
- *
- * \param dp Change in the global position of the superordinate union.
- * \return void
- *
- * This update function is triggered by the superordinate body in case of a translational
- * movement. This movement involves a change in the global position and the axis-aligned
- * bounding box.
- */
-template <typename... BodyTypes>
-void Union<BodyTypes...>::update( const Vec3& dp )
-{
- // Checking the state of the union
- WALBERLA_ASSERT( checkInvariants(), "Invalid union state detected" );
- WALBERLA_ASSERT( hasSuperBody(), "Invalid superordinate body detected" );
-
- // Updating the global position
- gpos_ += dp;
-
- // Updating the contained bodies
- for( auto& b : bodies_ )
- b.update( dp );
-
- // Setting the axis-aligned bounding box
- Union<BodyTypes...>::calcBoundingBox();
-
- // Checking the state of the union
- WALBERLA_ASSERT( checkInvariants(), "Invalid union state detected" );
-}
-//*************************************************************************************************
-
-
-//*************************************************************************************************
-/*!\brief Rotation update of a subordinate union.
- *
- * \param dq Change in the orientation of the superordinate union.
- * \return void
- *
- * This update function is triggered by the superordinate body in case of a rotational movement.
- * This movement involves a change in the global position, the orientation/rotation and the
- * axis-aligned bounding box of the box.
- */
-template <typename... BodyTypes>
-void Union<BodyTypes...>::update( const Quat& dq )
-{
- // Checking the state of the union
- WALBERLA_ASSERT( checkInvariants(), "Invalid union state detected" );
- WALBERLA_ASSERT( hasSuperBody(), "Invalid superordinate body detected" );
-
- // Calculating the new global position
- gpos_ = sb_->getPosition() + ( sb_->getRotation() * rpos_ );
-
- // Calculating the new orientation and rotation
- q_ = dq * q_;
- R_ = q_.toRotationMatrix();
-
- // Updating the contained bodies
- for( auto& b : bodies_ )
- b.update( dq );
-
- // Setting the axis-aligned bounding box
- Union<BodyTypes...>::calcBoundingBox();
-
- // Checking the state of the union
- WALBERLA_ASSERT( checkInvariants(), "Invalid union state detected" );
-}
-//*************************************************************************************************
-
-
-
-
//=================================================================================================
//
// RIGID BODY MANAGER FUNCTIONS
@@ -660,9 +559,9 @@ void Union<BodyTypes...>::update( const Quat& dq )
* The union takes full responsibility for the newly added body, including the necessary
* memory management. After adding the body to the union, the body is considered part of the
* union. All functions called on the union (as for instance all kinds of set functions,
- * translation or rotation functions) additionally affect the rigid body. However, the body
- * can still individually be positioned, oriented, translated, rotated, or made (in-)visible
- * within the union.\n\n
+ * translation or rotation functions) additionally define the values for the rigid body
+ * according to its relative position and orientation. Do not call any Set-Functions on the
+ * single rigid body.\n
*
*
* \section union_add_infinite Adding infinite rigid bodies
@@ -704,16 +603,23 @@ RigidBody& Union<BodyTypes...>::add( std::unique_ptr<RigidBody>&& body )
if( body->isGlobal() ^ global_ )
throw std::logic_error( "Global flags of body and union do not match" );
+ // Only add non-moving particles to a union.
+ WALBERLA_ASSERT_EQUAL(body->getLinearVel(), Vec3());
+ WALBERLA_ASSERT_EQUAL(body->getAngularVel(), Vec3());
+
// Registering the rigid body
auto& bd = bodies_.add( std::move(body) );
- Vec3 oldCenterOfMass = getPosition();
- Vec3 oldImpulse = getLinearVel() * getMass();
+ Vec3 bdglobalPos = bd.getPosition();
+ Quat bdglobalRot = bd.getQuaternion();
- Vec3 bodyCenterOfMass = bd.getPosition();
- Vec3 bodyImpulse = bd.getLinearVel() * bd.getMass();
+ bd.setSB(this);
+ //having a superbody will forward all getVel/Pos/Rot calls to the superbody from now !
- bd.setSB(this); //having a superbody will forward all getVel calls to superbody!!!
+ bd.setRelPosition(getQuaternion().getInverse().rotate(bdglobalPos - getPosition()));
+ bd.setRelOrientation(bdglobalRot * getQuaternion().getInverse());
+ // Update mass, COG and relative positions
+ calcCenterOfMass();
// Updating the axis-aligned bounding box
if( bodies_.size() == 1 )
@@ -721,24 +627,10 @@ RigidBody& Union<BodyTypes...>::add( std::unique_ptr<RigidBody>&& body )
else
aabb_.merge( bd.getAABB() );
- // Setting the union's total mass and center of mass
- calcCenterOfMass();
-
- // Setting the contained primitives' relative position in reference to the center of mass
- for( auto& b : bodies_ )
- b.calcRelPosition();
// Setting the moment of inertia
calcInertia();
- //moving impuls to union
- setLinearVel( Vec3(0,0,0) );
- setAngularVel( Vec3(0,0,0) );
- addImpulseAtPos( oldImpulse, oldCenterOfMass);
- addImpulseAtPos( bodyImpulse, bodyCenterOfMass);
- bd.setLinearVel( Vec3(0,0,0) );
- bd.setAngularVel( Vec3(0,0,0) );
-
// Signaling the internal modification to the superordinate body
signalModification();
@@ -773,10 +665,6 @@ void Union<BodyTypes...>::translateImpl( real_t dx, real_t dy, real_t dz )
// Changing the global position/reference point
gpos_ += dp;
- // Updating the contained bodies
- for( auto& b : bodies_ )
- b.update( dp );
-
Union<BodyTypes...>::calcBoundingBox(); // Setting the axis-aligned bounding box
wake(); // Waking the union from sleep mode
signalTranslation(); // Signaling the position change to the superordinate body
@@ -813,10 +701,6 @@ void Union<BodyTypes...>::rotateImpl( const Quat& dq )
q_ = dq * q_; // Updating the orientation of the union
R_ = q_.toRotationMatrix(); // Updating the rotation of the union
- // Updating the contained bodies
- for( auto& b : bodies_ )
- b.update( dq );
-
Union<BodyTypes...>::calcBoundingBox(); // Setting the axis-aligned bounding box
wake(); // Waking the union from sleep mode
signalRotation(); // Signaling the change of orientation to the superordinate body
@@ -845,10 +729,6 @@ void Union<BodyTypes...>::rotateAroundOriginImpl( const Quat& dq )
R_ = q_.toRotationMatrix(); // Updating the rotation of the union
gpos_ = dq.rotate( gpos_ ); // Updating the global position of the union
- // Updating the contained bodies
- for( auto& b : bodies_ )
- b.update( dq );
-
Union<BodyTypes...>::calcBoundingBox(); // Setting the axis-aligned bounding box
wake(); // Waking the union from sleep mode
signalTranslation(); // Signaling the position change to the superordinate body
@@ -878,10 +758,6 @@ void Union<BodyTypes...>::rotateAroundPointImpl( const Vec3& point, const Quat &
R_ = q_.toRotationMatrix(); // Updating the rotation of the union
gpos_ = point + dq.rotate( dp ); // Updating the global position of the union
- // Updating the contained bodies
- for( auto& b : bodies_ )
- b.update( dq );
-
Union<BodyTypes...>::calcBoundingBox(); // Setting the axis-aligned bounding box
wake(); // Waking the union from sleep mode
signalTranslation(); // Signaling the position change to the superordinate body
@@ -977,10 +853,6 @@ void Union<BodyTypes...>::handleModification()
// Setting the union's total mass and center of mass
calcCenterOfMass();
- // Setting the contained primitives' relative position in reference to the center of mass
- for( auto& b : bodies_ )
- b.calcRelPosition();
-
// Setting the moment of inertia
calcInertia();
@@ -1007,10 +879,6 @@ void Union<BodyTypes...>::handleTranslation()
// Setting the union's total mass and center of mass
calcCenterOfMass();
- // Setting the contained bodies' relative position in reference to the center of mass
- for( auto& b : bodies_ )
- b.calcRelPosition();
-
// Setting the moment of inertia
calcInertia();
@@ -1074,7 +942,6 @@ void Union<BodyTypes...>::print( std::ostream& os, const char* tab ) const
else os << "*infinite*\n";
os << tab << " Global position = " << gpos_ << "\n"
- << tab << " Relative position = " << rpos_ << "\n"
<< tab << " Linear velocity = " << v_ << "\n"
<< tab << " Angular velocity = " << w_ << "\n";
diff --git a/src/pe/rigidbody/UnionFactory.h b/src/pe/rigidbody/UnionFactory.h
index 62fe5828e..6ad1d934e 100644
--- a/src/pe/rigidbody/UnionFactory.h
+++ b/src/pe/rigidbody/UnionFactory.h
@@ -82,7 +82,7 @@ Union<BodyTypes...>* createUnion( BodyStorage& globalStorage, BlockStorage& bl
const id_t sid = UniqueID<RigidBody>::createGlobal();
WALBERLA_ASSERT_EQUAL(communicating, false);
WALBERLA_ASSERT_EQUAL(infiniteMass, true);
- auto un = std::make_unique<Union<BodyTypes...>>(sid, uid, gpos, Vec3(0,0,0), Quat(), global, false, true);
+ auto un = std::make_unique<Union<BodyTypes...>>(sid, uid, gpos, Quat(), global, false, true);
bd = static_cast<Union<BodyTypes...>*>(&globalStorage.add(std::move(un)));
} else
{
@@ -92,7 +92,7 @@ Union<BodyTypes...>* createUnion( BodyStorage& globalStorage, BlockStorage& bl
const id_t sid( UniqueID<RigidBody>::create() );
BodyStorage& bs = (*block.getData<Storage>(storageID))[0];
- auto un = std::make_unique<Union<BodyTypes...>>(sid, uid, gpos, Vec3(0,0,0), Quat(), global, communicating, infiniteMass);
+ auto un = std::make_unique<Union<BodyTypes...>>(sid, uid, gpos, Quat(), global, communicating, infiniteMass);
un->MPITrait.setOwner(Owner(MPIManager::instance()->rank(), block.getId().getID()));
bd = static_cast<Union<BodyTypes...>*>(&bs.add(std::move(un)));
}
@@ -156,7 +156,7 @@ BoxID createBox( Union<BodyTypes...>* un,
sid = UniqueID<RigidBody>::create();
}
- std::unique_ptr<Box> box = std::make_unique<Box>(sid, uid, gpos, Vec3(0,0,0), Quat(), lengths, material, global, communicating, infiniteMass);
+ std::unique_ptr<Box> box = std::make_unique<Box>(sid, uid, gpos, Quat(), lengths, material, global, communicating, infiniteMass);
box->MPITrait.setOwner( un->MPITrait.getOwner() );
if (box != NULL)
@@ -226,7 +226,7 @@ CapsuleID createCapsule( Union<BodyTypes...>* un,
sid = UniqueID<RigidBody>::create();
}
- std::unique_ptr<Capsule> capsule = std::make_unique<Capsule>(sid, uid, gpos, Vec3(0,0,0), Quat(), radius, length, material, global, communicating, infiniteMass);
+ std::unique_ptr<Capsule> capsule = std::make_unique<Capsule>(sid, uid, gpos, Quat(), radius, length, material, global, communicating, infiniteMass);
capsule->MPITrait.setOwner( un->MPITrait.getOwner() );
if (capsule != NULL)
@@ -284,7 +284,7 @@ SphereID createSphere( Union<BodyTypes...>* un,
sid = UniqueID<RigidBody>::create();
}
- std::unique_ptr<Sphere> sphere = std::make_unique<Sphere>(sid, uid, gpos, Vec3(0,0,0), Quat(), radius, material, global, communicating, infiniteMass);
+ std::unique_ptr<Sphere> sphere = std::make_unique<Sphere>(sid, uid, gpos, Quat(), radius, material, global, communicating, infiniteMass);
sphere->MPITrait.setOwner( un->MPITrait.getOwner() );
if (sphere != NULL)
diff --git a/tests/mesh/MeshMarshalling.cpp b/tests/mesh/MeshMarshalling.cpp
index 15e6bd675..cc62c6a63 100644
--- a/tests/mesh/MeshMarshalling.cpp
+++ b/tests/mesh/MeshMarshalling.cpp
@@ -62,8 +62,8 @@ std::vector<Vector3<real_t>> generateOctahedron( const real_t radius)
void checkMeshEquals(const mesh::TriangleMesh &m1, const mesh::TriangleMesh &m2){
// Very basic checks
WALBERLA_CHECK_FLOAT_EQUAL(mesh::computeVolume(m1), mesh::computeVolume(m2));
- WALBERLA_CHECK_EQUAL(mesh::computeCentroid(m1), mesh::computeCentroid(m2));
- WALBERLA_CHECK_EQUAL(mesh::computeInertiaTensor(m1), mesh::computeInertiaTensor(m2));
+ WALBERLA_CHECK_FLOAT_EQUAL(mesh::toWalberla(mesh::computeCentroid(m1)), mesh::toWalberla(mesh::computeCentroid(m2)));
+ WALBERLA_CHECK_FLOAT_EQUAL(mesh::computeInertiaTensor(m1), mesh::computeInertiaTensor(m2));
}
// Checks two convexPolyhedrons for equality
@@ -71,8 +71,10 @@ void checkConvexPolyhedronEquals(const mesh::pe::ConvexPolyhedron &b1, const mes
WALBERLA_CHECK_FLOAT_EQUAL(b1.getPosition(), b2.getPosition());
WALBERLA_CHECK_FLOAT_EQUAL(b1.getLinearVel(), b2.getLinearVel());
WALBERLA_CHECK_FLOAT_EQUAL(b1.getAngularVel(), b2.getAngularVel());
- WALBERLA_CHECK_EQUAL(b1.getInertia(), b2.getInertia());
+ WALBERLA_CHECK_FLOAT_EQUAL(b1.getInertia(), b2.getInertia());
WALBERLA_CHECK_EQUAL(b1.getMaterial(), b2.getMaterial());
+ WALBERLA_CHECK_FLOAT_EQUAL(b1.getQuaternion(), b2.getQuaternion());
+ WALBERLA_CHECK_FLOAT_EQUAL(b1.getRotation(), b2.getRotation());
// Check equality of the meshes
checkMeshEquals(b1.getMesh(), b2.getMesh());
WALBERLA_CHECK_EQUAL(b1.getID(), b2.getID());
@@ -90,7 +92,7 @@ void testConvexPolyhedron()
MaterialID iron = Material::find("iron");
- mesh::pe::ConvexPolyhedron b1(759846, 1234794, Vec3(real_c(1), real_c(2), real_c(3)), Vec3(0,0,0), Quat(), *octamesh, iron, false, true, false);
+ mesh::pe::ConvexPolyhedron b1(759846, 1234794, Vec3(real_c(1), real_c(2), real_c(3)), Quat(), *octamesh, iron, false, true, false);
b1.setLinearVel(Vec3(real_c(5.2), real_c(6.3), real_c(7.4)));
b1.setAngularVel(Vec3(real_c(1.2), real_c(2.3), real_c(3.4)));
@@ -99,7 +101,7 @@ void testConvexPolyhedron()
mpi::RecvBuffer rb(sb);
auto bPtr = UnmarshalDynamically<BodyTuple>::execute(rb, mesh::pe::ConvexPolyhedron::getStaticTypeID(), math::AABB(Vec3(-100,-100,-100), Vec3(100,100,100)), math::AABB(Vec3(-100,-100,-100), Vec3(100,100,100)));
- mesh::pe::ConvexPolyhedronID b2 = static_cast<mesh::pe::ConvexPolyhedronID>(bPtr.get());
+ auto b2 = dynamic_cast<mesh::pe::ConvexPolyhedronID>(bPtr.get());
checkConvexPolyhedronEquals(b1, *b2);
}
@@ -114,11 +116,12 @@ void testUnion()
qhull.run();
MaterialID iron = Material::find("iron");
-
- UnionT u1(159, 423, Vec3(real_c(1), real_c(2), real_c(3)), Vec3(0,0,0), Quat(), false, false, false);
- u1.add(std::make_unique<mesh::pe::ConvexPolyhedron>(753326, 1267824, Vec3(real_c(2), real_c(2), real_c(3)), Vec3(0,0,0), Quat(), *octamesh, iron, false, true, false));
- u1.add(std::make_unique<mesh::pe::ConvexPolyhedron>(753246, 1233424, Vec3(real_c(-1), real_c(4), real_c(-2)), Vec3(0,0,0), Quat(), *octamesh, iron, false, true, false));
-
+ Quat q(real_t(0.3), real_t(0.1), real_t(0.9), real_t(0.3));
+ Quat w(real_t(0.1), real_t(0.3), real_t(0.9), real_t(0.3));
+ UnionT u1(159, 423, Vec3(real_c(1), real_c(2), real_c(3)), q, false, false, false);
+ u1.add(std::make_unique<mesh::pe::ConvexPolyhedron>(753326, 1267824, Vec3(real_c(2), real_c(2), real_c(3)), w, *octamesh, iron, false, true, false));
+ u1.add(std::make_unique<mesh::pe::ConvexPolyhedron>(753246, 1233424, Vec3(real_c(-1), real_c(4), real_c(-2)), w, *octamesh, iron, false, true, false));
+
u1.setLinearVel(Vec3(real_c(5.2), real_c(6.3), real_c(7.4)));
u1.setAngularVel(Vec3(real_c(1.2), real_c(2.3), real_c(3.4)));
@@ -127,23 +130,24 @@ void testUnion()
mpi::RecvBuffer rb(sb);
auto uPtr = UnmarshalDynamically<BodyTuple>::execute(rb, UnionT::getStaticTypeID(), math::AABB(Vec3(-100,-100,-100), Vec3(100,100,100)), math::AABB(Vec3(-100,-100,-100), Vec3(100,100,100)));
- UnionID u2 = static_cast<UnionID>(uPtr.get());
+ auto u2 = dynamic_cast<UnionID>(uPtr.get());
WALBERLA_CHECK_NOT_NULLPTR( u2 );
-
WALBERLA_CHECK_EQUAL(u1.size(), 2);
WALBERLA_CHECK_EQUAL(u1.size(), u2->size());
- WALBERLA_CHECK_EQUAL(u1.getInertia(), u2->getInertia());
- WALBERLA_CHECK_EQUAL(u1.getPosition(), u2->getPosition());
+ WALBERLA_CHECK_FLOAT_EQUAL(u1.getInertia(), u2->getInertia());
+ WALBERLA_CHECK_FLOAT_EQUAL(u1.getPosition(), u2->getPosition());
WALBERLA_CHECK_FLOAT_EQUAL(u1.getLinearVel(), u2->getLinearVel());
WALBERLA_CHECK_FLOAT_EQUAL(u1.getAngularVel(), u2->getAngularVel());
-
+ WALBERLA_CHECK_FLOAT_EQUAL(u1.getQuaternion(), u2->getQuaternion());
+ WALBERLA_CHECK_FLOAT_EQUAL(u1.getRotation(), u2->getRotation());
+
//getting polyhedrons of first union
- mesh::pe::ConvexPolyhedronID p11 = static_cast<mesh::pe::ConvexPolyhedronID > (u1.begin().getBodyID());
- mesh::pe::ConvexPolyhedronID p21 = static_cast<mesh::pe::ConvexPolyhedronID > ((++(u1.begin())).getBodyID());
+ mesh::pe::ConvexPolyhedronID p11 = dynamic_cast<mesh::pe::ConvexPolyhedronID > (u1.begin().getBodyID());
+ mesh::pe::ConvexPolyhedronID p21 = dynamic_cast<mesh::pe::ConvexPolyhedronID > ((++(u1.begin())).getBodyID());
//getting polyhedrons of second union
- mesh::pe::ConvexPolyhedronID p12 = static_cast<mesh::pe::ConvexPolyhedronID > (u2->begin().getBodyID());
- mesh::pe::ConvexPolyhedronID p22 = static_cast<mesh::pe::ConvexPolyhedronID > ((++(u2->begin())).getBodyID());
+ mesh::pe::ConvexPolyhedronID p12 = dynamic_cast<mesh::pe::ConvexPolyhedronID > (u2->begin().getBodyID());
+ mesh::pe::ConvexPolyhedronID p22 = dynamic_cast<mesh::pe::ConvexPolyhedronID > ((++(u2->begin())).getBodyID());
checkConvexPolyhedronEquals(*p11, *p12);
checkConvexPolyhedronEquals(*p21, *p22);
diff --git a/tests/mesh/MeshPeRaytracing.cpp b/tests/mesh/MeshPeRaytracing.cpp
index ce8f27992..a31dd985a 100644
--- a/tests/mesh/MeshPeRaytracing.cpp
+++ b/tests/mesh/MeshPeRaytracing.cpp
@@ -58,9 +58,9 @@ int CpRayIntersectionTest(const int resolution = 10)
const Vec3 center(1,2,3);
- ConvexPolyhedron cp(0, 0, center, Vec3(0,0,0), Quat(), *mesh, Material::find("iron"), false, true, true);
+ ConvexPolyhedron cp(0, 0, center,Quat(), *mesh, Material::find("iron"), false, true, true);
cp.rotate(real_t(1), real_t(2), real_t(3));
- Box bx(0, 0, center, Vec3(0,0,0), Quat(), Vec3(2,2,2), Material::find("iron"), false, true, true);
+ Box bx(0, 0, center, Quat(), Vec3(2,2,2), Material::find("iron"), false, true, true);
bx.rotate(real_t(1), real_t(2), real_t(3));
real_t dx = real_t(1.0) / static_cast<real_t>(resolution);
diff --git a/tests/pe/BodyFlags.cpp b/tests/pe/BodyFlags.cpp
index fa6db7e1c..cd4f10542 100644
--- a/tests/pe/BodyFlags.cpp
+++ b/tests/pe/BodyFlags.cpp
@@ -70,12 +70,12 @@ int main( int argc, char ** argv )
MaterialID iron = Material::find("iron");
- Sphere refGlobalSphere(1, 0, Vec3(9, 9, 9), Vec3(0,0,0), Quat(), 3, iron, true, false, true);
+ Sphere refGlobalSphere(1, 0, Vec3(9, 9, 9), Quat(), 3, iron, true, false, true);
refGlobalSphere.setLinearVel(Vec3(2,2,2));
SphereID globalSphere = createSphere( *globalStorage, forest->getBlockStorage(), storageID, 0, Vec3(9,9,9), 3, iron, true, false, true);
globalSphere->setLinearVel(Vec3(2,2,2));
- Sphere refFixedSphere(2, 0, Vec3(9,9,14), Vec3(0,0,0), Quat(), 3, iron, false, false, true);
+ Sphere refFixedSphere(2, 0, Vec3(9,9,14), Quat(), 3, iron, false, false, true);
SphereID fixedSphere = createSphere( *globalStorage, forest->getBlockStorage(), storageID, 0, Vec3(9,9,14), 3, iron, false, false, true);
walberla::id_t fixedSphereID = 0;
if (fixedSphere != nullptr) fixedSphereID = fixedSphere->getSystemID();
diff --git a/tests/pe/BodyStorage.cpp b/tests/pe/BodyStorage.cpp
index 730639317..1ebdc53fd 100644
--- a/tests/pe/BodyStorage.cpp
+++ b/tests/pe/BodyStorage.cpp
@@ -31,14 +31,14 @@ using namespace walberla::pe;
class Body1 : public Sphere {
public:
static int refCount;
- Body1(walberla::id_t id, MaterialID matID) : Sphere(id, id, Vec3(0,0,0), Vec3(0,0,0), Quat(), 1, matID, false, true, false) {++refCount;}
+ Body1(walberla::id_t id, MaterialID matID) : Sphere(id, id, Vec3(0,0,0), Quat(), 1, matID, false, true, false) {++refCount;}
~Body1() override {--refCount;}
};
class Body2 : public Sphere {
public:
static int refCount;
- Body2(walberla::id_t id, MaterialID matID) : Sphere(id, id, Vec3(0,0,0), Vec3(0,0,0), Quat(), 1, matID, false, true, false) {++refCount;}
+ Body2(walberla::id_t id, MaterialID matID) : Sphere(id, id, Vec3(0,0,0), Quat(), 1, matID, false, true, false) {++refCount;}
~Body2() override {--refCount;}
};
diff --git a/tests/pe/CMakeLists.txt b/tests/pe/CMakeLists.txt
index ec3aab3f8..ca5346fa4 100644
--- a/tests/pe/CMakeLists.txt
+++ b/tests/pe/CMakeLists.txt
@@ -138,6 +138,11 @@ waLBerla_execute_test( NAME PE_STATICTYPEIDS )
waLBerla_compile_test( NAME PE_UNION FILES Union.cpp DEPENDS core )
waLBerla_execute_test( NAME PE_UNION )
+waLBerla_compile_test( NAME PE_UNIONBEHAVIOR FILES UnionBehavior.cpp DEPENDS core )
+if( WALBERLA_DOUBLE_ACCURACY )
+waLBerla_execute_test( NAME PE_UNIONBEHAVIOR )
+endif()
+
waLBerla_compile_test( NAME PE_RAYTRACING FILES Raytracing.cpp DEPENDS core )
waLBerla_execute_test( NAME PE_RAYTRACING )
diff --git a/tests/pe/CheckVitalParameters.cpp b/tests/pe/CheckVitalParameters.cpp
index 063b1b2b3..77325b726 100644
--- a/tests/pe/CheckVitalParameters.cpp
+++ b/tests/pe/CheckVitalParameters.cpp
@@ -34,7 +34,7 @@ int main( int argc, char ** argv )
WALBERLA_UNUSED(env);
MaterialID iron = Material::find("iron");
- Sphere sphere(walberla::UniqueID<RigidBody>::create(), 0, Vec3(15, 15, 15), Vec3(0,0,0), Quat(), 3, iron, false, true, false );
+ Sphere sphere(walberla::UniqueID<RigidBody>::create(), 0, Vec3(15, 15, 15), Quat(), 3, iron, false, true, false );
sphere.MPITrait.setOwner( Owner(0, walberla::blockforest::BlockID().getID() ) );
checkVitalParameters( &sphere, &sphere );
diff --git a/tests/pe/Collision.cpp b/tests/pe/Collision.cpp
index 76f054583..a69762d7f 100644
--- a/tests/pe/Collision.cpp
+++ b/tests/pe/Collision.cpp
@@ -55,10 +55,10 @@ void checkContact(const Contact& c1, const Contact& c2)
void SphereTest()
{
MaterialID iron = Material::find("iron");
- Sphere sp1(123, 1, Vec3(0,0,0), Vec3(0,0,0), Quat(), 1, iron, false, true, false);
- Sphere sp2(124, 2, Vec3(real_t(1.5),0,0), Vec3(0,0,0), Quat(), 1, iron, false, true, false);
- Sphere sp3(125, 3, Vec3(real_t(3.0),0,0), Vec3(0,0,0), Quat(), 1, iron, false, true, false);
- Sphere sp4(124, 2, Vec3(0,real_t(1.5),0), Vec3(0,0,0), Quat(), 1, iron, false, true, false);
+ Sphere sp1(123, 1, Vec3(0,0,0), Quat(), 1, iron, false, true, false);
+ Sphere sp2(124, 2, Vec3(real_t(1.5),0,0), Quat(), 1, iron, false, true, false);
+ Sphere sp3(125, 3, Vec3(real_t(3.0),0,0), Quat(), 1, iron, false, true, false);
+ Sphere sp4(124, 2, Vec3(0,real_t(1.5),0), Quat(), 1, iron, false, true, false);
Plane pl1(223, 1, Vec3(0,0,0), Vec3(1,1,1).getNormalized(), 0, iron);
CylindricalBoundary cb1(333, 0, Vec3(-100,0,0), 2, iron);
@@ -112,13 +112,13 @@ void SphereTest()
void BoxTest()
{
MaterialID iron = Material::find("iron");
- Box b1(123, 0, Vec3(0,0,0), Vec3(0,0,0), Quat(), Vec3(2,2,2), iron, false, true, false);
- Box b2(124, 0, Vec3(real_t(1.5),0,0), Vec3(0,0,0), Quat(), Vec3(2,2,2), iron, false, true, false);
- Box b3(125, 0, Vec3(real_t(3.0),0,0), Vec3(0,0,0), Quat(), Vec3(2,2,2), iron, false, true, false);
- Box b4(123, 0, Vec3(0,0,0), Vec3(0,0,0), Quat(), Vec3(2,2,2), iron, false, true, false);
+ Box b1(123, 0, Vec3(0,0,0), Quat(), Vec3(2,2,2), iron, false, true, false);
+ Box b2(124, 0, Vec3(real_t(1.5),0,0), Quat(), Vec3(2,2,2), iron, false, true, false);
+ Box b3(125, 0, Vec3(real_t(3.0),0,0), Quat(), Vec3(2,2,2), iron, false, true, false);
+ Box b4(123, 0, Vec3(0,0,0), Quat(), Vec3(2,2,2), iron, false, true, false);
b4.rotate( Vec3(1,1,0), real_t(atan(sqrt(2))) );
- Box b5(123, 0, Vec3(0,0,0), Vec3(0,0,0), Quat(), Vec3(2,2,2), iron, false, true, false);
+ Box b5(123, 0, Vec3(0,0,0), Quat(), Vec3(2,2,2), iron, false, true, false);
b5.rotate( Vec3(0,0,1), real_t(math::M_PI * 0.25) );
b5.rotate( Vec3(1,0,0), real_t(math::M_PI * 0.25) );
@@ -157,7 +157,7 @@ void BoxTest()
b5.setPosition( (Vec3(0,0,1) * real_t(sqrt(3)) + Vec3(0,0,1)) * 1.01);
WALBERLA_CHECK( !collideFunc(&b1, &b5) );
- Sphere s1(126, 0, Vec3(real_t(1.5), real_t(1.5), real_t(1.5)), Vec3(0,0,0), Quat(), 1, iron, false, true, false);
+ Sphere s1(126, 0, Vec3(real_t(1.5), real_t(1.5), real_t(1.5)), Quat(), 1, iron, false, true, false);
WALBERLA_CHECK( collideFunc(&b1, &s1) );
// WALBERLA_LOG_WARNING("contactPoint : " << contactPoint);
// WALBERLA_LOG_WARNING("contactNormal : " << contactNormal);
@@ -167,8 +167,8 @@ void BoxTest()
void CapsuleTest()
{
MaterialID iron = Material::find("iron");
- Capsule c1(100, 100, Vec3(0,0,0), Vec3(0,0,0), Quat(), 1, 2, iron, false, true, false);
- Sphere sp1(123, 123, Vec3(0,0,0), Vec3(0,0,0), Quat(), 1, iron, false, true, false);
+ Capsule c1(100, 100, Vec3(0,0,0), Quat(), 1, 2, iron, false, true, false);
+ Sphere sp1(123, 123, Vec3(0,0,0), Quat(), 1, iron, false, true, false);
std::vector<Contact> contacts;
fcd::AnalyticCollideFunctor< std::vector<Contact> > collideFunc(contacts);
@@ -210,9 +210,9 @@ void CapsuleTest2()
const real_t dynamic_cof ( static_cof ); // Coefficient of dynamic friction. Similar to static friction for low speed friction.
MaterialID material = createMaterial( "granular", real_t( 1.0 ), 0, static_cof, dynamic_cof, real_t( 0.5 ), 1, 1, 0, 0 );
//create obstacle
- Capsule c1(100, 100, Vec3(10,10,0), Vec3(0,0,0), Quat(), 3, 40, material, false, true, false);
+ Capsule c1(100, 100, Vec3(10,10,0), Quat(), 3, 40, material, false, true, false);
c1.rotate( Vec3(0,1,0), math::M_PI * real_t(0.5) );
- Sphere sp1(123, 123, Vec3(real_t(6.5316496854295262864), real_t(10.099999999999999645), real_t(0.46999999991564372914) ), Vec3(0,0,0), Quat(), real_t(0.47), material, false, true, false);
+ Sphere sp1(123, 123, Vec3(real_t(6.5316496854295262864), real_t(10.099999999999999645), real_t(0.46999999991564372914) ), Quat(), real_t(0.47), material, false, true, false);
std::vector<Contact> contacts;
@@ -236,8 +236,8 @@ void CapsuleTest2()
void UnionTest()
{
using UnionT = Union<Sphere>;
- UnionT un1(120, 0, Vec3(0,0,0), Vec3(0,0,0), Quat(), false, true, false);
- UnionT un2(121, 0, Vec3(real_t(1.5),0,0), Vec3(0,0,0), Quat(), false, true, false);
+ UnionT un1(120, 0, Vec3(0,0,0), Quat(), false, true, false);
+ UnionT un2(121, 0, Vec3(real_t(1.5),0,0), Quat(), false, true, false);
auto sp1 = createSphere(&un1, 123, Vec3(0,0,0), 1);
auto sp2 = createSphere(&un2, 124, Vec3(real_t(1.5),0,0), 1);
diff --git a/tests/pe/CollisionTobiasGJK.cpp b/tests/pe/CollisionTobiasGJK.cpp
index 6473222d0..41bdf1f69 100644
--- a/tests/pe/CollisionTobiasGJK.cpp
+++ b/tests/pe/CollisionTobiasGJK.cpp
@@ -160,9 +160,9 @@ void MainTest()
MaterialID iron = Material::find("iron");
// Original SPHERE <-> SPHERE
- Sphere sp1(123, 1, Vec3(0,0,0), Vec3(0,0,0), Quat(), 1, iron, false, true, false);
- Sphere sp2(124, 2, Vec3(real_t(1.5),0,0), Vec3(0,0,0), Quat(), 1, iron, false, true, false);
- Sphere sp3(125, 3, Vec3(real_t(3.0),0,0), Vec3(0,0,0), Quat(), 1, iron, false, true, false);
+ Sphere sp1(123, 1, Vec3(0,0,0), Quat(), 1, iron, false, true, false);
+ Sphere sp2(124, 2, Vec3(real_t(1.5),0,0), Quat(), 1, iron, false, true, false);
+ Sphere sp3(125, 3, Vec3(real_t(3.0),0,0), Quat(), 1, iron, false, true, false);
Vec3 normal;
Vec3 contactPoint;
@@ -179,8 +179,8 @@ void MainTest()
WALBERLA_LOG_INFO("Test 01: BOX <-> SPHERE");
real_t sqr3_inv = real_t(1.0)/std::sqrt(real_t(3.0));
real_t coordinate= real_t(5.0)* sqr3_inv + real_t(5.0); // 5*(1+ (1/sqrt(3)))
- Box box1_1(127, 5, Vec3(0, 0, 0), Vec3(0,0,0), Quat(), Vec3(10, 10, 10), iron, false, true, false);
- Sphere sphere1_2(130, 8, Vec3(coordinate, coordinate, coordinate), Vec3(0,0,0), Quat(), 5, iron, false, true, false);
+ Box box1_1(127, 5, Vec3(0, 0, 0), Quat(), Vec3(10, 10, 10), iron, false, true, false);
+ Sphere sphere1_2(130, 8, Vec3(coordinate, coordinate, coordinate), Quat(), 5, iron, false, true, false);
Vec3 wp1(real_t(5.0), real_t(5.0), real_t(5.0));
Vec3 wpm1(sqr3_inv*real_t(-0.5), sqr3_inv*real_t(-0.5), sqr3_inv*real_t(-0.5));
Vec3 axis1(-sqr3_inv, -sqr3_inv, -sqr3_inv);
@@ -189,7 +189,7 @@ void MainTest()
//Testcase 02 Box LongBox (touching plane)
//Reuse box1_1
WALBERLA_LOG_INFO("Test 02: BOX <-> LONG BOX");
- Box box2_1(131, 9, Vec3(real_t(20.0),0,0), Vec3(0,0,0), Quat(), Vec3(real_t(30.0),1,1), iron, false, true, false);
+ Box box2_1(131, 9, Vec3(real_t(20.0),0,0), Quat(), Vec3(real_t(30.0),1,1), iron, false, true, false);
Vec3 wp2(5, 0, 0);
Vec3 wpm2(real_t(-0.5),0,0);
Vec3 axis2(-1,0,0);
@@ -197,8 +197,8 @@ void MainTest()
//Testcase 03 Sphere Sphere
WALBERLA_LOG_INFO("Test 03: SPHERE <-> SPHERE");
- Sphere sphere3_1(129, 7, Vec3(0,0,0), Vec3(0,0,0), Quat(), 5, iron, false, true, false);
- Sphere sphere3_2(128, 6, Vec3(real_t(10.0),0,0), Vec3(0,0,0), Quat(), 5, iron, false, true, false);
+ Sphere sphere3_1(129, 7, Vec3(0,0,0), Quat(), 5, iron, false, true, false);
+ Sphere sphere3_2(128, 6, Vec3(real_t(10.0),0,0), Quat(), 5, iron, false, true, false);
Vec3 wp3(5, 0, 0);
Vec3 wpm3(real_t(-0.5),0,0);
Vec3 axis3(-1,0,0);
@@ -213,8 +213,8 @@ void MainTest()
//create turned box
real_t sqr2 = std::sqrt(real_t(2.0));
- Box box4_1(132, 10, Vec3(real_t(5.0)*(real_t(1.0)+sqr2), real_t(-5.0), 0), Vec3(0,0,0), q4, Vec3(10, 10, 10), iron, false, true, false);
- Box box4_2(133, 11, Vec3(0, 0, 0), Vec3(0,0,0), Quat(), Vec3(10, 10, 10), iron, false, true, false);
+ Box box4_1(132, 10, Vec3(real_t(5.0)*(real_t(1.0)+sqr2), real_t(-5.0), 0), q4, Vec3(10, 10, 10), iron, false, true, false);
+ Box box4_2(133, 11, Vec3(0, 0, 0), Quat(), Vec3(10, 10, 10), iron, false, true, false);
Vec3 wp4(5, -5, 0);
Vec3 wpm4(real_t(-0.25),real_t(+0.25),0);
Vec3 collision_axis4(-sqr2/real_t(2.0),+sqr2/real_t(2.0),0);
@@ -224,8 +224,8 @@ void MainTest()
//Testcase 05 Cube and Long Box non-centric (touching plane)
WALBERLA_LOG_INFO("Test 05: CUBE <-> LONG BOX (NON_CENTRIC)");
- Box box5_1(133, 12, Vec3(0, 0, 0), Vec3(0,0,0), Quat(), Vec3(10, 10, 10), iron, false, true, false);
- Box box5_2(134, 13, Vec3(real_t(15.0),real_t(5.5), 0), Vec3(0,0,0), Quat(), Vec3(real_t(30.0),1,1), iron, false, true, false);
+ Box box5_1(133, 12, Vec3(0, 0, 0), Quat(), Vec3(10, 10, 10), iron, false, true, false);
+ Box box5_2(134, 13, Vec3(real_t(15.0),real_t(5.5), 0), Quat(), Vec3(real_t(30.0),1,1), iron, false, true, false);
Vec3 wp5(real_t(3.75), 5, 0);
Vec3 wpm5(0, real_t(-0.5), 0);
Vec3 axis5(0, -1, 0);
@@ -243,8 +243,8 @@ void MainTest()
Quat q6(rot_axis6, angle6);
//create turned box with pos = (5*(1+sqrt(3)), 0, 0)
- Box box6_1(136, 14, Vec3(real_t(5.0)*(real_t(1.0)+sqr6_3), 0, 0), Vec3(0,0,0), q6, Vec3(10, 10, 10), iron, false, true, false);
- Box box6_2(136, 15, Vec3(0, 0, 0), Vec3(0,0,0), Quat(), Vec3(10, 10, 10), iron, false, true, false);
+ Box box6_1(136, 14, Vec3(real_t(5.0)*(real_t(1.0)+sqr6_3), 0, 0), q6, Vec3(10, 10, 10), iron, false, true, false);
+ Box box6_2(136, 15, Vec3(0, 0, 0), Quat(), Vec3(10, 10, 10), iron, false, true, false);
Vec3 wp6(5, 0, 0);
Vec3 wpm6(real_t(-0.5), 0, 0);
Vec3 axis6(-1, 0, 0);
@@ -253,8 +253,8 @@ void MainTest()
//Testcase 07:
// BOX <-> SPHERE
WALBERLA_LOG_INFO("Test 07: BOX <-> SPHERE");
- Sphere sphere7_1(137, 16, Vec3(0,0,0), Vec3(0,0,0), Quat(), 5, iron, false, true, false);
- Box box7_2(138, 17, Vec3(0, 0,real_t(7.5)), Vec3(0,0,0), Quat(), Vec3(5, 5, 5), iron, false, true, false);
+ Sphere sphere7_1(137, 16, Vec3(0,0,0), Quat(), 5, iron, false, true, false);
+ Box box7_2(138, 17, Vec3(0, 0,real_t(7.5)), Quat(), Vec3(5, 5, 5), iron, false, true, false);
Vec3 wpm7(0, 0, real_t(-0.5));
Vec3 wp7(0, 0, real_t(5.0));
Vec3 axis7(0, 0, real_t(-1.0));
@@ -264,8 +264,8 @@ void MainTest()
// CAPSULE <-> CAPSULE
WALBERLA_LOG_INFO("Test 08: CAPSULE <-> CAPSULE");
Quat q8(Vec3(0,1,0), walberla::math::M_PI/real_t(2.0)); //creates a y-axis aligned capsule
- Capsule cap8_1(139, 18, Vec3(0,0,0), Vec3(0,0,0), Quat(), real_t(4.0), real_t(10.0), iron, false, true, false);
- Capsule cap8_2(140, 19, Vec3(0,0, real_t(13.0)), Vec3(0,0,0), q8, real_t(4.0), real_t(10.0), iron, false, true, false);
+ Capsule cap8_1(139, 18, Vec3(0,0,0), Quat(), real_t(4.0), real_t(10.0), iron, false, true, false);
+ Capsule cap8_2(140, 19, Vec3(0,0, real_t(13.0)), q8, real_t(4.0), real_t(10.0), iron, false, true, false);
Vec3 wpm8(0, 0, real_t(-0.5));
Vec3 wp8(0, 0, real_t(4.0));
Vec3 axis8(0, 0, real_t(-1.0));
@@ -274,8 +274,8 @@ void MainTest()
//Testcase 09:
// ELLIPSOID <-> ELLIPSOID
WALBERLA_LOG_INFO("Test 09: ELLIPSOID <-> ELLIPSOID");
- Ellipsoid ell9_1(141, 20, Vec3(0,0,0), Vec3(0,0,0), Quat(), Vec3(10,5,5), iron, false, true, false);
- Ellipsoid ell9_2(142, 21, Vec3(15,0,0), Vec3(0,0,0), Quat(), Vec3(5,10,5), iron, false, true, false);
+ Ellipsoid ell9_1(141, 20, Vec3(0,0,0), Quat(), Vec3(10,5,5), iron, false, true, false);
+ Ellipsoid ell9_2(142, 21, Vec3(15,0,0), Quat(), Vec3(5,10,5), iron, false, true, false);
Vec3 wpm9(real_t(-0.5), 0, 0);
Vec3 wp9(real_t(10), 0, 0);
Vec3 axis9(real_t(-1.0), 0, 0);
@@ -292,8 +292,8 @@ void PlaneTest()
fcd::GenericFCD<BodyTuple, fcd::GJKEPACollideFunctor> testFCD;
Plane pl(1, 1, Vec3(0, 1, 0), Vec3(0, 1, 0), real_t(1.0), iron );
- Sphere sphere(2, 2, Vec3(0, real_t(1.9), 0), Vec3(0,0,0), Quat(), 1, iron, false, true, false);
- Sphere sphere2(3, 3, Vec3(0, real_t(0.1), 0), Vec3(0,0,0), Quat(), 1, iron, false, true, false);
+ Sphere sphere(2, 2, Vec3(0, real_t(1.9), 0), Quat(), 1, iron, false, true, false);
+ Sphere sphere2(3, 3, Vec3(0, real_t(0.1), 0), Quat(), 1, iron, false, true, false);
PossibleContacts pcs;
@@ -353,17 +353,17 @@ void UnionTest(){
fcd::GenericFCD<BodyTuple, fcd::GJKEPACollideFunctor> testFCD;
//A recursive union of three spheres is dropped on a box.
- Box box(179, 179, Vec3(0,0,0), Vec3(0,0,0), Quat(), Vec3(real_t(10),real_t(2), real_t(10)), iron, false, true, false);
+ Box box(179, 179, Vec3(0,0,0), Quat(), Vec3(real_t(10),real_t(2), real_t(10)), iron, false, true, false);
using UnionT = Union<Sphere>;
- auto unsub = std::make_unique<UnionT>(192, 192, Vec3(0,real_t(3.8),0), Vec3(0,0,0), Quat(), false, true, false);
+ auto unsub = std::make_unique<UnionT>(192, 192, Vec3(0,real_t(3.8),0), Quat(), false, true, false);
auto sp1 = createSphere(unsub.get(), 180, Vec3(-3,real_t(3.8),0), real_t(3.0));
auto sp2 = createSphere(unsub.get(), 181, Vec3(3,real_t(3.8),0), real_t(3.0));
//Create another union, and add sub union
- Union<Sphere, Union<Sphere>> un(193, 193, Vec3(0, 0, 0), Vec3(0,0,0), Quat(), false, true, false);
+ Union<Sphere, Union<Sphere>> un(193, 193, Vec3(0,0,0), Quat(), false, true, false);
createSphere(&un, 182, Vec3(0,real_t(6),0), real_t(3.0));
un.add(std::move(unsub));
diff --git a/tests/pe/Marshalling.cpp b/tests/pe/Marshalling.cpp
index fc18f19eb..52702aa49 100644
--- a/tests/pe/Marshalling.cpp
+++ b/tests/pe/Marshalling.cpp
@@ -50,7 +50,7 @@ void testBox()
MaterialID iron = Material::find("iron");
- Box b1(759846, 1234794, Vec3(real_c(1), real_c(2), real_c(3)), Vec3(0,0,0), Quat(), Vec3(1,2,3), iron, false, true, false);
+ Box b1(759846, 1234794, Vec3(real_c(1), real_c(2), real_c(3)), Quat(), Vec3(1,2,3), iron, false, true, false);
b1.setLinearVel(Vec3(real_c(5.2), real_c(6.3), real_c(7.4)));
b1.setAngularVel(Vec3(real_c(1.2), real_c(2.3), real_c(3.4)));
@@ -75,7 +75,7 @@ void testCapsule()
MaterialID iron = Material::find("iron");
- Capsule c1(759846, 1234794, Vec3(real_c(1), real_c(2), real_c(3)), Vec3(0,0,0), Quat(), 5, 7, iron, false, false, false);
+ Capsule c1(759846, 1234794, Vec3(real_c(1), real_c(2), real_c(3)), Quat(), 5, 7, iron, false, false, false);
c1.setLinearVel(Vec3(real_c(5.2), real_c(6.3), real_c(7.4)));
c1.setAngularVel(Vec3(real_c(1.2), real_c(2.3), real_c(3.4)));
@@ -101,7 +101,7 @@ void testSphere()
MaterialID iron = Material::find("iron");
- Sphere s1(759846, 1234794, Vec3(real_c(1), real_c(2), real_c(3)), Vec3(0,0,0), Quat(), 5, iron, false, false, false);
+ Sphere s1(759846, 1234794, Vec3(real_c(1), real_c(2), real_c(3)), Quat(), 5, iron, false, false, false);
s1.setLinearVel(Vec3(real_c(5.2), real_c(6.3), real_c(7.4)));
s1.setAngularVel(Vec3(real_c(1.2), real_c(2.3), real_c(3.4)));
@@ -126,7 +126,7 @@ void testSquirmer()
MaterialID iron = Material::find("iron");
- Squirmer s1(759846, 1234794, Vec3(real_c(1), real_c(2), real_c(3)), Vec3(0,0,0), Quat(), real_c(5), real_c(0.1), real_c(4.93), iron, false, false, false);
+ Squirmer s1(759846, 1234794, Vec3(real_c(1), real_c(2), real_c(3)), Quat(), real_c(5), real_c(0.1), real_c(4.93), iron, false, false, false);
mpi::SendBuffer sb;
WALBERLA_ASSERT_UNEQUAL(Sphere::getStaticTypeID(), Squirmer::getStaticTypeID(), "Squirmer did not get its own type ID");
@@ -147,7 +147,7 @@ void testEllipsoid()
MaterialID iron = Material::find("iron");
- Ellipsoid e1(759847, 1234795, Vec3(real_c(1), real_c(2), real_c(3)), Vec3(0,0,0), Quat(), Vec3(3,1,5), iron, false, false, false);
+ Ellipsoid e1(759847, 1234795, Vec3(real_c(1), real_c(2), real_c(3)), Quat(), Vec3(3,1,5), iron, false, false, false);
e1.setLinearVel(Vec3(real_c(5.2), real_c(6.3), real_c(7.4)));
e1.setAngularVel(Vec3(real_c(1.2), real_c(2.3), real_c(3.4)));
@@ -169,7 +169,7 @@ void testEllipsoid()
void testUnion()
{
WALBERLA_LOG_INFO_ON_ROOT("*** testUnion ***");
- UnionT u1(159, 423, Vec3(real_c(1), real_c(2), real_c(3)), Vec3(0,0,0), Quat(), false, false, false);
+ UnionT u1(159, 423, Vec3(real_c(1), real_c(2), real_c(3)), Quat(), false, false, false);
SphereID s11 = createSphere(&u1, 1234794, Vec3(real_c(1), real_c(2), real_c(3)), 2);
SphereID s21 = createSphere(&u1, 4567789, Vec3(real_c(3), real_c(2), real_c(3)), real_c(1.5));
WALBERLA_CHECK_NOT_NULLPTR( s11 );
@@ -185,6 +185,8 @@ void testUnion()
WALBERLA_CHECK_EQUAL(u1.size(), 2);
WALBERLA_CHECK_EQUAL(u1.size(), u2->size());
+ // More exhaustive tests (with inertia, rotation, etc.)
+ // can be found in tests/mesh/MeshMarshalling.cpp
//getting spheres of second union
SphereID s12 = static_cast<SphereID> (u2->begin().getBodyID());
diff --git a/tests/pe/Raytracing.cpp b/tests/pe/Raytracing.cpp
index 25e925b8a..15dd6ae60 100644
--- a/tests/pe/Raytracing.cpp
+++ b/tests/pe/Raytracing.cpp
@@ -43,7 +43,7 @@ typedef std::tuple<Box, Plane, Sphere, Capsule, Ellipsoid> BodyTuple ;
void SphereIntersectsTest()
{
MaterialID iron = Material::find("iron");
- Sphere sp1(123, 1, Vec3(3,3,3), Vec3(0,0,0), Quat(), 2, iron, false, true, false);
+ Sphere sp1(123, 1, Vec3(3,3,3), Quat(), 2, iron, false, true, false);
real_t t;
Vec3 n;
@@ -62,11 +62,11 @@ void SphereIntersectsTest()
WALBERLA_CHECK(intersects(&sp1, ray2, t, n));
// sphere behind ray origin
- Sphere sp2(123, 1, Vec3(3,-8,3), Vec3(0,0,0), Quat(), 2, iron, false, true, false);
+ Sphere sp2(123, 1, Vec3(3,-8,3), Quat(), 2, iron, false, true, false);
WALBERLA_CHECK(!intersects(&sp2, ray1, t, n));
// sphere around ray origin
- Sphere sp3(123, 1, Vec3(3,-5,3), Vec3(0,0,0), Quat(), 2, iron, false, true, false);
+ Sphere sp3(123, 1, Vec3(3,-5,3), Quat(), 2, iron, false, true, false);
WALBERLA_CHECK(intersects(&sp3, ray1, t, n));
WALBERLA_CHECK_FLOAT_EQUAL(t, real_t(2));
}
@@ -128,15 +128,15 @@ void BoxIntersectsTest() {
real_t t;
Vec3 n;
- Box box1(127, 5, Vec3(0, -15, 0), Vec3(0, 0, 0), Quat(), Vec3(10, 10, 10), iron, false, true, false);
+ Box box1(127, 5, Vec3(0, -15, 0), Quat(), Vec3(10, 10, 10), iron, false, true, false);
Ray ray1(Vec3(3,-5,3), Vec3(0,1,0));
WALBERLA_CHECK(!intersects(&box1, ray1, t, n));
- Box box2(128, 5, Vec3(0, -2, 0), Vec3(0, 0, 0), Quat(), Vec3(10, 10, 10), iron, false, true, false);
+ Box box2(128, 5, Vec3(0, -2, 0), Quat(), Vec3(10, 10, 10), iron, false, true, false);
WALBERLA_CHECK(intersects(&box2, ray1, t, n));
WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(t, real_t(8), real_t(1e-7));
- Box box3(128, 5, Vec3(0, 5, 0), Vec3(0, 0, 0), Quat(), Vec3(10, 10, 10), iron, false, true, false);
+ Box box3(128, 5, Vec3(0, 5, 0), Quat(), Vec3(10, 10, 10), iron, false, true, false);
WALBERLA_CHECK(intersects(&box3, ray1, t, n));
WALBERLA_CHECK_FLOAT_EQUAL(t, real_t(5));
@@ -160,14 +160,14 @@ void BoxIntersectsTest() {
WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(t, real_t(7.0710), real_t(1e-4));
Ray ray3(Vec3(3,-5,3), Vec3(2, real_t(-1.5), real_t(0.5)).getNormalized());
- Box box4(128, 5, Vec3(0, 8, 0), Vec3(0, 0, 0), Quat(), Vec3(10, 10, 10), iron, false, true, false);
+ Box box4(128, 5, Vec3(0, 8, 0), Quat(), Vec3(10, 10, 10), iron, false, true, false);
WALBERLA_CHECK(!intersects(&box4, ray3, t, n));
Ray ray4(Vec3(3,-5,3), Vec3(-2, 3, real_t(0.5)).getNormalized());
WALBERLA_CHECK(intersects(&box4, ray4, t, n));
WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(t, real_t(9.7068), real_t(1e-4));
- Box box5(128, 5, Vec3(4, 0, 0), Vec3(0, 0, 0), Quat(), Vec3(4, 4, 4), iron, false, true, false);
+ Box box5(128, 5, Vec3(4, 0, 0), Quat(), Vec3(4, 4, 4), iron, false, true, false);
box5.rotate(0,0,math::M_PI/4);
Ray ray5(Vec3(0,1.5,0), Vec3(1,0,0));
WALBERLA_CHECK(intersects(&box5, ray5, t, n));
@@ -202,7 +202,7 @@ void CapsuleIntersectsTest() {
real_t t;
Vec3 n;
- Capsule cp1(0, 0, Vec3(2,3,3), Vec3(0,0,0), Quat(), real_t(2), real_t(2), iron, false, true, false);
+ Capsule cp1(0, 0, Vec3(2,3,3), Quat(), real_t(2), real_t(2), iron, false, true, false);
// ray through the center
Ray ray1(Vec3(3,-5,3), Vec3(0,1,0));
@@ -229,7 +229,7 @@ void EllipsoidTest() {
real_t t;
Vec3 n;
- Ellipsoid el1(0,0, Vec3(2,3,3), Vec3(0,0,0), Quat(), Vec3(2,3,1), iron, false, true, false);
+ Ellipsoid el1(0,0, Vec3(2,3,3), Quat(), Vec3(2,3,1), iron, false, true, false);
Ray ray1(Vec3(-2,3,3), Vec3(1,0,0).getNormalized());
WALBERLA_CHECK(intersects(&el1, ray1, t, n));
diff --git a/tests/pe/RigidBody.cpp b/tests/pe/RigidBody.cpp
index f72d9ce0b..083c6ab43 100644
--- a/tests/pe/RigidBody.cpp
+++ b/tests/pe/RigidBody.cpp
@@ -82,10 +82,10 @@ void move( BodyStorage& storage, real_t dt )
void checkRotationFunctions()
{
MaterialID iron = Material::find("iron");
- auto sp1 = std::make_shared<Sphere>( 0, 0, Vec3(0,0,0), Vec3(0,0,0), Quat(), real_t(1), iron, false, true, false );
- auto sp2 = std::make_shared<Sphere>( 0, 0, Vec3(0,0,0), Vec3(0,0,0), Quat(), real_t(1), iron, false, true, false );
- auto sp3 = std::make_shared<Sphere>( 0, 0, Vec3(0,0,0), Vec3(0,0,0), Quat(), real_t(1), iron, false, true, false );
- auto sp4 = std::make_shared<Sphere>( 0, 0, Vec3(0,0,0), Vec3(0,0,0), Quat(), real_t(1), iron, false, true, false );
+ auto sp1 = std::make_shared<Sphere>( 0, 0, Vec3(0,0,0), Quat(), real_t(1), iron, false, true, false );
+ auto sp2 = std::make_shared<Sphere>( 0, 0, Vec3(0,0,0), Quat(), real_t(1), iron, false, true, false );
+ auto sp3 = std::make_shared<Sphere>( 0, 0, Vec3(0,0,0), Quat(), real_t(1), iron, false, true, false );
+ auto sp4 = std::make_shared<Sphere>( 0, 0, Vec3(0,0,0), Quat(), real_t(1), iron, false, true, false );
sp1->rotate( 1, 0, 0, math::M_PI * real_t(0.5));
sp1->rotate( 0, 1, 0, math::M_PI * real_t(0.5));
@@ -121,7 +121,7 @@ void checkRotationFunctions()
void checkPointFunctions()
{
MaterialID iron = Material::find("iron");
- auto sp1 = std::make_shared<Sphere>( 0, 0, Vec3(10,10,10), Vec3(0,0,0), Quat(), real_t(1), iron, false, true, false );
+ auto sp1 = std::make_shared<Sphere>( 0, 0, Vec3(10,10,10), Quat(), real_t(1), iron, false, true, false );
WALBERLA_CHECK( sp1->containsPoint( 10, 10, 10 ) );
WALBERLA_CHECK( sp1->containsPoint( real_c(10.9), 10, 10 ) );
@@ -148,7 +148,7 @@ int main( int argc, char** argv )
MaterialID iron = Material::find("iron");
BodyStorage storage;
- SpherePtr spPtr = std::make_unique<Sphere>(0, 0, Vec3(0,0,0), Vec3(0,0,0), Quat(), real_t(1), iron, false, true, false);
+ SpherePtr spPtr = std::make_unique<Sphere>(0, 0, Vec3(0,0,0), Quat(), real_t(1), iron, false, true, false);
SphereID sphere = static_cast<SphereID>(&storage.add(std::move(spPtr)));
Vec3 x0 = Vec3(-2,2,0);
diff --git a/tests/pe/SetBodyTypeIDs.cpp b/tests/pe/SetBodyTypeIDs.cpp
index 21e2fcd5c..f6e1482a0 100644
--- a/tests/pe/SetBodyTypeIDs.cpp
+++ b/tests/pe/SetBodyTypeIDs.cpp
@@ -132,8 +132,8 @@ int main( int argc, char** argv )
WALBERLA_CHECK_UNEQUAL(Plane::getStaticTypeID(), Capsule::getStaticTypeID());
SingleTypeFunctor singleFunc;
- Box bx (0, 0, Vec3(0), Vec3(0), Quat(), Vec3(1), Material::find("iron"), false, false, false);
- Capsule cap(0, 0, Vec3(0), Vec3(0), Quat(), 1, 1, Material::find("iron"), false, false, false);
+ Box bx (0, 0, Vec3(0), Quat(), Vec3(1), Material::find("iron"), false, false, false);
+ Capsule cap(0, 0, Vec3(0), Quat(), 1, 1, Material::find("iron"), false, false, false);
SingleCast<BodyTuple2, SingleTypeFunctor, void>::execute(Box::getStaticTypeID(), singleFunc);
SingleCast<BodyTuple2, SingleTypeFunctor, void>::execute(Capsule::getStaticTypeID(), singleFunc);
diff --git a/tests/pe/SimpleCCD.cpp b/tests/pe/SimpleCCD.cpp
index 3d0d5bc7e..63e9dc239 100644
--- a/tests/pe/SimpleCCD.cpp
+++ b/tests/pe/SimpleCCD.cpp
@@ -60,7 +60,7 @@ int main( int argc, char** argv )
math::seedRandomGenerator(1337);
for (uint_t i = 0; i < 100; ++i)
- storage[0].add( std::make_unique<Sphere>(UniqueID<Sphere>::createGlobal(), 0, Vec3( math::realRandom(real_c(0), real_c(10)), math::realRandom(real_c(0), real_c(10)), math::realRandom(real_c(0), real_c(10))), Vec3(0,0,0), Quat(), real_t(1), iron, false, false, false) );
+ storage[0].add( std::make_unique<Sphere>(UniqueID<Sphere>::createGlobal(), 0, Vec3( math::realRandom(real_c(0), real_c(10)), math::realRandom(real_c(0), real_c(10)), math::realRandom(real_c(0), real_c(10))), Quat(), real_t(1), iron, false, false, false) );
sccd.generatePossibleContacts();
@@ -84,14 +84,14 @@ int main( int argc, char** argv )
bs.clear();
- bs.add( std::make_unique<Sphere>(UniqueID<Sphere>::createGlobal(), 0, Vec3( math::realRandom(real_c(0), real_c(10)), math::realRandom(real_c(0), real_c(10)), math::realRandom(real_c(0), real_c(10))), Vec3(0,0,0), Quat(), real_t(1), iron, false, false, false) );
+ bs.add( std::make_unique<Sphere>(UniqueID<Sphere>::createGlobal(), 0, Vec3( math::realRandom(real_c(0), real_c(10)), math::realRandom(real_c(0), real_c(10)), math::realRandom(real_c(0), real_c(10))), Quat(), real_t(1), iron, false, false, false) );
WcTimingPool pool;
for (int runs = 0; runs < 10; ++runs)
{
auto oldSize = bs.size();
for (uint_t i = 0; i < oldSize; ++i)
- bs.add( std::make_unique<Sphere>(UniqueID<Sphere>::createGlobal(), 0, Vec3( math::realRandom(real_c(0), real_c(10)), math::realRandom(real_c(0), real_c(10)), math::realRandom(real_c(0), real_c(10))), Vec3(0,0,0), Quat(), real_t(0.5), iron, false, false, false) );
+ bs.add( std::make_unique<Sphere>(UniqueID<Sphere>::createGlobal(), 0, Vec3( math::realRandom(real_c(0), real_c(10)), math::realRandom(real_c(0), real_c(10)), math::realRandom(real_c(0), real_c(10))), Quat(), real_t(0.5), iron, false, false, false) );
pool["SCCD"].start();
sccd.generatePossibleContacts();
pool["SCCD"].end();
diff --git a/tests/pe/Synchronization.cpp b/tests/pe/Synchronization.cpp
index 3693e9ab8..67dd06832 100644
--- a/tests/pe/Synchronization.cpp
+++ b/tests/pe/Synchronization.cpp
@@ -116,7 +116,7 @@ int main( int argc, char ** argv )
MaterialID iron = Material::find("iron");
walberla::id_t sid = 123;
- Sphere refSphere(1, 0, Vec3(15, 15, 15), Vec3(0,0,0), Quat(), 3, iron, false, true, false);
+ Sphere refSphere(1, 0, Vec3(15, 15, 15), Quat(), 3, iron, false, true, false);
refSphere.setLinearVel(4, 5, 6);
refSphere.setAngularVel( 1, 2, 3);
Vec3 gpos = Vec3(15, 15, 15);
diff --git a/tests/pe/SynchronizationLargeBody.cpp b/tests/pe/SynchronizationLargeBody.cpp
index 28fdc39d1..94e205db7 100644
--- a/tests/pe/SynchronizationLargeBody.cpp
+++ b/tests/pe/SynchronizationLargeBody.cpp
@@ -160,7 +160,7 @@ int main( int argc, char ** argv )
walberla::id_t sid = 123;
Vec3 gpos = Vec3(3.5, 3.5, 3.5);
const real_t r = real_c(1.6);
- Sphere refSphere(1, 0, gpos, Vec3(0,0,0), Quat(), r, iron, false, true, false);
+ Sphere refSphere(1, 0, gpos, Quat(), r, iron, false, true, false);
refSphere.setLinearVel(4, 5, 6);
refSphere.setAngularVel( 1, 2, 3);
diff --git a/tests/pe/Union.cpp b/tests/pe/Union.cpp
index 49fd6e66b..de959af08 100644
--- a/tests/pe/Union.cpp
+++ b/tests/pe/Union.cpp
@@ -13,7 +13,7 @@
// You should have received a copy of the GNU General Public License along
// with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
//
-//! \file ParallelEquivalence.cpp
+//! \file Union.cpp
//! \author Sebastian Eibl <sebastian.eibl@fau.de>
//
//======================================================================================================================
@@ -49,6 +49,7 @@ typedef std::tuple<Sphere, Plane, UnionType> BodyTuple ;
void SnowManFallingOnPlane()
{
+ WALBERLA_LOG_INFO("- Performing snowman on plane test");
shared_ptr<BodyStorage> globalBodyStorage = make_shared<BodyStorage>();
// create blocks
@@ -99,25 +100,178 @@ void SnowManFallingOnPlane()
//WALBERLA_LOG_DEVEL(un);
}
-void ImpulsCarryover()
+void UnionConstruction()
{
+
MaterialID iron = Material::find("iron");
- auto un = std::make_unique<UnionType>(12, 0, Vec3(0,0,0), Vec3(0,0,0), Quat(), false, true, false);
- auto sp1 = std::make_unique<Sphere>( 10, 0, Vec3( 1,0,0), Vec3(0,0,0), Quat(), real_t(1), iron, false, true, false );
- auto sp2 = std::make_unique<Sphere>( 11, 0, Vec3(-1,0,0), Vec3(0,0,0), Quat(), real_t(1), iron, false, true, false );
+ // Construct a union of two spheres in a non-rotated state and check all its parameters
+ auto un = std::make_unique<UnionType>(12, 0, Vec3(0,0,0), Quat(), false, true, false);
+ auto sp1 = std::make_unique<Sphere>( 10, 0, Vec3( 1,0,0), Quat(real_t(0.5), real_t(0.5),real_t(0.5),real_t(0.5)), real_t(1), iron, false, true, false );
+ auto sp2 = std::make_unique<Sphere>( 11, 0, Vec3(-1,0,0), Quat(), real_t(1), iron, false, true, false );
+
+ WALBERLA_LOG_INFO("- Adding first body.");
+ // Add first body
+ un->add( std::move(sp1) );
+ const RigidBody& subb1 = *un->begin();
+
+ // Check relative and global position and rotation
+ WALBERLA_CHECK_FLOAT_EQUAL(un->getPosition(),Vec3(1,0,0));
+ WALBERLA_CHECK_FLOAT_EQUAL(un->getRelPosition(), Vec3(0,0,0));
+
+ WALBERLA_CHECK_FLOAT_EQUAL(subb1.getRelPosition(), Vec3(0, 0,0));
+ WALBERLA_CHECK_FLOAT_EQUAL(subb1.getPosition(), Vec3(1,0,0));
+ // Global quaternion
+ WALBERLA_CHECK_FLOAT_EQUAL(subb1.getQuaternion(), Quat(real_t(0.5), real_t(0.5),real_t(0.5),real_t(0.5)));
+ WALBERLA_CHECK_FLOAT_EQUAL(subb1.getRelQuaternion(), Quat(real_t(0.5), real_t(0.5),real_t(0.5),real_t(0.5)));
+
+ // Check mass volume and inertia
+ WALBERLA_CHECK_FLOAT_EQUAL(un->getVolume(), (real_t(4./3.)*real_t(math::M_PI)));
+ WALBERLA_CHECK_FLOAT_EQUAL(un->getMass(), un->getVolume()*Material::getDensity(iron));
+ real_t scalar_inertia = real_t(0.4)*un->getMass(); // for sphere: I = 2/5*m*r*r
+ WALBERLA_CHECK_EQUAL(un->getInertia(), Mat3(scalar_inertia,0,0,0,scalar_inertia,0,0,0,scalar_inertia));
+
+ WALBERLA_LOG_INFO("- Adding second body.");
+ un->add( std::move(sp2) );
+ const RigidBody & subb2 = *(un->begin()+1);
+ // Check relative and global position
+ WALBERLA_CHECK_FLOAT_EQUAL(un->getPosition(),Vec3(0,0,0));
+ WALBERLA_CHECK_FLOAT_EQUAL(un->getRelPosition(), Vec3(0,0,0));
+ WALBERLA_CHECK_FLOAT_EQUAL(subb1.getRelPosition(), Vec3(1,0,0));
+ WALBERLA_CHECK_FLOAT_EQUAL(subb1.getPosition(), Vec3(1,0,0));
+ WALBERLA_CHECK_FLOAT_EQUAL(subb2.getRelPosition(), Vec3(-1,0,0));
+ WALBERLA_CHECK_FLOAT_EQUAL(subb2.getPosition(), Vec3(-1,0,0));
+
+ // Check mass volume and inertia
+ WALBERLA_CHECK_FLOAT_EQUAL(un->getVolume(), (real_t(8./3.)*real_t(math::M_PI)));
+ WALBERLA_CHECK_FLOAT_EQUAL(un->getMass(), un->getVolume()*Material::getDensity(iron));
+ // Mass of one sphere
+ real_t masssphere = real_t(4./3.)*real_t(math::M_PI)*Material::getDensity(iron);
+ Mat3 bodyinertia(real_t(2.0)*scalar_inertia, 0, 0, 0, real_t(2.0)*(scalar_inertia + masssphere),0, 0, 0, real_t(2.0)*(scalar_inertia + masssphere));
+ WALBERLA_CHECK_FLOAT_EQUAL(un->getInertia(), bodyinertia);
+
+ WALBERLA_CHECK_FLOAT_EQUAL( un->getLinearVel(), Vec3(0,0,0));
+ WALBERLA_CHECK_FLOAT_EQUAL( un->getAngularVel(), Vec3(0,0,0));
+ WALBERLA_CHECK_FLOAT_EQUAL( un->getAngularVel(), Vec3(0,0,0));
+
+ WALBERLA_LOG_INFO("- Performing Rotation.");
+ //Check values for rotated union
+ Quat rotz30(Vec3(0,0,1), real_t(math::M_PI/6.0)); // rotate by 30 deg via z axis
+ real_t sin30 = real_t(0.5);
+ real_t cos30 = real_t(sqrt(3.0)/2.0);
+ un->setOrientation(rotz30);
+ WALBERLA_CHECK_FLOAT_EQUAL(un->getQuaternion(), rotz30);
+ WALBERLA_CHECK_FLOAT_EQUAL(subb1.getRelPosition(), Vec3(1,0,0));
+ WALBERLA_CHECK_FLOAT_EQUAL(subb1.getPosition(), Vec3(cos30,sin30,0));
+ WALBERLA_CHECK_FLOAT_EQUAL(subb2.getRelPosition(), Vec3(-1,0,0));
+ WALBERLA_CHECK_FLOAT_EQUAL(subb2.getPosition(), Vec3(-cos30,-sin30,0));
+ WALBERLA_CHECK_FLOAT_EQUAL(subb1.getRelQuaternion(), Quat(real_t(0.5), real_t(0.5),real_t(0.5),real_t(0.5)));
+ WALBERLA_CHECK_FLOAT_EQUAL(subb1.getQuaternion(), Quat(real_t(0.5), real_t(0.5),real_t(0.5),real_t(0.5))*rotz30 );
+ WALBERLA_CHECK_FLOAT_EQUAL(un->getInertia(), rotz30.toRotationMatrix()*bodyinertia*rotz30.toRotationMatrix().getTranspose());
+
+ WALBERLA_LOG_INFO("- Applying velocities.");
+ //Apply a linear velocity to the union
+ un->setLinearVel(Vec3(1,0,0));
+ WALBERLA_CHECK_FLOAT_EQUAL(subb1.getLinearVel(), Vec3(1,0,0));
+ WALBERLA_CHECK_FLOAT_EQUAL(subb2.getLinearVel(), Vec3(1,0,0));
+
+ un->setAngularVel(Vec3(0,0,1));
+ WALBERLA_CHECK_FLOAT_EQUAL(subb1.getAngularVel(), Vec3(0,0,1));
+ WALBERLA_CHECK_FLOAT_EQUAL(subb2.getAngularVel(), Vec3(0,0,1));
+ WALBERLA_CHECK_FLOAT_EQUAL(subb1.getLinearVel(), Vec3(real_t(1.0-sin30),cos30,0));
+ WALBERLA_CHECK_FLOAT_EQUAL(subb2.getLinearVel(), Vec3(real_t(1.0+sin30),-cos30,0));
+
+ WALBERLA_LOG_INFO("- Constructing rotated union.");
+ //Part 2: Construct exactly the same union, but now in a rotated state.
+ auto un2 = std::make_unique<UnionType>(12, 0, Vec3(0,0,0), rotz30, false, true, false);
+ auto sp12 = std::make_unique<Sphere>( 10, 0, Vec3( 1,0,0), Quat(real_t(0.5), real_t(0.5),real_t(0.5),real_t(0.5)), real_t(1), iron, false, true, false );
+ auto sp22 = std::make_unique<Sphere>( 11, 0, Vec3(-1,0,0), Quat(real_t(0.3), real_t(0.9),real_t(0.1),real_t(0.3)), real_t(1), iron, false, true, false );
+
+ un2->add( std::move(sp12) );
+ un2->add( std::move(sp22) );
+ const RigidBody & subb12 = *un2->begin();
+ const RigidBody & subb22 = *(un2->begin()+1);
+ // Check relative and global position and rotation
+ WALBERLA_CHECK_FLOAT_EQUAL(un2->getPosition(),Vec3(0,0,0));
+ WALBERLA_CHECK_FLOAT_EQUAL(un2->getRelPosition(), Vec3(0,0,0));
+
+ WALBERLA_CHECK_FLOAT_EQUAL(subb12.getRelPosition(), Vec3(cos30,-sin30, 0));
+ WALBERLA_CHECK_FLOAT_EQUAL(subb12.getPosition(), Vec3(1,0,0));
+
+ WALBERLA_CHECK_FLOAT_EQUAL(subb22.getRelPosition(), Vec3(-cos30,sin30, 0));
+ WALBERLA_CHECK_FLOAT_EQUAL(subb22.getPosition(), Vec3(-1,0,0));
- sp1->setLinearVel(Vec3(0,real_c(+1),0));
- sp2->setLinearVel(Vec3(0,real_c(-1),0));
+ // Global quaternion and local
+ WALBERLA_CHECK_FLOAT_EQUAL(subb12.getQuaternion(), Quat(real_t(0.5), real_t(0.5),real_t(0.5),real_t(0.5)));
+ WALBERLA_CHECK_FLOAT_EQUAL(subb12.getRelQuaternion(), subb12.getQuaternion()*rotz30.getInverse());
+ WALBERLA_CHECK_FLOAT_EQUAL(subb22.getQuaternion(), Quat(real_t(0.3), real_t(0.9),real_t(0.1),real_t(0.3)));
+ WALBERLA_CHECK_FLOAT_EQUAL(subb22.getRelQuaternion(), subb22.getQuaternion()*rotz30.getInverse());
+ // Inertia tensor
+ WALBERLA_CHECK_FLOAT_EQUAL(un2->getInertia(), bodyinertia);
+
+ //Rotate again by 30 deg and perform the same checks as before
+ un2->rotate(rotz30);
+ WALBERLA_CHECK_FLOAT_EQUAL(subb12.getPosition(), Vec3(cos30,sin30,0));
+ WALBERLA_CHECK_FLOAT_EQUAL(subb22.getPosition(), Vec3(-cos30,-sin30,0));
+ WALBERLA_CHECK_FLOAT_EQUAL(subb12.getQuaternion(), Quat(real_t(0.5), real_t(0.5),real_t(0.5),real_t(0.5))*rotz30 );
+ WALBERLA_CHECK_FLOAT_EQUAL(subb22.getQuaternion(), Quat(real_t(0.3), real_t(0.9),real_t(0.1),real_t(0.3))*rotz30 );
+ WALBERLA_CHECK_FLOAT_EQUAL(un->getInertia(), rotz30.toRotationMatrix()*bodyinertia*rotz30.toRotationMatrix().getTranspose());
+
+ WALBERLA_LOG_INFO("- Applying velocities.");
+ //Apply a linear velocity to the union
+ un2->setLinearVel(Vec3(1,0,0));
+ WALBERLA_CHECK_FLOAT_EQUAL(subb12.getLinearVel(), Vec3(1,0,0));
+ WALBERLA_CHECK_FLOAT_EQUAL(subb22.getLinearVel(), Vec3(1,0,0));
+
+ un2->setAngularVel(Vec3(0,0,1));
+ WALBERLA_CHECK_FLOAT_EQUAL(subb12.getAngularVel(), Vec3(0,0,1));
+ WALBERLA_CHECK_FLOAT_EQUAL(subb22.getAngularVel(), Vec3(0,0,1));
+ WALBERLA_CHECK_FLOAT_EQUAL(subb12.getLinearVel(), Vec3(real_t(1.0-sin30), cos30, 0));
+ WALBERLA_CHECK_FLOAT_EQUAL(subb22.getLinearVel(), Vec3(real_t(1.0+sin30), -cos30, 0));
+
+}
+
+void checkAABB(const AABB &a1, const AABB& a2){
+ WALBERLA_CHECK_FLOAT_EQUAL(a1.xMin(), a2.xMin());
+ WALBERLA_CHECK_FLOAT_EQUAL(a1.yMin(), a2.yMin());
+ WALBERLA_CHECK_FLOAT_EQUAL(a1.zMin(), a2.zMin());
+ WALBERLA_CHECK_FLOAT_EQUAL(a1.xMax(), a2.xMax());
+ WALBERLA_CHECK_FLOAT_EQUAL(a1.yMax(), a2.yMax());
+ WALBERLA_CHECK_FLOAT_EQUAL(a1.zMax(), a2.zMax());
+}
+
+void UnionAABB() {
+ WALBERLA_LOG_INFO("- Performing AABB test");
+ MaterialID iron = Material::find("iron");
+ // Construct a union of two spheres in a non-rotated state and check all its parameters
+ auto un = std::make_unique<UnionType>(12, 0, Vec3(0,0,0), Quat(), false, true, false);
+ auto sp1 = std::make_unique<Sphere>( 10, 0, Vec3( 1,0,0), Quat(real_t(0.5), real_t(0.5),real_t(0.5),real_t(0.5)), real_t(1), iron, false, true, false );
+ auto sp2 = std::make_unique<Sphere>( 11, 0, Vec3(-1,0,0), Quat(), real_t(1), iron, false, true, false );
un->add( std::move(sp1) );
un->add( std::move(sp2) );
+ AABB aabb = un->getAABB();
+ WALBERLA_CHECK_FLOAT_EQUAL(aabb.minCorner(), Vec3(-2,-1,-1));
+ WALBERLA_CHECK_FLOAT_EQUAL(aabb.maxCorner(), Vec3(2, 1, 1));
+
+ Vec3 shift(10,10,10);
+ un->setPosition(shift);
+ aabb = un->getAABB();
+ checkAABB(aabb, AABB(real_t(8),real_t(9),real_t(9),real_t(12),real_t(11),real_t(11)));
+
+ Quat rotz30(Vec3(0,0,1), real_t(math::M_PI/6.0)); // rotate by 30 deg via z axis
+ real_t sin30 = real_t(0.5);
+ real_t cos30 = real_t(sqrt(3.0)/2.0);
+ un->setOrientation(rotz30);
+ aabb = un->getAABB();
+ checkAABB(aabb, AABB(real_t(9.0-cos30),real_t(9.0-sin30),real_t(9), real_t(11.0+cos30), real_t(11.0+sin30), real_t(11)));
+
+ aabb = un->begin()->getAABB();
+ checkAABB(aabb, AABB(real_t(9.0+cos30),real_t(9.0+sin30),real_t(9), real_t(11.0+cos30), real_t(11.0+sin30), real_t(11)));
+
+ aabb = (un->begin()+1)->getAABB();
+ checkAABB(aabb, AABB(real_t(9.0-cos30),real_t(9.0-sin30),real_t(9), real_t(11.0-cos30), real_t(11.0-sin30), real_t(11)));
- WALBERLA_CHECK_FLOAT_EQUAL( un->getPosition(), Vec3(0,0,0) );
- WALBERLA_CHECK_FLOAT_EQUAL( un->getLinearVel(), Vec3(0,0,0) );
- WALBERLA_CHECK_FLOAT_EQUAL( un->getAngularVel() * Vec3(1,0,0), real_t(0) );
- WALBERLA_CHECK_FLOAT_EQUAL( un->getAngularVel() * Vec3(0,1,0), real_t(0) );
- WALBERLA_CHECK_GREATER( un->getAngularVel() * Vec3(0,0,1), real_t(0) );
}
int main( int argc, char ** argv )
@@ -126,8 +280,11 @@ int main( int argc, char ** argv )
walberla::MPIManager::instance()->initializeMPI( &argc, &argv );
+ UnionConstruction();
+ UnionAABB();
SnowManFallingOnPlane();
- ImpulsCarryover();
+
+
return EXIT_SUCCESS;
}
@@ -136,4 +293,4 @@ int main( int argc, char ** argv )
int main( int argc, char* argv[] )
{
return walberla::main( argc, argv );
-}
\ No newline at end of file
+}
diff --git a/tests/pe/UnionBehavior.cpp b/tests/pe/UnionBehavior.cpp
new file mode 100644
index 000000000..51037cf48
--- /dev/null
+++ b/tests/pe/UnionBehavior.cpp
@@ -0,0 +1,259 @@
+//======================================================================================================================
+//
+// This file is part of waLBerla. waLBerla is free software: you can
+// redistribute it and/or modify it under the terms of the GNU General Public
+// License as published by the Free Software Foundation, either version 3 of
+// the License, or (at your option) any later version.
+//
+// waLBerla is distributed in the hope that it will be useful, but WITHOUT
+// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+// FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+// for more details.
+//
+// You should have received a copy of the GNU General Public License along
+// with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file UnionBehavior.cpp
+//! \author Tobias Leemann <tobias.leemann@fau.de>
+//! \brief Testcase checking whether a single body and a union representing the same geometry behave equally.
+//======================================================================================================================
+
+//! [Includes]
+#include <pe/basic.h>
+
+#include <core/Environment.h>
+#include <core/grid_generator/HCPIterator.h>
+#include <core/grid_generator/SCIterator.h>
+#include <core/logging/Logging.h>
+#include <core/math/Random.h>
+
+#include "blockforest/StructuredBlockForest.h"
+#include "blockforest/Initialization.h"
+#include "pe/fcd/GenericFCD.h"
+
+#include "pe/rigidbody/BoxFactory.h"
+#include "vtk/VTKOutput.h"
+#include "pe/vtk/BodyVtkOutput.h"
+
+
+#include <iostream>
+#include <sstream>
+#include <fstream>
+#include <cstring>
+#include <tuple>
+
+
+//! [Includes]
+
+namespace walberla {
+using namespace walberla::pe;
+
+
+//! [BodyTypeTuple]
+using BoxUnion = Union<Box>;
+using BodyTypeTuple = std::tuple<Plane, Box, BoxUnion>;
+//! [BodyTypeTuple]
+
+// Override the collide function, so that the collisions between the two bodies themselves are not tracked,
+// only collisions with the surrounding planes are handled
+class PlaneOnlyFCD : public fcd::GenericFCD<BodyTypeTuple, fcd::AnalyticCollideFunctor>{
+public:
+ Contacts& generateContacts(PossibleContacts& possibleContacts) override
+ {
+ contacts_.clear();
+ fcd::AnalyticCollideFunctor<decltype(contacts_)> func(contacts_);
+ for (auto &c : possibleContacts)
+ {
+ if(c.first->getTypeID() == Plane::getStaticTypeID() || c.second->getTypeID() == Plane::getStaticTypeID()) {
+ DoubleCast<BodyTypeTuple, BodyTypeTuple, fcd::AnalyticCollideFunctor<decltype(contacts_)>, bool>::execute(
+ c.first, c.second, func);
+ }
+ }
+ return contacts_;
+ }
+};
+
+
+int main( int argc, char ** argv )
+{
+
+ Environment env(argc, argv);
+ WALBERLA_UNUSED(env);
+
+ // Parameters
+ auto dt = real_t(0.005);
+ int simulationSteps = 4000;
+ auto domainsize = real_t(100);
+ // Size of the box
+ auto boxlenght = real_t(6.0);
+ auto boxwidth = real_t(2.0);
+
+ Vec3 pos(real_t(domainsize/2.0),real_t(domainsize/2.0),real_t(0.8*domainsize));
+ Vec3 v0(10,2,3); // Linear velocity of the objects
+ Vec3 w0(4,15,1); // Some angular velocity of the objects (quite high)
+ int visSpacing = 10;
+
+
+ // Simulation part
+ math::seedRandomGenerator( static_cast<unsigned int>(1337 * mpi::MPIManager::instance()->worldRank()) );
+
+
+
+ //! [Parameters]
+
+ WALBERLA_LOG_INFO_ON_ROOT("*** GLOBALBODYSTORAGE ***");
+ //! [GlobalBodyStorage]
+ shared_ptr<BodyStorage> globalBodyStorage = make_shared<BodyStorage>();
+ //! [GlobalBodyStorage]
+
+ WALBERLA_LOG_INFO_ON_ROOT("*** BLOCKFOREST ***");
+ // create forest
+ //! [BlockForest]
+ auto forest = blockforest::createBlockForest( AABB(0,0,0,domainsize,domainsize,domainsize), // simulation domain
+ Vector3<uint_t>(3,1,1), // blocks in each direction
+ Vector3<bool>(true, false, false) // periodicity
+ );
+
+ //! [BlockForest]
+ if (!forest)
+ {
+ WALBERLA_LOG_INFO_ON_ROOT( "No BlockForest created ... exiting!");
+ return EXIT_SUCCESS;
+ }
+
+ WALBERLA_LOG_INFO_ON_ROOT("*** STORAGEDATAHANDLING ***");
+ // add block data
+ //! [StorageDataHandling]
+ auto storageID = forest->addBlockData(createStorageDataHandling<BodyTypeTuple>(), "Storage");
+ //! [StorageDataHandling]
+ //! [AdditionalBlockData]
+ auto ccdID = forest->addBlockData(ccd::createHashGridsDataHandling( globalBodyStorage, storageID ), "CCD");
+
+ auto sharedFCDBDH = make_shared<blockforest::AlwaysCreateBlockDataHandling<PlaneOnlyFCD> >( );
+
+ auto fcdID = forest->addBlockData(sharedFCDBDH, "FCD");
+ //! [AdditionalBlockData]
+
+ WALBERLA_LOG_INFO_ON_ROOT("*** INTEGRATOR ***");
+ //! [Integrator]
+ cr::DEM cr(globalBodyStorage, forest, storageID, ccdID, fcdID);
+ cr.setGlobalLinearAcceleration( Vec3(0,0,-6) );
+ //! [Integrator]
+
+ WALBERLA_LOG_INFO_ON_ROOT("*** BodyTypeTuple ***");
+ // initialize body type ids
+ //! [SetBodyTypeIDs]
+ SetBodyTypeIDs<BodyTypeTuple>::execute();
+ //! [SetBodyTypeIDs]
+
+
+ // UNCOMMENT THIS BLOCK FOR VTK OUTPUT
+ /*WALBERLA_LOG_INFO_ON_ROOT("*** VTK OUTPUT ***");
+ //! [VTK Domain Output]
+ auto vtkDomainOutput = vtk::createVTKOutput_DomainDecomposition( forest, "domain_decomposition", 1, "VTK", "simulation_step" );
+ vtkDomainOutput->write(true);
+ auto vtkSphereHelper = make_shared<DefaultBodyVTKOutput>(storageID, *forest) ;
+ auto vtkSphereOutput = vtk::createVTKOutput_PointData(vtkSphereHelper, "Bodies", 1, "VTK", "simulation_step", false, false);*/
+ //! [VTK Domain Output]
+
+ WALBERLA_LOG_INFO_ON_ROOT("*** SETUP - START ***");
+ //! [Material]
+ const real_t static_cof ( real_c(1.2) / 2 ); // Coefficient of static friction. Note: pe doubles the input coefficient of friction for material-material contacts.
+ const real_t dynamic_cof ( static_cof ); // Coefficient of dynamic friction. Similar to static friction for low speed friction.
+ MaterialID material = createMaterial( "granular", real_t( 1.0 ), 0, static_cof, dynamic_cof, real_t( 0.5 ), 1, real_t(8.11e5), real_t(6.86e1), real_t(6.86e1) );
+ //! [Material]
+
+ // Create surronding planes in y,z directions, but not in x.
+ auto simulationDomain = forest->getDomain();
+ //! [Planes]
+ createPlane(*globalBodyStorage, 100, Vec3(0,1,0), simulationDomain.minCorner(), material );
+ createPlane(*globalBodyStorage, 101, Vec3(0,-1,0), simulationDomain.maxCorner(), material );
+ createPlane(*globalBodyStorage, 102, Vec3(0,0,1), simulationDomain.minCorner(), material );
+ createPlane(*globalBodyStorage, 103, Vec3(0,0,-1), simulationDomain.maxCorner(), material );
+ //! [Planes]
+
+ //! [Gas]
+ uint_t numParticles = uint_c(0);
+
+ // Create a union of a two half boxes with the center at pos
+ BoxUnion* particleU = createUnion<Box>( *globalBodyStorage, *forest, storageID, 0, Vec3(10,10,10));
+ // add 2 parts to the union
+ createBox(particleU, 4, pos-Vec3(real_t(boxlenght/4.0),0,0), Vec3(real_t(boxlenght/2.0), boxwidth, boxwidth), material);
+ createBox(particleU, 5, pos+Vec3(real_t(boxlenght/4.0),0,0), Vec3(real_t(boxlenght/2.0), boxwidth, boxwidth), material);
+
+ if (particleU != nullptr) particleU->setLinearVel(v0);
+ if (particleU != nullptr) particleU->setAngularVel(w0);
+ if (particleU != nullptr) ++numParticles;
+
+ // Generate the same box as one particle at pos
+ BoxID particleS = createBox( *globalBodyStorage, *forest, storageID, 1, pos, Vec3(boxlenght, boxwidth, boxwidth), material );
+ if (particleS != nullptr) particleS->setLinearVel(v0);
+ if (particleS != nullptr) particleS->setAngularVel(w0);
+ if (particleS != nullptr) ++numParticles;
+
+ WALBERLA_LOG_INFO_ON_ROOT("#particles created per process: " << numParticles);
+ syncNextNeighbors<BodyTypeTuple>(*forest, storageID);
+ //! [Gas]
+
+ WALBERLA_LOG_INFO_ON_ROOT("*** SETUP - END ***");
+
+ WALBERLA_LOG_INFO_ON_ROOT("*** SIMULATION - START ***");
+ //vtkDomainOutput->write(true);
+ //! [GameLoop]
+ for (int i=0; i < simulationSteps; ++i)
+ {
+ if( i % visSpacing == 0 )
+ {
+ WALBERLA_LOG_INFO_ON_ROOT( "Timestep " << i << " / " << simulationSteps );
+ // UNCOMMENT FOR OUTPUT
+ //vtkSphereOutput->write(true);
+ }
+
+ cr.timestep( real_c(dt) );
+
+ int bdcount = 0;
+ RigidBody* refU = nullptr;
+ RigidBody* refS = nullptr;
+ for (auto blockIt = forest->begin(); blockIt != forest->end(); ++blockIt)
+ {
+ for (auto bodyIt = LocalBodyIterator::begin(*blockIt, storageID); bodyIt != LocalBodyIterator::end(); ++bodyIt)
+ {
+ if(bodyIt->getID() == 0){
+ refU = &*bodyIt;
+ }
+ if(bodyIt->getID() == 1){
+ refS = &*bodyIt;
+ }
+ bdcount ++;
+ }
+ }
+
+ // Compare with an quite large epsilon, as round of errors in collisions will amplify themselves over time
+ auto eps = real_t(1e-4);
+ WALBERLA_CHECK_GREATER_EQUAL(bdcount, 2, "Lost a body.");
+ WALBERLA_CHECK_NOT_NULLPTR(refS, "Single body was not found.");
+ WALBERLA_CHECK_NOT_NULLPTR(refU, "Union was not found.");
+ // norm(pos1-pos2) < 1e-4 -> norm*norm < 1e-8
+ WALBERLA_CHECK_LESS((refU->getPosition()-refS->getPosition()).sqrLength(), eps*eps );
+ WALBERLA_CHECK_LESS((refU->getLinearVel()-refS->getLinearVel()).sqrLength(), eps*eps );
+ WALBERLA_CHECK_LESS((refU->getAngularVel()-refS->getAngularVel()).sqrLength(), eps*eps);
+
+ // Sums of matrix differences...
+ WALBERLA_CHECK_LESS(std::abs(Vec3(1,1,1) * ((refU->getRotation()-refS->getRotation())*Vec3(1,1,1))), real_t(9.0*eps));
+ WALBERLA_CHECK_LESS(std::abs(Vec3(1,1,1) * ((refU->getInertia()-refS->getInertia())*Vec3(1,1,1))), real_t(9.0*eps));
+ //WALBERLA_LOG_INFO_ON_ROOT("Timestep " << i << ": " <<refU->getPosition() << "/" << refS->getPosition() << (Vec3(1,1,1) * ((refU->getInertia()-refS->getInertia())*Vec3(1,1,1)))<< " " << abs(Vec3(1,1,1) * ((refU->getRotation()-refS->getRotation())*Vec3(1,1,1))));
+
+
+ syncNextNeighbors<BodyTypeTuple>(*forest, storageID);
+ }
+ //! [GameLoop]
+ WALBERLA_LOG_INFO_ON_ROOT("*** SIMULATION - END ***");
+
+ return EXIT_SUCCESS;
+}
+} // namespace walberla
+
+int main( int argc, char* argv[] )
+{
+ return walberla::main( argc, argv );
+}
diff --git a/tests/pe/VolumeInertia.cpp b/tests/pe/VolumeInertia.cpp
index 6deebadad..2ab7738bb 100644
--- a/tests/pe/VolumeInertia.cpp
+++ b/tests/pe/VolumeInertia.cpp
@@ -97,7 +97,7 @@ int main( int argc, char ** argv )
Vec3 COM;
Mat3 inertia;
- Sphere sp(0, 0, Vec3(0,0,0), Vec3(0,0,0), Quat(), real_t(2.34), material, false, true, false);
+ Sphere sp(0, 0, Vec3(0,0,0), Quat(), real_t(2.34), material, false, true, false);
calcNumeric(sp, sp.getAABB(), real_t(0.01), volume, COM, inertia);
WALBERLA_CHECK( walberla::debug::check_functions_detail::check_float_equal_eps(sp.getVolume(), volume, real_t(10e-4)) );
WALBERLA_CHECK( walberla::debug::check_functions_detail::check_float_equal_eps(Sphere::calcVolume( real_t(2.34) ), volume, real_t(10e-4)) );
@@ -105,7 +105,7 @@ int main( int argc, char ** argv )
WALBERLA_CHECK( walberla::debug::check_functions_detail::check_float_equal_eps(sp.getMass(), volume * Material::getDensity(material), real_t(10e-4)) );
WALBERLA_CHECK( walberla::debug::check_functions_detail::check_float_equal_eps(COM, Vec3(0), real_t(10e-4)) );
- Box bx(0, 0, Vec3(0,0,0), Vec3(0,0,0), Quat(), Vec3(real_t(1.5), real_t(2.5), real_t(3.5)), material, false, true, false);
+ Box bx(0, 0, Vec3(0,0,0), Quat(), Vec3(real_t(1.5), real_t(2.5), real_t(3.5)), material, false, true, false);
calcNumeric(bx, bx.getAABB(), real_t(0.01), volume, COM, inertia);
WALBERLA_CHECK( walberla::debug::check_functions_detail::check_float_equal_eps(bx.getVolume(), volume, real_t(10e-4)) );
WALBERLA_CHECK( walberla::debug::check_functions_detail::check_float_equal_eps(Box::calcVolume( Vec3(real_t(1.5), real_t(2.5), real_t(3.5)) ), volume, real_t(10e-4)) );
@@ -113,7 +113,7 @@ int main( int argc, char ** argv )
WALBERLA_CHECK( walberla::debug::check_functions_detail::check_float_equal_eps(bx.getMass(), volume * Material::getDensity(material), real_t(10e-4)) );
WALBERLA_CHECK( walberla::debug::check_functions_detail::check_float_equal_eps(COM, Vec3(0), real_t(10e-4)) );
- Ellipsoid el(0, 0, Vec3(0,0,0), Vec3(0,0,0), Quat(), Vec3(real_t(1.5), real_t(2.5), real_t(3.5)), material, false, true, false);
+ Ellipsoid el(0, 0, Vec3(0,0,0), Quat(), Vec3(real_t(1.5), real_t(2.5), real_t(3.5)), material, false, true, false);
calcNumeric(el, el.getAABB(), real_t(0.01), volume, COM, inertia);
WALBERLA_CHECK( walberla::debug::check_functions_detail::check_float_equal_eps(el.getVolume(), volume, real_t(10e-4)) );
WALBERLA_CHECK( walberla::debug::check_functions_detail::check_float_equal_eps(Ellipsoid::calcVolume( Vec3(real_t(1.5), real_t(2.5), real_t(3.5)) ), volume, real_t(10e-4)) );
@@ -121,7 +121,7 @@ int main( int argc, char ** argv )
WALBERLA_CHECK( walberla::debug::check_functions_detail::check_float_equal_eps(el.getMass(), volume * Material::getDensity(material), real_t(10e-4)) );
WALBERLA_CHECK( walberla::debug::check_functions_detail::check_float_equal_eps(COM, Vec3(0), real_t(10e-4)) );
- Capsule cp(0, 0, Vec3(0,0,0), Vec3(0,0,0), Quat(), real_t(1.5), real_t(2.5), material, false, true, false);
+ Capsule cp(0, 0, Vec3(0,0,0), Quat(), real_t(1.5), real_t(2.5), material, false, true, false);
calcNumeric(cp, cp.getAABB(), real_t(0.01), volume, COM, inertia);
WALBERLA_CHECK( walberla::debug::check_functions_detail::check_float_equal_eps(cp.getVolume(), volume, real_t(10e-4)) );
WALBERLA_CHECK( walberla::debug::check_functions_detail::check_float_equal_eps(Capsule::calcVolume( real_t(1.5), real_t(2.5) ), volume, real_t(10e-4)) );
diff --git a/tests/pe_coupling/geometry/PeIntersectionRatioTest.cpp b/tests/pe_coupling/geometry/PeIntersectionRatioTest.cpp
index 496f7db51..76e9a79ba 100644
--- a/tests/pe_coupling/geometry/PeIntersectionRatioTest.cpp
+++ b/tests/pe_coupling/geometry/PeIntersectionRatioTest.cpp
@@ -60,8 +60,7 @@ int main( int argc, char **argv )
walberla::id_t sid = 0;
walberla::id_t uid = 0;
-
- Vector3<real_t> rPos( real_t(0));
+
Vector3<real_t> rotationAngles( real_t(0));
Quaternion<real_t> quat( rotationAngles );
pe::MaterialID material = pe::Material::find("iron");
@@ -74,7 +73,7 @@ int main( int argc, char **argv )
Vector3<real_t> bodyPos(real_t(1), real_t(0), real_t(0));
real_t radius = real_t(1);
- pe::Sphere sphere(++sid, ++uid, bodyPos, rPos, quat, radius, material, false, false, false);
+ pe::Sphere sphere(++sid, ++uid, bodyPos, quat, radius, material, false, false, false);
pe::RigidBody & rb = sphere; // otherwise not the pe_coupling/geometry version is matched
@@ -123,7 +122,7 @@ int main( int argc, char **argv )
Vector3<real_t> bodyPos(real_t(1), real_t(0), real_t(0));
Vector3<real_t> semiAxes1(real_t(1), real_t(1), real_t(1));
- pe::Ellipsoid ellip1(++sid, ++uid, bodyPos, rPos, quat, semiAxes1, material, false, false, false);
+ pe::Ellipsoid ellip1(++sid, ++uid, bodyPos, quat, semiAxes1, material, false, false, false);
pe::RigidBody & rb1 = ellip1; // otherwise not the pe_coupling/geometry version is matched
@@ -138,7 +137,7 @@ int main( int argc, char **argv )
WALBERLA_CHECK_FLOAT_EQUAL(delta2, (std::sqrt(2) - real_t(1)) / std::sqrt(2), "Intersection ratio with ellipsoid wrong!");
Vector3<real_t> semiAxes2(real_t(2), real_t(0.5), real_t(2));
- pe::Ellipsoid ellip2(++sid, ++uid, bodyPos, rPos, quat, semiAxes2, material, false, false, false);
+ pe::Ellipsoid ellip2(++sid, ++uid, bodyPos, quat, semiAxes2, material, false, false, false);
pe::RigidBody & rb2 = ellip2; // otherwise not the pe_coupling/geometry version is matched
--
GitLab