From 1b558276f45997a9c658b219353e4317ebc98a22 Mon Sep 17 00:00:00 2001 From: Sebastian Eibl <sebastian.eibl@fau.de> Date: Mon, 2 Oct 2017 18:03:04 +0200 Subject: [PATCH] added cylindrical boundary --- src/pe/Types.h | 5 + src/pe/basic.h | 1 + src/pe/fcd/AnalyticCollisionDetection.h | 48 +++- src/pe/rigidbody/CylindricalBoundary.cpp | 239 ++++++++++++++++++ src/pe/rigidbody/CylindricalBoundary.h | 186 ++++++++++++++ .../rigidbody/CylindricalBoundaryFactory.cpp | 51 ++++ src/pe/rigidbody/CylindricalBoundaryFactory.h | 56 ++++ 7 files changed, 580 insertions(+), 6 deletions(-) create mode 100644 src/pe/rigidbody/CylindricalBoundary.cpp create mode 100644 src/pe/rigidbody/CylindricalBoundary.h create mode 100644 src/pe/rigidbody/CylindricalBoundaryFactory.cpp create mode 100644 src/pe/rigidbody/CylindricalBoundaryFactory.h diff --git a/src/pe/Types.h b/src/pe/Types.h index d37661d54..35d559f4f 100644 --- a/src/pe/Types.h +++ b/src/pe/Types.h @@ -65,6 +65,7 @@ class Box; class Capsule; class Contact; class Cylinder; +class CylindricalBoundary; class FixedJoint; class ForceGenerator; class GeomPrimitive; @@ -106,6 +107,10 @@ typedef Cylinder CylinderType; //!< Type of the cylinder geo typedef Cylinder* CylinderID; //!< Handle for a cylinder primitive. typedef const Cylinder* ConstCylinderID; //!< Handle for a constant cylinder primitive. +typedef CylindricalBoundary CylindricalBoundaryType; //!< Type of the cylindrical boundary geometric primitive. +typedef CylindricalBoundary* CylindricalBoundaryID; //!< Handle for a cylindrical boundary primitive. +typedef const CylindricalBoundary* ConstCylindricalBoundaryID; //!< Handle for a constant cylindrical boundary primitive. + typedef Plane PlaneType; //!< Type of the plane geometric primitive. typedef Plane* PlaneID; //!< Handle for a plane primitive. typedef const Plane* ConstPlaneID; //!< Handle for a constant plane primitive. diff --git a/src/pe/basic.h b/src/pe/basic.h index 25bf84ca0..70bbdf442 100644 --- a/src/pe/basic.h +++ b/src/pe/basic.h @@ -40,6 +40,7 @@ #include "pe/rigidbody/BodyIterators.h" #include "pe/rigidbody/BoxFactory.h" #include "pe/rigidbody/CapsuleFactory.h" +#include "pe/rigidbody/CylindricalBoundaryFactory.h" #include "pe/rigidbody/PlaneFactory.h" #include "pe/rigidbody/SphereFactory.h" #include "pe/rigidbody/UnionFactory.h" diff --git a/src/pe/fcd/AnalyticCollisionDetection.h b/src/pe/fcd/AnalyticCollisionDetection.h index 477779436..8b6fa2853 100644 --- a/src/pe/fcd/AnalyticCollisionDetection.h +++ b/src/pe/fcd/AnalyticCollisionDetection.h @@ -26,6 +26,7 @@ #include "pe/contact/Contact.h" #include "pe/rigidbody/Box.h" #include "pe/rigidbody/Capsule.h" +#include "pe/rigidbody/CylindricalBoundary.h" #include "pe/rigidbody/Plane.h" #include "pe/rigidbody/Sphere.h" #include "pe/rigidbody/Union.h" @@ -54,34 +55,42 @@ bool collide( SphereID s1, SphereID s2, Container& container ); template <typename Container> inline bool collide( SphereID s, PlaneID p, Container& container ); - template <typename Container> inline bool collide( PlaneID p, SphereID s, Container& container ); + template <typename Container> inline -bool collide( SphereID s, BoxID b, Container& container ); +bool collide( SphereID s, CylindricalBoundaryID cb, Container& container ); +template <typename Container> +inline +bool collide( CylindricalBoundaryID cb, SphereID s, Container& container ); +template <typename Container> +inline +bool collide( SphereID s, BoxID b, Container& container ); template <typename Container> inline bool collide( BoxID b, SphereID s, Container& container ); + template <typename Container> inline bool collide( BoxID b1, BoxID b2, Container& container ); + template <typename Container> inline bool collide( BoxID b, PlaneID p, Container& container ); - template <typename Container> inline bool collide( PlaneID p, BoxID b, Container& container ); + template <typename Container> inline bool collide( CapsuleID c1, CapsuleID c2, Container& container ); + template <typename Container> inline bool collide( CapsuleID c, PlaneID p, Container& container ); - template <typename Container> inline bool collide( PlaneID p, CapsuleID c, Container& container ); @@ -89,7 +98,6 @@ bool collide( PlaneID p, CapsuleID c, Container& container ); template <typename Container> inline bool collide( SphereID s, CapsuleID c, Container& container ); - template <typename Container> inline bool collide( CapsuleID c, SphereID s, Container& container ); @@ -98,7 +106,6 @@ bool collide( CapsuleID c, SphereID s, Container& container ); template <typename Container> inline bool collide( BoxID b, CapsuleID c, Container& container ); - template <typename Container> inline bool collide( CapsuleID c, BoxID b, Container& container ); @@ -219,6 +226,35 @@ bool collide( PlaneID p, SphereID s, Container& container ) return collide(s, p, container); } +template <typename Container> +inline +bool collide( SphereID s, CylindricalBoundaryID cb, Container& container ) +{ + WALBERLA_ASSERT_GREATER( cb->getRadius(), s->getRadius() ); + + const Vec3 dist = (s->getPosition() - cb->getPosition()) - ((s->getPosition() - cb->getPosition()) * cb->getAxis()) * cb->getAxis(); + const real_t effRadius = cb->getRadius() - s->getRadius(); + if( effRadius * effRadius - dist.sqrLength() < contactThreshold ) { + const Vec3 contactNormal = -dist.getNormalized(); + const real_t penetrationDepth = effRadius - dist.length(); + const Vec3 contactPoint = ( s->getPosition() - ( s->getRadius() + penetrationDepth ) * contactNormal ); + + WALBERLA_LOG_DETAIL( " Contact created between sphere " << s->getID() + << " and cylindrical boundary " << cb->getID() << " (dist=" << penetrationDepth << ")"); + + container.push_back( Contact(s, cb, contactPoint, contactNormal, penetrationDepth) ); + return true; + } + return false; +} + +template <typename Container> +inline +bool collide( CylindricalBoundaryID cb, SphereID s, Container& container ) +{ + return collide(s, cb, container); +} + //************************************************************************************************* /*!\brief Contact generation between a Sphere and a Box. * diff --git a/src/pe/rigidbody/CylindricalBoundary.cpp b/src/pe/rigidbody/CylindricalBoundary.cpp new file mode 100644 index 000000000..1a9638551 --- /dev/null +++ b/src/pe/rigidbody/CylindricalBoundary.cpp @@ -0,0 +1,239 @@ +//====================================================================================================================== +// +// 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 CylindricalBoundary.cpp +//! \author Sebastian Eibl <sebastian.eibl@fau.de> +// +//====================================================================================================================== + +#include "CylindricalBoundary.h" + +//************************************************************************************************* +// Includes +//************************************************************************************************* + +#include <iomanip> +#include <iostream> +#include <stdexcept> +#include "core/math/Shims.h" +#include <pe/Materials.h> +#include <core/math/Matrix3.h> +#include <core/math/Limits.h> +#include <core/debug/Debug.h> + +namespace walberla { +namespace pe { + +//================================================================================================= +// +// CONSTRUCTOR +// +//================================================================================================= + +//************************************************************************************************* +/*!\brief Constructor for the CylindricalBoundary class. + * + * \param sid Unique system-specific ID for the cylindrical boundary. + * \param uid User-specific ID for the cylindrical boundary. + * \param gpos Global geometric center of the cylindrical boundary. + * \param radius The radius of the cylinder. + * \param material The material of the cylindrical boundary. + * + * The cylindrical boundary is created lying along the x-axis. + */ +CylindricalBoundary::CylindricalBoundary( id_t sid, id_t uid, const Vec3& gpos, const real_t radius, + MaterialID material ) + : GeomPrimitive( getStaticTypeID(), sid, uid, material ) // Initializing the base object + , radius_(radius) // Radius of the cylinder // Length of the capsule +{ + //boundaries are always considered locally and have infinite mass + setGlobal( true ); + setCommunicating( false ); + setMass( true ); + setFinite( false ); + + // Checking the radius + // Since the constructor is never directly called but only used in a small number + // of functions that already check the capsule arguments, only asserts are used here to + // double check the arguments. + WALBERLA_ASSERT_GREATER( radius, real_t(0), "Invalid capsule radius" ); + + // Initializing the instantiated capsule + gpos_ = gpos; + q_ = Quat(); // Setting the orientation + R_ = q_.toRotationMatrix(); // Setting the rotation matrix + + // Setting the axis-aligned bounding box + CylindricalBoundary::calcBoundingBox(); +} +//************************************************************************************************* + + + + +//************************************************************************************************* +/*!\brief Destructor for the cylindrical boundary class. + */ +CylindricalBoundary::~CylindricalBoundary() +{ + WALBERLA_LOG_DETAIL( "Destroyed CylindricalBoundary " << sid_ ); +} +//************************************************************************************************* + + + + +//================================================================================================= +// +// UTILITY FUNCTIONS +// +//================================================================================================= + +//************************************************************************************************* +/*!\brief Checks, whether a point in body relative coordinates lies inside the cylindrical boundary. + * + * \param px The x-component of the relative coordinate. + * \param py The y-component of the relative coordinate. + * \param pz The z-component of the relative coordinate. + * \return \a true if the point lies inside the cylindrical boundary, \a false if not. + */ +bool CylindricalBoundary::containsRelPointImpl( real_t /*px*/, real_t py, real_t pz ) const +{ + return ( math::sq(py) + math::sq(pz) ) >= ( radius_ * radius_ ); +} +//************************************************************************************************* + + +//************************************************************************************************* +/*!\brief Checks, whether a point in body relative coordinates lies on the surface of the cylindrical boundary. + * + * \param px The x-component of the relative coordinate. + * \param py The y-component of the relative coordinate. + * \param pz The z-component of the relative coordinate. + * \return \a true if the point lies on the surface of the cylindrical boundary, \a false if not. + * + * The tolerance level of the check is pe::surfaceThreshold. + */ +bool CylindricalBoundary::isSurfaceRelPointImpl( real_t /*px*/, real_t py, real_t pz ) const +{ + return ( std::fabs( ( math::sq(py) + math::sq(pz) ) - ( radius_ * radius_ ) ) <= surfaceThreshold*surfaceThreshold ); +} +//************************************************************************************************* + +//************************************************************************************************* +/*!\brief Calculation of the bounding box of the cylindrical boundary. + * + * \return (-inf, -inf, -inf, +inf, +inf, +inf ) + */ +void CylindricalBoundary::calcBoundingBox() +{ + aabb_ = math::AABB( + -math::Limits<real_t>::inf(), + -math::Limits<real_t>::inf(), + -math::Limits<real_t>::inf(), + +math::Limits<real_t>::inf(), + +math::Limits<real_t>::inf(), + +math::Limits<real_t>::inf()); + + WALBERLA_ASSERT( aabb_.contains( gpos_ ), "Invalid bounding box detected" ); +} +//************************************************************************************************* + + +//================================================================================================= +// +// OUTPUT FUNCTIONS +// +//================================================================================================= + +//************************************************************************************************* +/*!\brief Output of the current state of a cylindrical boundary. + * + * \param os Reference to the output stream. + * \param tab Indentation in front of every line of the cylindrical boundary output. + * \return void + */ +void CylindricalBoundary::print( std::ostream& os, const char* tab ) const +{ + using std::setw; + + os << tab << " Cylindrical Boundary " << getID() << " with radius " << getRadius() << "\n"; + + os << tab << " Fixed: " << isFixed() << " , sleeping: " << !isAwake() << "\n"; + + os << tab << " System ID = " << getSystemID() << "\n" + << tab << " Total mass = " << getMass() << "\n" + << tab << " Material = " << Material::getName( material_ ) << "\n" + << tab << " Global position = " << getPosition() << "\n" + << tab << " Linear velocity = " << getLinearVel() << "\n" + << tab << " Angular velocity = " << getAngularVel() << "\n"; + + 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"; +} +//************************************************************************************************* + +//================================================================================================= +// +// GLOBAL OPERATORS +// +//================================================================================================= + +//************************************************************************************************* +/*!\brief Global output operator for cylindrical boundaryies. + * + * \param os Reference to the output stream. + * \param c Reference to a cylindrical boundary object. + * \return Reference to the output stream. + */ +std::ostream& operator<<( std::ostream& os, const CylindricalBoundary& c ) +{ + os << "--" << "CYLINDRICAL BOUNDARY PARAMETERS" + << "------------------------------------------------------------\n"; + c.print( os, "" ); + os << "--------------------------------------------------------------------------------\n" + << std::endl; + return os; +} +//************************************************************************************************* + + +//************************************************************************************************* +/*!\brief Global output operator for cylindrical boundaryies handles. + * + * \param os Reference to the output stream. + * \param c Constant cylindrical boundary handle. + * \return Reference to the output stream. + */ +std::ostream& operator<<( std::ostream& os, CylindricalBoundaryID c ) +{ + os << "--" << "CYLINDRICAL BOUNDARY PARAMETERS" + << "------------------------------------------------------------\n"; + c->print( os, "" ); + os << "--------------------------------------------------------------------------------\n" + << std::endl; + return os; +} +//************************************************************************************************* + +id_t CylindricalBoundary::staticTypeID_ = std::numeric_limits<id_t>::max(); + +} // namespace pe +} // namespace walberla + + diff --git a/src/pe/rigidbody/CylindricalBoundary.h b/src/pe/rigidbody/CylindricalBoundary.h new file mode 100644 index 000000000..4c095c1a2 --- /dev/null +++ b/src/pe/rigidbody/CylindricalBoundary.h @@ -0,0 +1,186 @@ +//====================================================================================================================== +// +// 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 CylindricalBoundary.h +//! \author Sebastian Eibl <sebastian.eibl@fau.de> +// +//====================================================================================================================== + +#pragma once + + +//************************************************************************************************* +// Includes +//************************************************************************************************* + +#include <pe/rigidbody/GeomPrimitive.h> +#include <pe/Types.h> +#include <core/math/Constants.h> +#include <core/math/Matrix3.h> +#include <core/math/Vector3.h> +#include <core/DataTypes.h> +#include <core/logging/Logging.h> +#include <pe/Config.h> + + +namespace walberla { +namespace pe { + +//================================================================================================= +// +// CLASS DEFINITION +// +//================================================================================================= + +//************************************************************************************************* +/** + * \ingroup pe + */ +class CylindricalBoundary : public GeomPrimitive +{ +public: + //**Constructors******************************************************************************** + /*!\name Constructors */ + //@{ + explicit CylindricalBoundary( id_t sid, id_t uid, const Vec3& gpos, const real_t radius, + MaterialID material ); + //@} + //********************************************************************************************** + + //**Destructor********************************************************************************** + /*!\name Destructor */ + //@{ + virtual ~CylindricalBoundary(); + //@} + //********************************************************************************************** + +public: + //**Get functions******************************************************************************* + /*!\name Get functions */ + //@{ + inline real_t getRadius() const; + inline Vec3 getAxis() const; + //@} + //********************************************************************************************** + + //**Utility functions*************************************************************************** + /*!\name Utility functions */ + //@{ + static inline id_t getStaticTypeID(); + //@} + //********************************************************************************************** + + //**Output functions**************************************************************************** + /*!\name Output functions */ + //@{ + virtual void print( std::ostream& os, const char* tab ) const; + //@} + //********************************************************************************************** + +protected: + //**Utility functions*************************************************************************** + /*!\name Utility functions */ + //@{ + virtual bool containsRelPointImpl ( real_t px, real_t py, real_t pz ) const; + virtual bool isSurfaceRelPointImpl( real_t px, real_t py, real_t pz ) const; + + virtual void calcBoundingBox(); // Calculation of the axis-aligned bounding box + //@} + //********************************************************************************************** + + //**Member variables**************************************************************************** + /*!\name Member variables */ + //@{ + real_t radius_; //!< The radius of the cylinder. + //@} + //********************************************************************************************** + +private: + static id_t staticTypeID_; //< type id of sphere, will be set by SetBodyTypeIDs + static void setStaticTypeID(id_t typeID) {staticTypeID_ = typeID;} + + //** friend declaration + /// needed to be able to set static type ids with setStaticTypeID + template <class T> + friend struct SetBodyTypeIDs; +}; +//************************************************************************************************* + + + +//================================================================================================= +// +// GET FUNCTIONS +// +//================================================================================================= + +//************************************************************************************************* +/*!\brief Returns the radius of the cylinder. + * + * \return The radius of the cylinder. + */ +inline real_t CylindricalBoundary::getRadius() const +{ + return radius_; +} +//************************************************************************************************* + +//************************************************************************************************* +/*!\brief Returns the main axis of the cylinder. + * + * \return The main axis of the cylinder. + */ +inline Vec3 CylindricalBoundary::getAxis() const +{ + return getRotation() * Vec3(1,0,0); +} +//************************************************************************************************* + + +//================================================================================================= +// +// UTILITY FUNCTIONS +// +//================================================================================================= + +//************************************************************************************************* +/*!\brief Returns unique type id of this type + * + * \return geometry specific type id + */ +inline id_t CylindricalBoundary::getStaticTypeID() +{ + return staticTypeID_; +} + + + + +//================================================================================================= +// +// GLOBAL OPERATORS +// +//================================================================================================= + +//************************************************************************************************* +/*!\name Capsule operators */ +//@{ +std::ostream& operator<<( std::ostream& os, const CylindricalBoundary& cb ); +std::ostream& operator<<( std::ostream& os, CylindricalBoundaryID cb ); +//@} +//************************************************************************************************* + +} // namespace pe +} // namespace walberla diff --git a/src/pe/rigidbody/CylindricalBoundaryFactory.cpp b/src/pe/rigidbody/CylindricalBoundaryFactory.cpp new file mode 100644 index 000000000..a82ea0c68 --- /dev/null +++ b/src/pe/rigidbody/CylindricalBoundaryFactory.cpp @@ -0,0 +1,51 @@ +//====================================================================================================================== +// +// 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 CylindricalBoundaryFactory.cpp +//! \author Sebastian Eibl <sebastian.eibl@fau.de> +// +//====================================================================================================================== + +#include "CylindricalBoundaryFactory.h" + +#include "pe/Materials.h" +#include "pe/rigidbody/BodyStorage.h" +#include "pe/rigidbody/CylindricalBoundary.h" + +#include "core/logging/Logging.h" + +namespace walberla { +namespace pe { + +CylindricalBoundaryID createCylindricalBoundary( BodyStorage& globalStorage, + id_t uid, const Vec3& gpos, const real_t radius, + MaterialID material) +{ + WALBERLA_ASSERT_UNEQUAL( CylindricalBoundary::getStaticTypeID(), std::numeric_limits<id_t>::max(), "CylindricalBoundary TypeID not initalized!"); + + const id_t sid( UniqueID<RigidBody>::createGlobal() ); + + CylindricalBoundaryID cb = new CylindricalBoundary( sid, uid, gpos, radius, material ); + + globalStorage.add(cb); + + // Logging the successful creation of the plane + WALBERLA_LOG_DETAIL( "Created cylindrical boundary " << sid << "\n" << cb); + + return cb; +} + +} // namespace pe +} // namespace walberla diff --git a/src/pe/rigidbody/CylindricalBoundaryFactory.h b/src/pe/rigidbody/CylindricalBoundaryFactory.h new file mode 100644 index 000000000..9d516f478 --- /dev/null +++ b/src/pe/rigidbody/CylindricalBoundaryFactory.h @@ -0,0 +1,56 @@ +//====================================================================================================================== +// +// 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 CylindricalBoundaryFactory.h +//! \author Sebastian Eibl <sebastian.eibl@fau.de> +// +//====================================================================================================================== + +#pragma once + +#include "pe/rigidbody/BodyStorage.h" +#include "pe/Materials.h" +#include "pe/Types.h" + +#include "core/debug/Debug.h" + +namespace walberla { +namespace pe { + +//================================================================================================= +// +// CYLINDRICAL BOUNDARY SETUP FUNCTIONS +// +//================================================================================================= + +//************************************************************************************************* +/** + * \ingroup pe + * \brief Setup of a new Cylindrical Boundary. + * + * \param globalStorage process local global storage + * \param uid The user-specific ID. + * \param gpos One point located on the central axis. + * \param radius radius of the cylinder + * \param material The material of the boundary. + * \return Handle for the new boundary. + */ +CylindricalBoundaryID createCylindricalBoundary( BodyStorage& globalStorage, + id_t uid, const Vec3& gpos, const real_t radius, + MaterialID material = Material::find("iron")); +//************************************************************************************************* + +} // namespace pe +} // namespace walberla -- GitLab