diff --git a/src/pe/Types.h b/src/pe/Types.h
index d37661d54fb2d76466a4ae9b3f4cc8fe608456c2..35d559f4f7c19d81a0ba36d04482a9d7639c7c43 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 25bf84ca03855d3dce0a22e1099dc0860d6cfce5..70bbdf4426bbaeb38d75c05fb89da6bea1c7bc9a 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 4777794366f6d6cd9eda691b45b187785ed1da01..8b6fa28538b785034d8a1e94bf371f6940aeb403 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 0000000000000000000000000000000000000000..1a96385517c916ce10aafe20132d840269aa5178
--- /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 0000000000000000000000000000000000000000..4c095c1a20bae85959c756a856500ec7b4332b93
--- /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 0000000000000000000000000000000000000000..a82ea0c684ad7dbc7f60c2d02b99fb7eb6ded2f2
--- /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 0000000000000000000000000000000000000000..9d516f4789cd09dbaf292b21d2856be7de993bba
--- /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
diff --git a/src/pe_coupling/mapping/BodyBBMapping.cpp b/src/pe_coupling/mapping/BodyBBMapping.cpp
index 070434f31e6dfc97476ecaa7d4bdc79e197f46d5..9defb26ecd2fc688ff0432903e38dda314691a16 100644
--- a/src/pe_coupling/mapping/BodyBBMapping.cpp
+++ b/src/pe_coupling/mapping/BodyBBMapping.cpp
@@ -33,8 +33,13 @@ CellInterval getCellBB( const pe::ConstBodyID body, const IBlock & block, const
{
CellInterval cellBB;
- const auto & aabb = body->getAABB();
- blockStorage->getCellBBFromAABB( cellBB, aabb, blockStorage->getLevel(block) );
+ if (body->isFinite())
+ {
+ blockStorage->getCellBBFromAABB( cellBB, body->getAABB(), blockStorage->getLevel(block) );
+ } else
+ {
+ blockStorage->getCellBBFromAABB( cellBB, body->getAABB().getIntersection( blockStorage->getDomain() ), blockStorage->getLevel(block) );
+ }
cellBB.xMin() -= cell_idx_c(1); cellBB.yMin() -= cell_idx_c(1); cellBB.zMin() -= cell_idx_c(1);
cellBB.xMax() += cell_idx_c(1); cellBB.yMax() += cell_idx_c(1); cellBB.zMax() += cell_idx_c(1);
diff --git a/src/pe_coupling/mapping/BodyMapping.h b/src/pe_coupling/mapping/BodyMapping.h
index 8d6bfc09a6f4137bb6595dc2d0440b51d1b986b8..6df7a1341381f868b6c6d46e6134fa48f7fc1f81 100644
--- a/src/pe_coupling/mapping/BodyMapping.h
+++ b/src/pe_coupling/mapping/BodyMapping.h
@@ -48,7 +48,7 @@ void mapBody( const pe::BodyID & body, IBlock & block, const shared_ptr<Structur
WALBERLA_UNUSED(moBodiesOnly); // undo when other coupling algorithms are available
- if( ( fixedBodiesOnly && !body->isFixed() ) || !body->isFinite() /*|| ( moBodiesOnly && !isMOBody( body ) )*/ )
+ if( ( fixedBodiesOnly && !body->isFixed() ) /*|| ( moBodiesOnly && !isMOBody( body ) )*/ )
return;
BoundaryHandling_T * boundaryHandling = block.getData< BoundaryHandling_T >( boundaryHandlingID );
@@ -108,7 +108,24 @@ void mapGlobalBodies( const shared_ptr<StructuredBlockStorage> & blockStorage, c
for( auto blockIt = blockStorage->begin(); blockIt != blockStorage->end(); ++blockIt )
{
for( auto bodyIt = globalBodyStorage.begin(); bodyIt != globalBodyStorage.end(); ++bodyIt)
+ {
+ mapBody< BoundaryHandling_T >( *bodyIt, *blockIt, blockStorage, boundaryHandlingID, obstacle, fixedBodiesOnly, moBodiesOnly );
+ }
+ }
+}
+
+template< typename BoundaryHandling_T >
+void mapGlobalBody( const id_t globalBodySystemID,
+ const shared_ptr<StructuredBlockStorage> & blockStorage, const BlockDataID & boundaryHandlingID, pe::BodyStorage & globalBodyStorage, const FlagUID & obstacle,
+ const bool fixedBodiesOnly = true, const bool moBodiesOnly = true )
+{
+ for( auto blockIt = blockStorage->begin(); blockIt != blockStorage->end(); ++blockIt )
+ {
+ auto bodyIt = globalBodyStorage.find( globalBodySystemID );
+ if( bodyIt != globalBodyStorage.end() )
+ {
mapBody< BoundaryHandling_T >( *bodyIt, *blockIt, blockStorage, boundaryHandlingID, obstacle, fixedBodiesOnly, moBodiesOnly );
+ }
}
}
diff --git a/src/pe_coupling/momentum_exchange_method/BodyMapping.h b/src/pe_coupling/momentum_exchange_method/BodyMapping.h
index 4170551cb8a985d286006df571f30c521b927184..732b31fe67bc483c3a6e52ad6daf85560c3fc726 100644
--- a/src/pe_coupling/momentum_exchange_method/BodyMapping.h
+++ b/src/pe_coupling/momentum_exchange_method/BodyMapping.h
@@ -181,7 +181,7 @@ void mapMovingBody( const pe::BodyID body, IBlock & block, const shared_ptr<Stru
WALBERLA_ASSERT_EQUAL( &block.getBlockStorage(), &(blockStorage->getBlockStorage()) );
- if( body->isFixed() || !body->isFinite() )
+ if( body->isFixed() /*|| !body->isFinite()*/ )
return;
BoundaryHandling_T * boundaryHandling = block.getData< BoundaryHandling_T >( boundaryHandlingID );
@@ -240,6 +240,33 @@ void mapMovingBodies( const shared_ptr<StructuredBlockStorage> & blockStorage, c
}
}
+template< typename BoundaryHandling_T >
+void mapMovingGlobalBodies( const shared_ptr<StructuredBlockStorage> & blockStorage, const BlockDataID & boundaryHandlingID, pe::BodyStorage & globalBodyStorage,
+ const BlockDataID & bodyFieldID, const FlagUID & obstacle )
+{
+ for( auto blockIt = blockStorage->begin(); blockIt != blockStorage->end(); ++blockIt )
+ {
+ for( auto bodyIt = globalBodyStorage.begin(); bodyIt != globalBodyStorage.end(); ++bodyIt)
+ {
+ mapMovingBody< BoundaryHandling_T >( *bodyIt, *blockIt, blockStorage, boundaryHandlingID, bodyFieldID, obstacle );
+ }
+ }
+}
+
+template< typename BoundaryHandling_T >
+void mapMovingGlobalBody( const id_t globalBodySystemID,
+ const shared_ptr<StructuredBlockStorage> & blockStorage, const BlockDataID & boundaryHandlingID, pe::BodyStorage & globalBodyStorage,
+ const BlockDataID & bodyFieldID, const FlagUID & obstacle )
+{
+ for( auto blockIt = blockStorage->begin(); blockIt != blockStorage->end(); ++blockIt )
+ {
+ auto bodyIt = globalBodyStorage.find( globalBodySystemID );
+ if( bodyIt != globalBodyStorage.end() )
+ {
+ mapMovingBody< BoundaryHandling_T >( *bodyIt, *blockIt, blockStorage, boundaryHandlingID, bodyFieldID, obstacle );
+ }
+ }
+}
} // namespace pe_coupling
diff --git a/src/pe_coupling/momentum_exchange_method/boundary/SimpleBB.h b/src/pe_coupling/momentum_exchange_method/boundary/SimpleBB.h
index 161dfc872e1094df917d2394048790afdf52da92..0c7d3aae8c5b7352eb6426ef87d54e0986f54e1b 100644
--- a/src/pe_coupling/momentum_exchange_method/boundary/SimpleBB.h
+++ b/src/pe_coupling/momentum_exchange_method/boundary/SimpleBB.h
@@ -142,7 +142,7 @@ inline void SimpleBB< LatticeModel_T, FlagField_T >::treatDirection( const cell_
WALBERLA_ASSERT_UNEQUAL( mask & this->mask_, numeric_cast<flag_t>(0) );
WALBERLA_ASSERT_EQUAL ( mask & this->mask_, this->mask_ ); // only true if "this->mask_" only contains one single flag, which is the case for the
// current implementation of this boundary condition
- WALBERLA_ASSERT_NOT_NULLPTR( bodyField_->get(nx,ny,nz) );
+ WALBERLA_ASSERT_NOT_NULLPTR( bodyField_->get(nx,ny,nz), "(" << nx << ", " << ny << ", " << nz << ")" );
const real_t pdf_old = pdfField_->get( x, y, z, Stencil_T::idx[dir] );
diff --git a/tests/pe/Collision.cpp b/tests/pe/Collision.cpp
index 9abc2fb2e7080262dda47b15ca7b48c8f18dbd6c..98ce1fc4a43f00854911c54a7373f0c86e50cbf8 100644
--- a/tests/pe/Collision.cpp
+++ b/tests/pe/Collision.cpp
@@ -58,7 +58,9 @@ void SphereTest()
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);
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);
std::vector<Contact> contacts;
fcd::AnalyticCollideFunctor< std::vector<Contact> > collideFunc(contacts);
@@ -85,6 +87,26 @@ void SphereTest()
Contact( &sp2, &pl1, Vec3(1,real_t(-0.5),real_t(-0.5)), Vec3(1, 1, 1).getNormalized(), real_t(-0.133974596215561181)) );
WALBERLA_CHECK( !collideFunc(&sp3, &pl1) );
+
+ // SPHERE <-> CYLINDRICAL BOUNDARY
+ WALBERLA_LOG_INFO("SPHERE <-> CYLINDRICAL BOUNDARY");
+ contacts.clear();
+ WALBERLA_CHECK( !collideFunc(&sp1, &cb1) );
+ WALBERLA_CHECK( !collideFunc(&sp2, &cb1) );
+ WALBERLA_CHECK( collideFunc(&sp4, &cb1) );
+ checkContact( contacts.at(0),
+ Contact( &sp4, &cb1, Vec3(0,real_t(2),real_t(0)), Vec3(0, -1, 0).getNormalized(), real_t(-0.5)) );
+ cb1.rotateAroundOrigin( Vec3( 0,0,1), math::M_PI * real_t(0.25) );
+ WALBERLA_CHECK( !collideFunc(&sp1, &cb1) );
+ WALBERLA_CHECK( collideFunc(&sp2, &cb1) );
+ WALBERLA_CHECK( collideFunc(&sp4, &cb1) );
+ const real_t xPos = real_t(3) / real_t(4) + real_t(2) / real_c(sqrt(real_t(2)));
+ const real_t yPos = xPos - real_t(4) / real_c(sqrt(real_t(2)));
+ const real_t dist = real_c(sqrt((xPos - real_t(1.5)) * (xPos - real_t(1.5)) + yPos * yPos)) - sp4.getRadius();
+ checkContact( contacts.at(1),
+ Contact( &sp2, &cb1, Vec3(xPos, yPos, 0), Vec3(-1, +1, 0).getNormalized(), dist) );
+ checkContact( contacts.at(2),
+ Contact( &sp4, &cb1, Vec3(yPos, xPos, 0), Vec3(+1, -1, 0).getNormalized(), dist) );
}
void BoxTest()
diff --git a/tests/pe_coupling/CMakeLists.txt b/tests/pe_coupling/CMakeLists.txt
index e31ded9e826492cb2926ada9d7272dd4647bccd6..b956c82da8a313411a476bb030f0de629de0955a 100644
--- a/tests/pe_coupling/CMakeLists.txt
+++ b/tests/pe_coupling/CMakeLists.txt
@@ -99,6 +99,10 @@ waLBerla_execute_test( NAME LubricationCorrectionMEMSphereWallTest COMMAND $<T
waLBerla_compile_test( FILES momentum_exchange_method/PeriodicParticleChannelMEM.cpp DEPENDS blockforest pe timeloop )
waLBerla_execute_test( NAME PeriodicParticleChannelMEMTest COMMAND $<TARGET_FILE:PeriodicParticleChannelMEM> --shortrun PROCESSES 4 LABELS longrun )
+waLBerla_compile_test( FILES momentum_exchange_method/TaylorCouetteFlowMEM.cpp DEPENDS blockforest pe timeloop )
+waLBerla_execute_test( NAME TaylorCouetteFlowMEMTest COMMAND $<TARGET_FILE:TaylorCouetteFlowMEM> --shortrun )
+waLBerla_execute_test( NAME TaylorCouetteFlowMEMTest2 COMMAND $<TARGET_FILE:TaylorCouetteFlowMEM> PROCESSES 4 LABELS longrun)
+
waLBerla_compile_test( FILES momentum_exchange_method/SegreSilberbergMEM.cpp DEPENDS blockforest pe timeloop )
waLBerla_execute_test( NAME SegreSilberbergMEMSyncNextNeighbor COMMAND $<TARGET_FILE:SegreSilberbergMEM> --shortrun PROCESSES 9 )
waLBerla_execute_test( NAME SegreSilberbergMEMSyncShadowOwner COMMAND $<TARGET_FILE:SegreSilberbergMEM> --shortrun --syncShadowOwners PROCESSES 9 )
diff --git a/tests/pe_coupling/momentum_exchange_method/PeriodicParticleChannelMEM.cpp b/tests/pe_coupling/momentum_exchange_method/PeriodicParticleChannelMEM.cpp
index 51ae2623cf6413f54c9efa451cebcc94b7f64e18..46e97a1129656e65c26e5383cf9c471fa71f3b59 100644
--- a/tests/pe_coupling/momentum_exchange_method/PeriodicParticleChannelMEM.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/PeriodicParticleChannelMEM.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 PeriodicParticleChannelMEMPe.cpp
+//! \file PeriodicParticleChannelMEM.cpp
//! \ingroup pe_coupling
//! \author Florian Schornbaum <florian.schornbaum@fau.de>
//! \author Christoph Rettinger <christoph.rettinger@fau.de>
@@ -65,8 +65,9 @@
#include "vtk/all.h"
+#include <vector>
-namespace periodic_particle_channel_mem_pe
+namespace periodic_particle_channel_mem
{
///////////
@@ -423,7 +424,7 @@ int main( int argc, char **argv )
// create pe bodies
const auto material = pe::createMaterial( "granular", real_c(1.2), real_c(0.25), real_c(0.4), real_c(0.4), real_c(0.35), real_c(1.39e11), real_c(5.18e7), real_c(1.07e2), real_c(1.07e2) );
- // global (fixed) bodies
+ // global bodies
// bounding planes
pe::createPlane( *globalBodyStorage, 0, Vector3<real_t>(0,0,1), Vector3<real_t>(0,0,0), material );
pe::createPlane( *globalBodyStorage, 0, Vector3<real_t>(0,0,-1), Vector3<real_t>(0,0,real_c(width)), material );
@@ -431,14 +432,18 @@ int main( int argc, char **argv )
pe::createPlane( *globalBodyStorage, 0, Vector3<real_t>(0,-1,0), Vector3<real_t>(0,real_c(width),0), material );
// spheres as obstacles
- pe::createSphere( *globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0, Vector3<real_t>(real_c(length) / real_t(2), real_t(50), real_t(110)), real_t(60), material, false, true, true );
- pe::createSphere( *globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0, Vector3<real_t>( real_t(0), real_t(50), -real_t(60)), real_t(80), material, false, true, true );
- pe::createSphere( *globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0, Vector3<real_t>( real_c(length), real_t(50), -real_t(60)), real_t(80), material, false, true, true );
+ std::vector<walberla::id_t> globalBodiesToBeMapped;
+ auto globalSphere1 = pe::createSphere( *globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0, Vector3<real_t>(real_c(length) / real_t(2), real_t(50), real_t(110)), real_t(60), material, true, false, true );
+ globalBodiesToBeMapped.push_back(globalSphere1->getSystemID() );
+ auto globalSphere2 = pe::createSphere( *globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0, Vector3<real_t>( real_t(0), real_t(50), -real_t(60)), real_t(80), material, true, false, true );
+ globalBodiesToBeMapped.push_back(globalSphere2->getSystemID() );
+ auto globalSphere3 = pe::createSphere( *globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0, Vector3<real_t>( real_c(length), real_t(50), -real_t(60)), real_t(80), material, true, false, true );
+ globalBodiesToBeMapped.push_back(globalSphere3->getSystemID() );
// local bodies: moving spheres
const real_t radius = real_t(10);
- auto sphere = pe::createSphere( *globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0, Vector3<real_t>( real_t(15), real_t(50), real_t(35) ), radius, material ); //TODO change
+ auto sphere = pe::createSphere( *globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0, Vector3<real_t>( real_t(15), real_t(50), real_t(35) ), radius, material );
if( sphere != NULL ) sphere->setLinearVel( velocity, real_t(0), real_t(0) );
sphere = pe::createSphere( *globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0, Vector3<real_t>( real_t(15), real_t(35), real_t(50) ), radius, material );
@@ -494,7 +499,11 @@ int main( int argc, char **argv )
// map pe bodies into the LBM simulation
// global bodies act as no-slip obstacles and are not treated by the fluid-solid coupling
- pe_coupling::mapGlobalBodies< BoundaryHandling_T >( blocks, boundaryHandlingID, *globalBodyStorage, NoSlip_Flag );
+ // special care has to be taken here that only the global spheres (not the planes) are mapped since the planes would overwrite the already set boundary flags
+ for( auto globalBodyIt = globalBodiesToBeMapped.begin(); globalBodyIt != globalBodiesToBeMapped.end(); ++globalBodyIt )
+ {
+ pe_coupling::mapGlobalBody< BoundaryHandling_T >( *globalBodyIt, blocks, boundaryHandlingID, *globalBodyStorage, NoSlip_Flag, false );
+ }
// moving bodies are handled by the momentum exchange method
pe_coupling::mapMovingBodies< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MO_Flag );
@@ -608,8 +617,8 @@ int main( int argc, char **argv )
return 0;
}
-} //namespace periodic_particle_channel_mem_pe
+} //namespace periodic_particle_channel_mem
int main( int argc, char **argv ){
- periodic_particle_channel_mem_pe::main(argc, argv);
+ periodic_particle_channel_mem::main(argc, argv);
}
diff --git a/tests/pe_coupling/momentum_exchange_method/TaylorCouetteFlowMEM.cpp b/tests/pe_coupling/momentum_exchange_method/TaylorCouetteFlowMEM.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b7d4fc43c277766caaacf7f03a129577dbc79bd1
--- /dev/null
+++ b/tests/pe_coupling/momentum_exchange_method/TaylorCouetteFlowMEM.cpp
@@ -0,0 +1,421 @@
+//======================================================================================================================
+//
+// 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 TaylorCouetteFlowMEM.cpp
+//! \ingroup pe_coupling
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#include "blockforest/Initialization.h"
+#include "blockforest/communication/UniformBufferedScheme.h"
+
+#include "boundary/all.h"
+
+#include "core/DataTypes.h"
+#include "core/Environment.h"
+#include "core/SharedFunctor.h"
+#include "core/debug/Debug.h"
+#include "core/debug/TestSubsystem.h"
+#include "core/math/all.h"
+#include "core/timing/RemainingTimeLogger.h"
+
+#include "domain_decomposition/SharedSweep.h"
+
+#include "field/AddToStorage.h"
+#include "field/communication/PackInfo.h"
+#include "field/vtk/all.h"
+
+#include "lbm/boundary/all.h"
+#include "lbm/communication/PdfFieldPackInfo.h"
+#include "lbm/field/AddToStorage.h"
+#include "lbm/field/PdfField.h"
+#include "lbm/lattice_model/D3Q19.h"
+#include "lbm/sweeps/CellwiseSweep.h"
+#include "lbm/vtk/all.h"
+
+#include "timeloop/SweepTimeloop.h"
+
+#include "pe/basic.h"
+#include "pe/Types.h"
+
+#include "pe_coupling/mapping/all.h"
+#include "pe_coupling/momentum_exchange_method/all.h"
+
+#include "vtk/all.h"
+
+namespace taylor_coette_flow_mem
+{
+
+///////////
+// USING //
+///////////
+
+using namespace walberla;
+using walberla::uint_t;
+
+
+
+//////////////
+// TYPEDEFS //
+//////////////
+
+// pdf field & flag field
+typedef lbm::D3Q19< lbm::collision_model::TRT, false > LatticeModel_T;
+typedef LatticeModel_T::Stencil Stencil_T;
+typedef lbm::PdfField< LatticeModel_T > PdfField_T;
+
+typedef walberla::uint8_t flag_t;
+typedef FlagField< flag_t > FlagField_T;
+typedef GhostLayerField< pe::BodyID, 1 > BodyField_T;
+
+const uint_t FieldGhostLayers = 1;
+
+typedef pe_coupling::CurvedLinear< LatticeModel_T, FlagField_T > MO_T;
+
+typedef boost::tuples::tuple< MO_T > BoundaryConditions_T;
+typedef BoundaryHandling< FlagField_T, Stencil_T, BoundaryConditions_T > BoundaryHandling_T;
+
+typedef boost::tuple< pe::Capsule, pe::CylindricalBoundary > BodyTypeTuple;
+
+///////////
+// FLAGS //
+///////////
+
+const FlagUID Fluid_Flag( "fluid" );
+const FlagUID MO_Flag( "moving obstacle" );
+const FlagUID FormerMO_Flag( "former moving obstacle" );
+
+
+
+/////////////////////
+// CHECK FUNCTIONS //
+/////////////////////
+
+class VelocityChecker
+{
+public:
+ VelocityChecker( const shared_ptr< StructuredBlockStorage > & blocks, const BlockDataID & pdfFieldID,
+ real_t radius1, real_t radius2, real_t angularVel1, real_t angularVel2, real_t domainLength, real_t domainWidth ) :
+ blocks_( blocks ), pdfFieldID_( pdfFieldID ), radius1_( radius1 ), radius2_( radius2 ), angularVel1_( angularVel1 ), angularVel2_( angularVel2 ),
+ domainLength_( domainLength ), domainWidth_( domainWidth )
+ { }
+ real_t getMaximumRelativeError()
+ {
+ Vector3<real_t> midPoint( domainLength_ * real_t(0.5), domainWidth_ * real_t(0.5), domainWidth_ * real_t(0.5) );
+ Cell midCell = blocks_->getCell( midPoint );
+ CellInterval evaluationCells( midCell.x(), midCell.y(), midCell.z(), midCell.x(), blocks_->getDomainCellBB().yMax(), midCell.z() );
+
+ real_t maxError = real_t(0);
+
+ for( auto blockIt = blocks_->begin(); blockIt != blocks_->end(); ++blockIt )
+ {
+ auto pdfField = blockIt->getData< PdfField_T > ( pdfFieldID_ );
+ CellInterval intersection = blocks_->getBlockCellBB( *blockIt );
+ intersection.intersect( evaluationCells );
+ blocks_->transformGlobalToBlockLocalCellInterval( intersection, *blockIt );
+ for( auto fieldIt = pdfField->beginSliceXYZ(intersection); fieldIt != pdfField->end(); ++fieldIt )
+ {
+ const Vector3< real_t > cellCenter = blocks_->getBlockLocalCellCenter( *blockIt, fieldIt.cell() );
+ auto diff = cellCenter - midPoint;
+ real_t radius = diff.length();
+ if( radius > radius1_ && radius < radius2_ )
+ {
+ real_t velSimulation = pdfField->getVelocity( fieldIt.cell() )[2];
+ real_t velAnalytical = getAnalyticalVelocity( radius );
+ real_t relVelDiff = std::fabs(velSimulation - velAnalytical) / std::fabs(velAnalytical);
+ maxError = std::max( maxError, relVelDiff );
+ }
+ }
+ }
+ WALBERLA_MPI_SECTION()
+ {
+ mpi::allReduceInplace( maxError, mpi::MAX );
+ }
+ return maxError;
+ }
+private:
+ real_t getAnalyticalVelocity( real_t r )
+ {
+ real_t radius1sqr = radius1_ * radius1_;
+ real_t radius2sqr = radius2_ * radius2_;
+ real_t a = ( angularVel2_ * radius2sqr - angularVel1_ * radius1sqr ) / ( radius2sqr - radius1sqr );
+ real_t b = ( angularVel1_ - angularVel2_ ) * radius1sqr * radius2sqr / ( radius2sqr - radius1sqr );
+ real_t velInAgularDir = a * r + b / r;
+ return velInAgularDir;
+ }
+
+ shared_ptr< StructuredBlockStorage > blocks_;
+ const BlockDataID pdfFieldID_;
+ real_t radius1_, radius2_, angularVel1_, angularVel2_, domainLength_, domainWidth_;
+};
+
+/////////////////////////////////////
+// BOUNDARY HANDLING CUSTOMIZATION //
+/////////////////////////////////////
+
+// class MyBoundaryHandling is responsible for creating the boundary handling and setting up the geometry at the outer domain boundaries
+
+class MyBoundaryHandling
+{
+public:
+
+ MyBoundaryHandling( const BlockDataID & flagFieldID, const BlockDataID & pdfFieldID, const BlockDataID & bodyFieldID ) :
+ flagFieldID_( flagFieldID ), pdfFieldID_( pdfFieldID ), bodyFieldID_( bodyFieldID ) {}
+
+ BoundaryHandling_T * operator()( IBlock* const block, const StructuredBlockStorage* const storage ) const;
+
+private:
+
+ const BlockDataID flagFieldID_;
+ const BlockDataID pdfFieldID_;
+ const BlockDataID bodyFieldID_;
+
+}; // class MyBoundaryHandling
+
+
+
+BoundaryHandling_T * MyBoundaryHandling::operator()( IBlock * const block, const StructuredBlockStorage * const storage ) const
+{
+ WALBERLA_ASSERT_NOT_NULLPTR( block );
+ WALBERLA_ASSERT_NOT_NULLPTR( storage );
+
+ FlagField_T * flagField = block->getData< FlagField_T >( flagFieldID_ );
+ PdfField_T * pdfField = block->getData< PdfField_T > ( pdfFieldID_ );
+ BodyField_T * bodyField = block->getData< BodyField_T >( bodyFieldID_ );
+
+ const auto fluid = flagField->flagExists( Fluid_Flag ) ? flagField->getFlag( Fluid_Flag ) : flagField->registerFlag( Fluid_Flag );
+
+ BoundaryHandling_T * handling = new BoundaryHandling_T( "cf boundary handling", flagField, fluid,
+ boost::tuples::make_tuple( MO_T( "MO", MO_Flag, pdfField, flagField, bodyField, fluid, *storage, *block ) ) );
+
+ CellInterval domainBB = storage->getDomainCellBB();
+ domainBB.xMin() -= cell_idx_c( FieldGhostLayers );
+ domainBB.xMax() += cell_idx_c( FieldGhostLayers );
+
+ domainBB.yMin() -= cell_idx_c( FieldGhostLayers );
+ domainBB.yMax() += cell_idx_c( FieldGhostLayers );
+
+ domainBB.zMin() -= cell_idx_c( FieldGhostLayers );
+ domainBB.zMax() += cell_idx_c( FieldGhostLayers );
+
+ handling->fillWithDomain( FieldGhostLayers );
+
+ return handling;
+}
+
+
+//////////
+// MAIN //
+//////////
+/*
+ * This test features a so-called Taylor-Couette flow, i.e. a Couette flow between two rotating, coaxial cylinders
+ * For such a setup, the velocity in angular direction can be calculated analytically.
+ */
+int main( int argc, char **argv )
+{
+ debug::enterTestMode();
+
+ mpi::Environment env( argc, argv );
+
+ auto processes = MPIManager::instance()->numProcesses();
+ if( processes != 1 && processes % 4 != 0 )
+ {
+ std::cerr << "Number of processes must be equal to either 1 or a multiple of 4!" << std::endl;
+ return EXIT_FAILURE;
+ }
+
+ bool shortrun = false;
+ uint_t vtkIOFreq = 0;
+
+ for( int i = 1; i < argc; ++i )
+ {
+ if( std::strcmp( argv[i], "--shortrun" ) == 0 ) { shortrun = true; continue; }
+ if( std::strcmp( argv[i], "--vtkIOFreq" ) == 0 ) { vtkIOFreq = uint_c( std::atof( argv[++i] ) ); continue; }
+ WALBERLA_ABORT("Unrecognized command line argument found: " << argv[i]);
+ }
+
+ bool vtkIO = (vtkIOFreq == 0 ) ? false : true;
+
+ ///////////////////////////
+ // SIMULATION PROPERTIES //
+ ///////////////////////////
+
+ const uint_t length = uint_c(20);
+ const uint_t width = uint_c(101);
+ if( width % 2 == 0 ) WALBERLA_ABORT("Width has to be odd, since otherwise the current implementation of the velocity evaluation does not work correctly!");
+
+ const real_t omega = real_c(1.33);
+
+ const uint_t timesteps = shortrun ? uint_c(10) : uint_c(5000);
+
+ const real_t radius1 = real_c(width) * real_t(0.25); // radius of internal cylinder
+ const real_t radius2 = real_c(width) * real_t(0.45); // radius of external cylinder
+ const real_t angularVel1 = real_t(0.001); // angular velocity of internal cylinder
+ const real_t angularVel2 = real_t(-0.001); // angular velocity of external cylinder
+
+ ///////////////////////////
+ // BLOCK STRUCTURE SETUP //
+ ///////////////////////////
+
+ // (length x width x width) - periodic in x-direction!
+ Vector3<uint_t> blockDist = getFactors3D( uint_c(processes), Vector3<real_t>( real_c(length), real_c(width), real_c(width) ) );
+ if ( processes == 1 )
+ {
+ // due to periodicity, the pe needs at least three blocks in the periodic direction
+ blockDist[0] = uint_t(4);
+ }
+
+ const uint_t xCells = length / blockDist[0];
+ const uint_t yCells = width / blockDist[1];
+ const uint_t zCells = width / blockDist[2];
+
+ const real_t dx = real_t(1);
+
+ auto blocks = blockforest::createUniformBlockGrid( blockDist[0], blockDist[1], blockDist[2], xCells, yCells, zCells, dx,
+ ( processes != 1 ),
+ true, false, false );
+ if( vtkIO )
+ {
+ vtk::writeDomainDecomposition( blocks );
+ }
+ /////////////////
+ // PE COUPLING //
+ /////////////////
+
+ // set up pe functionality
+ shared_ptr<pe::BodyStorage> globalBodyStorage = make_shared<pe::BodyStorage>();
+ pe::SetBodyTypeIDs<BodyTypeTuple>::execute();
+ auto bodyStorageID = blocks->addBlockData(pe::createStorageDataHandling<BodyTypeTuple>(), "pe Body Storage");
+
+ // set up synchronization procedure
+ const real_t overlap = real_t( 1.5 ) * dx;
+ boost::function<void(void)> syncCall = boost::bind( pe::syncShadowOwners<BodyTypeTuple>, boost::ref(blocks->getBlockForest()), bodyStorageID, static_cast<WcTimingTree*>(NULL), overlap, false );
+
+ // create pe bodies
+ const auto material = pe::createMaterial( "granular", real_t(1.2), real_t(0.25), real_t(0.4), real_t(0.4), real_t(0.35), real_t(1.39e11), real_t(5.18e7), real_t(1.07e2), real_t(1.07e2) );
+
+ // instead of a cylinder, we use the capsule and make sure it extends the computational domain in x-direction to effectively get a cylinder
+ auto cap = pe::createCapsule( *globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0, pe::Vec3(real_c(length),real_c(width),real_c(width)) * real_t(0.5), radius1, real_c(length)*real_t(2), material, true, false, true );
+ cap->setAngularVel( pe::Vec3(1,0,0) * angularVel1 );
+
+ auto cb = pe::createCylindricalBoundary( *globalBodyStorage, 0, pe::Vec3(real_c(length),real_c(width),real_c(width))*real_t(0.5), radius2 );
+ cb->setAngularVel( pe::Vec3(1,0,0) * angularVel2 );
+
+ //synchronize the pe set up on all processes
+ syncCall();
+
+ ////////////////////////
+ // ADD DATA TO BLOCKS //
+ ////////////////////////
+
+ // add pdf field
+ LatticeModel_T latticeModel = LatticeModel_T( lbm::collision_model::TRT::constructWithMagicNumber( omega ) );
+
+ BlockDataID pdfFieldID = lbm::addPdfFieldToStorage( blocks, "pdf field (zyxf)", latticeModel, Vector3< real_t >( real_t(0) ), real_t(1), uint_t(1), field::zyxf );
+
+ // add flag field
+ BlockDataID flagFieldID = field::addFlagFieldToStorage<FlagField_T>( blocks, "flag field" );
+
+ // add body field
+ BlockDataID bodyFieldID = field::addToStorage<BodyField_T>( blocks, "body field", NULL, field::zyxf );
+
+ // add boundary handling & initialize outer domain boundaries (moving walls on the front, back, top, and bottom plane)
+ BlockDataID boundaryHandlingID = blocks->addStructuredBlockData< BoundaryHandling_T >(
+ MyBoundaryHandling( flagFieldID, pdfFieldID, bodyFieldID ), "boundary handling" );
+
+ // map pe bodies into the LBM simulation
+ // moving bodies are handled by the momentum exchange method, thus act here as velocity boundary condition
+ pe_coupling::mapMovingGlobalBodies< BoundaryHandling_T >( blocks, boundaryHandlingID, *globalBodyStorage, bodyFieldID, MO_Flag );
+
+ ///////////////
+ // TIME LOOP //
+ ///////////////
+
+ SweepTimeloop timeloop( blocks->getBlockStorage(), timesteps );
+
+ // setup of the LBM communication for synchronizing the pdf field between neighboring blocks
+ boost::function< void () > commFunction;
+ blockforest::communication::UniformBufferedScheme< Stencil_T > scheme( blocks );
+ scheme.addPackInfo( make_shared< lbm::PdfFieldPackInfo< LatticeModel_T > >( pdfFieldID ) );
+ commFunction = scheme;
+
+ // add LBM communication function and boundary handling sweep
+ timeloop.add()
+ << BeforeFunction( commFunction, "LBM Communication" )
+ << Sweep( BoundaryHandling_T::getBlockSweep( boundaryHandlingID ), "Boundary Handling" );
+
+ // LBM stream collide sweep
+ timeloop.add() << Sweep( makeSharedSweep( lbm::makeCellwiseSweep< LatticeModel_T, FlagField_T >( pdfFieldID, flagFieldID, Fluid_Flag ) ), "LBM DEFAULT" );
+
+ ////////////////
+ // VTK OUTPUT //
+ ////////////////
+
+ if( vtkIO )
+ {
+ const uint_t writeFrequency = 10; //vtk::determineWriteFrequency( dt_SI, uint_c(30) );
+
+ // flag field (written only once in the first time step, ghost layers are also written)
+ auto flagFieldVTK = vtk::createVTKOutput_BlockData( blocks, "flag_field", timesteps, FieldGhostLayers );
+ flagFieldVTK->addCellDataWriter( make_shared< field::VTKWriter< FlagField_T > >( flagFieldID, "FlagField" ) );
+ timeloop.addFuncAfterTimeStep( vtk::writeFiles( flagFieldVTK ), "VTK (flag field data)" );
+
+ // pdf field (ghost layers cannot be written because re-sampling/coarsening is applied)
+ auto pdfFieldVTK = vtk::createVTKOutput_BlockData( blocks, "fluid_field", writeFrequency );
+ pdfFieldVTK->setSamplingResolution( real_c(1) );
+
+ blockforest::communication::UniformBufferedScheme< stencil::D3Q27 > pdfGhostLayerSync( blocks );
+ pdfGhostLayerSync.addPackInfo( make_shared< field::communication::PackInfo< PdfField_T > >( pdfFieldID ) );
+ pdfFieldVTK->addBeforeFunction( pdfGhostLayerSync );
+
+ field::FlagFieldCellFilter< FlagField_T > fluidFilter( flagFieldID );
+ fluidFilter.addFlag( Fluid_Flag );
+ pdfFieldVTK->addCellInclusionFilter( fluidFilter );
+
+ pdfFieldVTK->addCellDataWriter( make_shared< lbm::VelocityVTKWriter< LatticeModel_T, float > >( pdfFieldID, "VelocityFromPDF" ) );
+ pdfFieldVTK->addCellDataWriter( make_shared< lbm::DensityVTKWriter < LatticeModel_T, float > >( pdfFieldID, "DensityFromPDF" ) );
+
+ timeloop.addFuncAfterTimeStep( vtk::writeFiles( pdfFieldVTK ), "VTK (fluid field data)" );
+ }
+
+ timeloop.addFuncAfterTimeStep( RemainingTimeLogger( timeloop.getNrOfTimeSteps() ), "Remaining Time Logger" );
+
+ ////////////////////////
+ // EXECUTE SIMULATION //
+ ////////////////////////
+
+ WcTimingPool timeloopTiming;
+ timeloop.run( timeloopTiming );
+
+ VelocityChecker velChecker(blocks, pdfFieldID, radius1, radius2, angularVel1, angularVel2, real_c(length), real_c(width) );
+ real_t maxError = velChecker.getMaximumRelativeError();
+ if( !shortrun )
+ {
+ timeloopTiming.logResultOnRoot();
+ WALBERLA_LOG_RESULT_ON_ROOT("Maximum relative error in velocity in angular direction: " << maxError );
+ // error has to be below 10%
+ WALBERLA_CHECK_LESS( maxError, real_t(0.1) );
+ }
+
+ return 0;
+}
+
+} //namespace taylor_coette_flow_mem
+
+int main( int argc, char **argv ){
+ taylor_coette_flow_mem::main(argc, argv);
+}