From 4395f3dff44045bb0484bedf24f060a83eb53cba Mon Sep 17 00:00:00 2001
From: Michael Kuron <mkuron@icp.uni-stuttgart.de>
Date: Wed, 14 Jun 2017 13:50:44 +0200
Subject: [PATCH] Squirmer rigid body type
---
src/pe/Types.h | 5 +
src/pe/communication/DynamicMarshalling.h | 1 +
src/pe/communication/rigidbody/Squirmer.cpp | 59 +++++++
src/pe/communication/rigidbody/Squirmer.h | 78 ++++++++++
src/pe/rigidbody/RigidBody.cpp | 75 +++++++++
src/pe/rigidbody/RigidBody.h | 76 +--------
src/pe/rigidbody/Sphere.cpp | 7 +-
src/pe/rigidbody/Sphere.h | 3 +
src/pe/rigidbody/Squirmer.cpp | 164 ++++++++++++++++++++
src/pe/rigidbody/Squirmer.h | 130 ++++++++++++++++
src/pe/rigidbody/SquirmerFactory.cpp | 84 ++++++++++
src/pe/rigidbody/SquirmerFactory.h | 74 +++++++++
tests/pe/Marshalling.cpp | 22 ++-
13 files changed, 704 insertions(+), 74 deletions(-)
create mode 100644 src/pe/communication/rigidbody/Squirmer.cpp
create mode 100644 src/pe/communication/rigidbody/Squirmer.h
create mode 100644 src/pe/rigidbody/Squirmer.cpp
create mode 100644 src/pe/rigidbody/Squirmer.h
create mode 100644 src/pe/rigidbody/SquirmerFactory.cpp
create mode 100644 src/pe/rigidbody/SquirmerFactory.h
diff --git a/src/pe/Types.h b/src/pe/Types.h
index 15d386d54..d37661d54 100644
--- a/src/pe/Types.h
+++ b/src/pe/Types.h
@@ -82,6 +82,7 @@ class Section;
class SliderJoint;
class Sphere;
class Spring;
+class Squirmer;
class TriangleMesh;
template <typename BodyTypeTuple>
class Union;
@@ -113,6 +114,10 @@ typedef Sphere SphereType; //!< Type of the sphere geome
typedef Sphere* SphereID; //!< Handle for a sphere primitive.
typedef const Sphere* ConstSphereID; //!< Handle for a constant sphere primitive.
+typedef Squirmer SquirmerType; //!< Type of the squirmer geometric primitive.
+typedef Squirmer* SquirmerID; //!< Handle for a squirmer primitive.
+typedef const Squirmer* ConstSquirmerID; //!< Handle for a constant squirmer primitive.
+
typedef TriangleMesh MeshType; //!< Type of the triangle mesh geometric primitive.
typedef TriangleMesh* MeshID; //!< Handle for a triangle mesh primitive.
typedef const TriangleMesh* ConstMeshID; //!< Handle for a constant triangle mesh primitive.
diff --git a/src/pe/communication/DynamicMarshalling.h b/src/pe/communication/DynamicMarshalling.h
index 121931250..735792226 100644
--- a/src/pe/communication/DynamicMarshalling.h
+++ b/src/pe/communication/DynamicMarshalling.h
@@ -32,6 +32,7 @@
#include "pe/communication/rigidbody/Box.h"
#include "pe/communication/rigidbody/Capsule.h"
#include "pe/communication/rigidbody/Sphere.h"
+#include "pe/communication/rigidbody/Squirmer.h"
#include "pe/communication/rigidbody/Union.h"
#include "core/Abort.h"
diff --git a/src/pe/communication/rigidbody/Squirmer.cpp b/src/pe/communication/rigidbody/Squirmer.cpp
new file mode 100644
index 000000000..d4b790b99
--- /dev/null
+++ b/src/pe/communication/rigidbody/Squirmer.cpp
@@ -0,0 +1,59 @@
+//======================================================================================================================
+//
+// 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 Squirmer.cpp
+//! \author Michael Kuron <mkuron@icp.uni-stuttgart.de>
+//! \brief Marshalling of objects for data transmission or storage.
+//
+//======================================================================================================================
+
+#include "Squirmer.h"
+
+namespace walberla {
+namespace pe {
+namespace communication {
+
+//*************************************************************************************************
+/*!\brief Marshalling a squirmer primitive.
+ *
+ * \param buffer The buffer to be filled.
+ * \param obj The object to be marshalled.
+ * \return void
+ */
+void marshal( mpi::SendBuffer& buffer, const Squirmer& obj ) {
+ marshal( buffer, static_cast<const Sphere&>( obj ) );
+ buffer << obj.getSquirmerVelocity();
+ buffer << obj.getSquirmerBeta();
+}
+//*************************************************************************************************
+
+
+//*************************************************************************************************
+/*!\brief Unmarshalling a squirmer primitive.
+ *
+ * \param buffer The buffer from where to read.
+ * \param objparam The object to be reconstructed.
+ * \return void
+ */
+void unmarshal( mpi::RecvBuffer& buffer, SquirmerParameters& objparam ) {
+ unmarshal( buffer, static_cast<SphereParameters&>( objparam ) );
+ buffer >> objparam.squirmerVelocity_;
+ buffer >> objparam.squirmerBeta_;
+}
+//*************************************************************************************************
+
+} // namespace communication
+} // namespace pe
+} // namespace walberla
diff --git a/src/pe/communication/rigidbody/Squirmer.h b/src/pe/communication/rigidbody/Squirmer.h
new file mode 100644
index 000000000..9029c9358
--- /dev/null
+++ b/src/pe/communication/rigidbody/Squirmer.h
@@ -0,0 +1,78 @@
+//======================================================================================================================
+//
+// 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 Squirmer.h
+//! \author Michael Kuron <mkuron@icp.uni-stuttgart.de>
+//! \brief Marshalling of objects for data transmission or storage.
+//
+//======================================================================================================================
+
+#pragma once
+
+
+//*************************************************************************************************
+// Includes
+//*************************************************************************************************
+
+#include "pe/communication/Instantiate.h"
+#include "pe/communication/Marshalling.h"
+#include "pe/rigidbody/Squirmer.h"
+#include "pe/communication/rigidbody/Sphere.h"
+
+namespace walberla {
+namespace pe {
+namespace communication {
+
+struct SquirmerParameters : public SphereParameters {
+ real_t squirmerVelocity_;
+ real_t squirmerBeta_;
+};
+
+//*************************************************************************************************
+/*!\brief Marshalling a squirmer primitive.
+ *
+ * \param buffer The buffer to be filled.
+ * \param obj The object to be marshalled.
+ * \return void
+ */
+void marshal( mpi::SendBuffer& buffer, const Squirmer& obj );
+//*************************************************************************************************
+
+//*************************************************************************************************
+/*!\brief Unmarshalling a squirmer primitive.
+ *
+ * \param buffer The buffer from where to read.
+ * \param objparam The object to be reconstructed.
+ * \return void
+ */
+void unmarshal( mpi::RecvBuffer& buffer, SquirmerParameters& objparam );
+//*************************************************************************************************
+
+
+inline SquirmerID instantiate( mpi::RecvBuffer& buffer, const math::AABB& domain, const math::AABB& block, SquirmerID& newBody )
+{
+ SquirmerParameters subobjparam;
+ unmarshal( buffer, subobjparam );
+ correctBodyPosition(domain, block.center(), subobjparam.gpos_);
+ newBody = new Squirmer( subobjparam.sid_, subobjparam.uid_, subobjparam.gpos_, subobjparam.rpos_, subobjparam.q_, subobjparam.radius_, subobjparam.squirmerVelocity_, subobjparam.squirmerBeta_, subobjparam.material_, false, subobjparam.communicating_, subobjparam.infiniteMass_ );
+ newBody->setLinearVel( subobjparam.v_ );
+ newBody->setAngularVel( subobjparam.w_ );
+ newBody->MPITrait.setOwner( subobjparam.mpiTrait_.owner_ );
+ return newBody;
+}
+
+} // namespace communication
+} // namespace pe
+} // namespace walberla
diff --git a/src/pe/rigidbody/RigidBody.cpp b/src/pe/rigidbody/RigidBody.cpp
index c67cac62c..45d40109e 100644
--- a/src/pe/rigidbody/RigidBody.cpp
+++ b/src/pe/rigidbody/RigidBody.cpp
@@ -183,6 +183,81 @@ bool RigidBody::checkInvariants()
+//=================================================================================================
+//
+// GET FUNCTIONS
+//
+//=================================================================================================
+
+
+//*************************************************************************************************
+/*!\brief Calculation of the global velocity of a relative point.
+ *
+ * \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 The global velocity.
+ *
+ * The function calculates the global velocity of a point relative to the body's center of mass.
+ */
+const Vec3 RigidBody::velFromBF( real_t px, real_t py, real_t pz ) const
+{
+ return velFromBF( Vec3( px, py, pz ) );
+}
+//*************************************************************************************************
+
+
+//*************************************************************************************************
+/*!\brief Calculation of the global velocity of a relative point.
+ *
+ * \param rpos The relative coordinate.
+ * \return The global velocity.
+ *
+ * The function calculates the global velocity of a point relative to the body's center of mass.
+ */
+const Vec3 RigidBody::velFromBF( const Vec3& rpos ) const
+{
+ if( !hasSuperBody() )
+ return v_ + w_ % ( R_ * rpos );
+ else return sb_->velFromWF( gpos_ + R_ * rpos );
+}
+//*************************************************************************************************
+
+
+//*************************************************************************************************
+/*!\brief Calculation of the global velocity of a point in global coordinates.
+ *
+ * \param px The x-component of the global coordinate.
+ * \param py The y-component of the global coordinate.
+ * \param pz The z-component of the global coordinate.
+ * \return The global velocity.
+ *
+ * The function calculates the global velocity of a point in global coordinates.
+ */
+const Vec3 RigidBody::velFromWF( real_t px, real_t py, real_t pz ) const
+{
+ return velFromWF( Vec3( px, py, pz ) );
+}
+//*************************************************************************************************
+
+
+//*************************************************************************************************
+/*!\brief Calculation of the global velocity of a point in global coordinates.
+ *
+ * \param gpos The global coordinate.
+ * \return The global velocity.
+ *
+ * The function calculates the global velocity of a point in global coordinates.
+ */
+const Vec3 RigidBody::velFromWF( const Vec3& gpos ) const
+{
+ if( !hasSuperBody() )
+ return v_ + w_ % ( gpos - gpos_ );
+ else return sb_->velFromWF( gpos );
+}
+//*************************************************************************************************
+
+
//=================================================================================================
//
diff --git a/src/pe/rigidbody/RigidBody.h b/src/pe/rigidbody/RigidBody.h
index 408750a52..c58df20ee 100644
--- a/src/pe/rigidbody/RigidBody.h
+++ b/src/pe/rigidbody/RigidBody.h
@@ -143,14 +143,14 @@ public:
inline const Vec3 vectorFromBFtoWF( const Vec3& v ) const;
inline const Vec3 pointFromBFtoWF ( real_t px, real_t py, real_t pz ) const;
inline const Vec3 pointFromBFtoWF ( const Vec3& rpos ) const;
- inline const Vec3 velFromBF ( real_t px, real_t py, real_t pz ) const;
- inline const Vec3 velFromBF ( const Vec3& rpos ) const;
+ virtual const Vec3 velFromBF ( real_t px, real_t py, real_t pz ) const;
+ virtual const Vec3 velFromBF ( const Vec3& rpos ) const;
inline const Vec3 vectorFromWFtoBF( real_t vx, real_t vy, real_t vz ) const;
inline const Vec3 vectorFromWFtoBF( const Vec3& v ) const;
inline const Vec3 pointFromWFtoBF ( real_t px, real_t py, real_t pz ) const;
inline const Vec3 pointFromWFtoBF ( const Vec3& gpos ) const;
- inline const Vec3 velFromWF ( real_t px, real_t py, real_t pz ) const;
- inline const Vec3 velFromWF ( const Vec3& gpos ) const;
+ virtual const Vec3 velFromWF ( real_t px, real_t py, real_t pz ) const;
+ virtual const Vec3 velFromWF ( const Vec3& gpos ) const;
inline const Vec3 accFromWF ( real_t px, real_t py, real_t pz ) const;
const Vec3 accFromWF ( const Vec3& gpos ) const;
@@ -1083,40 +1083,6 @@ inline const Vec3 RigidBody::pointFromBFtoWF( const Vec3& rpos ) const
//*************************************************************************************************
-//*************************************************************************************************
-/*!\brief Calculation of the global velocity of a relative point.
- *
- * \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 The global velocity.
- *
- * The function calculates the global velocity of a point relative to the body's center of mass.
- */
-inline const Vec3 RigidBody::velFromBF( real_t px, real_t py, real_t pz ) const
-{
- return velFromBF( Vec3( px, py, pz ) );
-}
-//*************************************************************************************************
-
-
-//*************************************************************************************************
-/*!\brief Calculation of the global velocity of a relative point.
- *
- * \param rpos The relative coordinate.
- * \return The global velocity.
- *
- * The function calculates the global velocity of a point relative to the body's center of mass.
- */
-inline const Vec3 RigidBody::velFromBF( const Vec3& rpos ) const
-{
- if( !hasSuperBody() )
- return v_ + w_ % ( R_ * rpos );
- else return sb_->velFromWF( gpos_ + R_ * rpos );
-}
-//*************************************************************************************************
-
-
//*************************************************************************************************
/*!\brief Transformation from a global to a relative vector.
*
@@ -1185,40 +1151,6 @@ inline const Vec3 RigidBody::pointFromWFtoBF( const Vec3& gpos ) const
//*************************************************************************************************
-//*************************************************************************************************
-/*!\brief Calculation of the global velocity of a point in global coordinates.
- *
- * \param px The x-component of the global coordinate.
- * \param py The y-component of the global coordinate.
- * \param pz The z-component of the global coordinate.
- * \return The global velocity.
- *
- * The function calculates the global velocity of a point in global coordinates.
- */
-inline const Vec3 RigidBody::velFromWF( real_t px, real_t py, real_t pz ) const
-{
- return velFromWF( Vec3( px, py, pz ) );
-}
-//*************************************************************************************************
-
-
-//*************************************************************************************************
-/*!\brief Calculation of the global velocity of a point in global coordinates.
- *
- * \param gpos The global coordinate.
- * \return The global velocity.
- *
- * The function calculates the global velocity of a point in global coordinates.
- */
-inline const Vec3 RigidBody::velFromWF( const Vec3& gpos ) const
-{
- if( !hasSuperBody() )
- return v_ + w_ % ( gpos - gpos_ );
- else return sb_->velFromWF( gpos );
-}
-//*************************************************************************************************
-
-
//*************************************************************************************************
/*!\brief Calculation of the global acceleration of a point in global coordinates.
*
diff --git a/src/pe/rigidbody/Sphere.cpp b/src/pe/rigidbody/Sphere.cpp
index 78ceb34de..8ea3a81a7 100644
--- a/src/pe/rigidbody/Sphere.cpp
+++ b/src/pe/rigidbody/Sphere.cpp
@@ -61,7 +61,12 @@ namespace pe {
Sphere::Sphere( id_t sid, id_t uid, const Vec3& gpos, const Vec3& rpos, const Quat& q,
real_t radius, MaterialID material,
const bool global, const bool communicating, const bool infiniteMass )
- : GeomPrimitive( getStaticTypeID(), sid, uid, material ) // Initialization of the parent class
+ : Sphere::Sphere( getStaticTypeID(), sid, uid, gpos, rpos, 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,
+ real_t radius, MaterialID material,
+ const bool global, const bool communicating, const bool infiniteMass )
+ : GeomPrimitive( typeId, sid, uid, material ) // Initialization of the parent class
, radius_(radius)
{
// Checking the radius
diff --git a/src/pe/rigidbody/Sphere.h b/src/pe/rigidbody/Sphere.h
index edfa5acfb..bc0e82566 100644
--- a/src/pe/rigidbody/Sphere.h
+++ b/src/pe/rigidbody/Sphere.h
@@ -66,6 +66,9 @@ public:
explicit Sphere( id_t sid, id_t uid, const Vec3& gpos, const Vec3& rpos, 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,
+ real_t radius, MaterialID material,
+ const bool global, const bool communicating, const bool infiniteMass );
//@}
//**********************************************************************************************
diff --git a/src/pe/rigidbody/Squirmer.cpp b/src/pe/rigidbody/Squirmer.cpp
new file mode 100644
index 000000000..9bd809802
--- /dev/null
+++ b/src/pe/rigidbody/Squirmer.cpp
@@ -0,0 +1,164 @@
+//======================================================================================================================
+//
+// 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 Squirmer.cpp
+//! \author Christian Burkard <buch@icp.uni-stuttgart.de>
+//! \author Michael Kuron <mkuron@icp.uni-stuttgart.de>
+//
+//======================================================================================================================
+
+#include "Squirmer.h"
+
+namespace walberla {
+namespace pe {
+
+//=================================================================================================
+//
+// CONSTRUCTOR
+//
+//=================================================================================================
+
+
+//*************************************************************************************************
+//*************************************************************************************************
+/*!\brief Constructor for the Squirmer class.
+ *
+ * \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 squirmerVelocity The velocity of the squirmer.
+ * \param squirmerBeta The dipolar characteristic of the squirmer.
+ * \param material The material of the sphere.
+ * \param global specifies if the sphere should be created in the global storage
+ * \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,
+ 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
+ , squirmerVelocity_(squirmerVelocity), squirmerBeta_(squirmerBeta)
+{
+}
+
+//=================================================================================================
+//
+// DESTRUCTOR
+//
+//=================================================================================================
+
+//*************************************************************************************************
+/*!\brief Destructor for the Squirmer class.
+ */
+Squirmer::~Squirmer()
+{
+ // Logging the destruction of the squirmer
+ WALBERLA_LOG_DETAIL( "Destroyed squirmer " << sid_ );
+}
+//*************************************************************************************************
+
+
+//=================================================================================================
+//
+// GET FUNCTIONS
+//
+//=================================================================================================
+
+//*************************************************************************************************
+/*!\brief Calculation of the global velocity of a relative point.
+ *
+ * \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 The global velocity.
+ *
+ * The function calculates the global velocity of a point relative to the body's center of mass.
+ */
+const Vec3 Squirmer::velFromBF( real_t px, real_t py, real_t pz ) const
+{
+ return velFromBF( Vec3( px, py, pz ) );
+}
+//*************************************************************************************************
+
+
+//*************************************************************************************************
+/*!\brief Calculation of the global velocity of a relative point.
+ *
+ * \param rpos The relative coordinate.
+ * \return The global velocity.
+ *
+ * The function calculates the global velocity of a point relative to the body's center of mass.
+ */
+const Vec3 Squirmer::velFromBF( const Vec3& rpos ) const
+{
+ return Sphere::velFromBF( rpos ) + getSquirmerVelocity( rpos );
+}
+//*************************************************************************************************
+
+//*************************************************************************************************
+/*!\brief Calculation of the global velocity of a point in global coordinates.
+ *
+ * \param px The x-component of the global coordinate.
+ * \param py The y-component of the global coordinate.
+ * \param pz The z-component of the global coordinate.
+ * \return The global velocity.
+ *
+ * The function calculates the global velocity of a point in global coordinates.
+ */
+const Vec3 Squirmer::velFromWF( real_t px, real_t py, real_t pz ) const
+{
+ return velFromWF( Vec3( px, py, pz ) );
+}
+//*************************************************************************************************
+
+
+//*************************************************************************************************
+/*!\brief Returns the squirmer's surface velocity at a given relative position.
+ *
+ * \return The local surface velocity of the squirmer.
+ */
+const Vec3 Squirmer::getSquirmerVelocity( const Vec3& rpos ) const
+{
+ const auto rs = rpos.getNormalized();
+ const auto e = getQuaternion().rotate(Vec3(0.0,0.0,1.0)).getNormalized();
+ const auto v0 = getSquirmerVelocity();
+ const auto beta = getSquirmerBeta();
+ return ((e * rs) * rs - e) * real_c(1.5) * v0 * (1 + beta * (e * rs));
+}
+//*************************************************************************************************
+
+
+//*************************************************************************************************
+/*!\brief Calculation of the global velocity of a point in global coordinates.
+ *
+ * \param gpos The global coordinate.
+ * \return The global velocity.
+ *
+ * The function calculates the global velocity of a point in global coordinates.
+ */
+const Vec3 Squirmer::velFromWF( const Vec3& gpos ) const
+{
+ const auto rpos = pointFromWFtoBF( gpos );
+ return Sphere::velFromWF( gpos ) + getSquirmerVelocity( rpos );
+}
+//*************************************************************************************************
+
+id_t Squirmer::staticTypeID_ = std::numeric_limits<id_t>::max();
+
+} // namespace pe
+} // namespace walberla
diff --git a/src/pe/rigidbody/Squirmer.h b/src/pe/rigidbody/Squirmer.h
new file mode 100644
index 000000000..3dc76acad
--- /dev/null
+++ b/src/pe/rigidbody/Squirmer.h
@@ -0,0 +1,130 @@
+//======================================================================================================================
+//
+// 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 Squirmer.h
+//! \author Michael Kuron <mkuron@icp.uni-stuttgart.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "Sphere.h"
+
+namespace walberla {
+namespace pe {
+
+class Squirmer: public Sphere
+{
+
+public:
+ //**Constructors********************************************************************************
+ /*!\name Constructors */
+ //@{
+ explicit Squirmer( id_t sid, id_t uid, const Vec3& gpos, const Vec3& rpos, const Quat& q,
+ real_t radius, real_t squirmerVelocity, real_t squirmerBeta, MaterialID material,
+ const bool global, const bool communicating, const bool infiniteMass );
+ //@}
+ //**********************************************************************************************
+
+ //**Destructor**********************************************************************************
+ /*!\name Destructor */
+ //@{
+ virtual ~Squirmer();
+ //@}
+ //**********************************************************************************************
+ //**********************************************************************************************
+
+public:
+ //**Get functions*******************************************************************************
+ /*!\name Get functions */
+ //@{
+ inline real_t getSquirmerVelocity() const;
+ inline real_t getSquirmerBeta() const;
+
+ inline const Vec3 velFromBF ( real_t px, real_t py, real_t pz ) const WALBERLA_OVERRIDE;
+ inline const Vec3 velFromBF ( const Vec3& rpos ) const WALBERLA_OVERRIDE;
+ inline const Vec3 velFromWF ( real_t px, real_t py, real_t pz ) const WALBERLA_OVERRIDE;
+ inline const Vec3 velFromWF ( const Vec3& gpos ) const WALBERLA_OVERRIDE;
+
+ static inline id_t getStaticTypeID() WALBERLA_OVERRIDE;
+ //@}
+ //**********************************************************************************************
+
+protected:
+
+ //**Get functions*******************************************************************************
+ /*!\name Get functions */
+ //@{
+ const Vec3 getSquirmerVelocity( const Vec3& rpos ) const;
+ //@}
+
+ //**Member variables****************************************************************************
+ /*!\name Member variables */
+ //@{
+ real_t squirmerVelocity_; //!< Velocity of the squirmer.
+ real_t squirmerBeta_; //!< Dipolar character of the squirmer.
+ //@}
+ //**********************************************************************************************
+
+private:
+ static id_t staticTypeID_; //< type id of sphere, will be set by SetBodyTypeIDs
+ static void setStaticTypeID(id_t typeID) WALBERLA_OVERRIDE {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 squirmer's velocity.
+ *
+ * \return The velocity of the squirmer.
+ */
+inline real_t Squirmer::getSquirmerVelocity() const
+{
+ return squirmerVelocity_;
+}
+//*************************************************************************************************
+
+//*************************************************************************************************
+/*!\brief Returns the squirmer's dipole characteristic.
+ *
+ * \return The dipole characteristic of the squirmer.
+ */
+inline real_t Squirmer::getSquirmerBeta() const
+{
+ return squirmerBeta_;
+}
+//*************************************************************************************************
+
+//*************************************************************************************************
+/*!\brief Returns unique type id of this type
+ *
+ * \return geometry specific type id
+ */
+inline id_t Squirmer::getStaticTypeID()
+{
+ return staticTypeID_;
+}
+
+} // namespace pe
+} // namespace walberla
diff --git a/src/pe/rigidbody/SquirmerFactory.cpp b/src/pe/rigidbody/SquirmerFactory.cpp
new file mode 100644
index 000000000..df95bd1c5
--- /dev/null
+++ b/src/pe/rigidbody/SquirmerFactory.cpp
@@ -0,0 +1,84 @@
+//======================================================================================================================
+//
+// 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 SquirmerFactory.cpp
+//! \author Michael Kuron <mkuron@icp.uni-stuttgart.de>
+//
+//======================================================================================================================
+
+#include "SquirmerFactory.h"
+
+#include "pe/Materials.h"
+#include "pe/rigidbody/BodyStorage.h"
+#include "pe/rigidbody/Squirmer.h"
+
+namespace walberla {
+namespace pe {
+
+SquirmerID createSquirmer( BodyStorage& globalStorage, BlockStorage& blocks, BlockDataID storageID,
+ id_t uid, const Vec3& gpos, real_t radius,
+ real_t squirmerVelocity, real_t squirmerBeta,
+ MaterialID material,
+ bool global, bool communicating, bool infiniteMass )
+{
+ WALBERLA_ASSERT_UNEQUAL( Squirmer::getStaticTypeID(), std::numeric_limits<id_t>::max(), "Squirmer TypeID not initalized!");
+ // Checking the radius
+ if( radius <= real_c(0) )
+ throw std::invalid_argument( "Invalid squirmer radius" );
+
+ SquirmerID squirmer = NULL;
+
+ if (global)
+ {
+ const id_t sid = UniqueID<RigidBody>::createGlobal();
+ WALBERLA_ASSERT_EQUAL(communicating, false);
+ WALBERLA_ASSERT_EQUAL(infiniteMass, true);
+ squirmer = new Squirmer(sid, uid, gpos, Vec3(0,0,0), Quat(), radius, squirmerVelocity, squirmerBeta, material, global, false, true);
+ globalStorage.add(squirmer);
+ } else
+ {
+ for (auto it = blocks.begin(); it != blocks.end(); ++it){
+ IBlock* block = (&(*it));
+ if (block->getAABB().contains(gpos))
+ {
+ const id_t sid( UniqueID<RigidBody>::create() );
+
+ Storage* bs = block->getData<Storage>(storageID);
+ squirmer = new Squirmer(sid, uid, gpos, Vec3(0,0,0), Quat(), radius, squirmerVelocity, squirmerBeta, material, global, communicating, infiniteMass);
+ squirmer->MPITrait.setOwner(Owner(MPIManager::instance()->rank(), block->getId().getID()));
+ (*bs)[0].add(squirmer);
+ }
+ }
+ }
+
+ if (squirmer != NULL)
+ {
+ // Logging the successful creation of the squirmer
+ WALBERLA_LOG_DETAIL(
+ "Created squirmer " << squirmer->getSystemID() << "\n"
+ << " User-ID = " << uid << "\n"
+ << " Global position = " << gpos << "\n"
+ << " Radius = " << radius << "\n"
+ << " LinVel = " << sphere->getLinearVel() << "\n"
+ << " Material = " << Material::getName( material )
+ );
+ }
+
+ return squirmer;
+}
+
+} // namespace pe
+} // namespace walberla
+
diff --git a/src/pe/rigidbody/SquirmerFactory.h b/src/pe/rigidbody/SquirmerFactory.h
new file mode 100644
index 000000000..01aca8784
--- /dev/null
+++ b/src/pe/rigidbody/SquirmerFactory.h
@@ -0,0 +1,74 @@
+//======================================================================================================================
+//
+// 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 SquirmerFactory.h
+//! \author Michael Kuron <mkuron@icp.uni-stuttgart.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "pe/Materials.h"
+#include "pe/Types.h"
+
+#include "blockforest/BlockForest.h"
+#include "core/debug/Debug.h"
+
+namespace walberla {
+namespace pe {
+
+//=================================================================================================
+//
+// SQUIRMER SETUP FUNCTIONS
+//
+//=================================================================================================
+
+//*************************************************************************************************
+/**
+ * \ingroup pe
+ * \brief Setup of a new Squirmer.
+ *
+ * \param globalStorage process local global storage
+ * \param blocks storage of all the blocks on this process
+ * \param storageID BlockDataID of the BlockStorage block datum
+ * \param uid The user-specific ID of the sphere.
+ * \param gpos The global position of the center of the sphere.
+ * \param radius The radius of the sphere \f$ (0..\infty) \f$.
+ * \param squirmerVelocity The velocity of the squirmer.
+ * \param squirmerBeta The dipolar characteristic of the squirmer.
+ * \param material The material of the sphere.
+ * \param global specifies if the sphere should be created in the global storage
+ * \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
+ * \return Handle for the new sphere.
+ * \exception std::invalid_argument Invalid sphere radius.
+ * \exception std::invalid_argument Invalid global sphere position.
+ *
+ * This function creates a squirmer primitive in the \b pe simulation system. The squirmer with
+ * user-specific ID \a uid is placed at the global position \a gpos, has the radius \a radius,
+ * and consists of the material \a material. Its self-propulsion velocity is \a squirmerVelocity
+ * and its dipolar characteristic is \a squirmerBeta.
+ *
+ */
+SquirmerID createSquirmer( BodyStorage& globalStorage, BlockStorage& blocks, BlockDataID storageID,
+ id_t uid, const Vec3& gpos, real_t radius,
+ real_t squirmerVelocity, real_t squirmerBeta,
+ MaterialID material = Material::find("iron"),
+ bool global = false, bool communicating = true, bool infiniteMass = false );
+//*************************************************************************************************
+
+} // namespace pe
+} // namespace walberla
+
diff --git a/tests/pe/Marshalling.cpp b/tests/pe/Marshalling.cpp
index aedbc3566..4139e331c 100644
--- a/tests/pe/Marshalling.cpp
+++ b/tests/pe/Marshalling.cpp
@@ -23,6 +23,7 @@
#include "pe/rigidbody/Box.h"
#include "pe/rigidbody/Capsule.h"
#include "pe/rigidbody/Sphere.h"
+#include "pe/rigidbody/Squirmer.h"
#include "pe/rigidbody/UnionFactory.h"
#include "pe/rigidbody/Union.h"
#include "pe/communication/DynamicMarshalling.h"
@@ -39,7 +40,7 @@ typedef boost::tuple<Sphere> UnionTypeTuple;
typedef Union< UnionTypeTuple > UnionT;
typedef UnionT* UnionID;
-typedef boost::tuple<Box, Capsule, Sphere, UnionT> BodyTuple ;
+typedef boost::tuple<Box, Capsule, Sphere, Squirmer, UnionT> BodyTuple ;
void testBox()
{
@@ -108,6 +109,24 @@ void testSphere()
WALBERLA_CHECK_EQUAL(s1.getSystemID(), s2->getSystemID());
}
+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);
+
+ mpi::SendBuffer sb;
+ WALBERLA_ASSERT_UNEQUAL(Sphere::getStaticTypeID(), Squirmer::getStaticTypeID(), "Squirmer did not get its own type ID");
+ WALBERLA_ASSERT_UNEQUAL(Sphere::getStaticTypeID(), s1.getStaticTypeID(), "Squirmer did not get its own type ID");
+ MarshalDynamically<BodyTuple>::execute(sb, s1);
+ mpi::RecvBuffer rb(sb);
+
+ SquirmerID s2 = static_cast<SquirmerID> (UnmarshalDynamically<BodyTuple>::execute(rb, Squirmer::getStaticTypeID(), math::AABB(Vec3(-100,-100,-100), Vec3(100,100,100)), math::AABB(Vec3(-100,-100,-100), Vec3(100,100,100))));
+
+ WALBERLA_CHECK_FLOAT_EQUAL(s1.getSquirmerVelocity(), s2->getSquirmerVelocity());
+ WALBERLA_CHECK_FLOAT_EQUAL(s1.getSquirmerBeta(), s2->getSquirmerBeta());
+}
+
void testUnion()
{
UnionT u1(159, 423, Vec3(real_c(1), real_c(2), real_c(3)), Vec3(0,0,0), Quat(), false, false, false);
@@ -158,6 +177,7 @@ int main( int argc, char** argv )
testBox();
testCapsule();
testUnion();
+ testSquirmer();
return EXIT_SUCCESS;
}
--
GitLab