diff --git a/src/pe/Types.h b/src/pe/Types.h
index 15d386d540f0854f7db0867611a611ea976f35a2..d37661d54fb2d76466a4ae9b3f4cc8fe608456c2 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/rigidbody/Squirmer.cpp b/src/pe/communication/rigidbody/Squirmer.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d4b790b99508f884f8f44e8982005a20d826bdfe
--- /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 0000000000000000000000000000000000000000..9029c93586fa36c5d3b4e29d2aaefe74a4d5a4d1
--- /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 c67cac62c77b044ed564cf86767229d862b80f89..45d40109e31a29020a55517e8c28b1b5fecddf2f 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 408750a521d7d321a4be45303b5114dee9cdd56e..c58df20ee46b0b2674b8e6a9be542ad813ac3e56 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 78ceb34de931b4e3df2b770d489e88c302dbd4fb..8ea3a81a768d586f47ab77b3016cead0e6b85fb3 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 edfa5acfb8fd94e7510125d1ce3a86b52912689e..bc0e82566c5f171eee4aa354a2dbedb7e75033b7 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 0000000000000000000000000000000000000000..99f650b00290b619fa411c1a74f77a8c193c6200
--- /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( 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 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.
+ *
+ * \param rpos The relative global coordinate.
+ * \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
+{
+   return Sphere::velFromWF( gpos ) + getSquirmerVelocity( gpos - gpos_ );
+}
+//*************************************************************************************************
+
+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 0000000000000000000000000000000000000000..3dc76acad0ecd56c46f9d17f3027fbd3f4fd3d7e
--- /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 0000000000000000000000000000000000000000..df95bd1c509ed05786b7c25bebd7b69190ccc1b1
--- /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 0000000000000000000000000000000000000000..01aca87847b4571e08f2e17db0919a98dd001e8b
--- /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/src/pe_coupling/momentum_exchange_method/boundary/SimpleBB.h b/src/pe_coupling/momentum_exchange_method/boundary/SimpleBB.h
index a6807a6ff0a78df2cf23f7831a298821cd1bb0eb..161dfc872e1094df917d2394048790afdf52da92 100644
--- a/src/pe_coupling/momentum_exchange_method/boundary/SimpleBB.h
+++ b/src/pe_coupling/momentum_exchange_method/boundary/SimpleBB.h
@@ -59,8 +59,6 @@ namespace pe_coupling {
 template< typename LatticeModel_T, typename FlagField_T >
 class SimpleBB : public Boundary< typename FlagField_T::flag_t >
 {
-   static_assert( LatticeModel_T::compressible == false, "Only works with incompressible models!" );
-
    typedef lbm::PdfField< LatticeModel_T >   PDFField_T;
    typedef typename LatticeModel_T::Stencil  Stencil_T;
    typedef typename FlagField_T::flag_t      flag_t;
@@ -172,9 +170,21 @@ inline void SimpleBB< LatticeModel_T, FlagField_T >::treatDirection( const cell_
    auto boundaryVelocity = body.velFromWF( boundaryPosition );
 
    // include effect of boundary velocity
-   pdf_new -= real_c(3.0) * alpha * LatticeModel_T::w[ Stencil_T::idx[dir] ] * ( real_c( stencil::cx[ dir ] ) * boundaryVelocity[0] +
-                                                                               real_c( stencil::cy[ dir ] ) * boundaryVelocity[1] +
-                                                                               real_c( stencil::cz[ dir ] ) * boundaryVelocity[2] );
+   if( LatticeModel_T::compressible )
+   {
+       const auto density  = pdfField_->getDensity(x,y,z);
+       pdf_new -= real_c(3.0) * alpha * density * LatticeModel_T::w[ Stencil_T::idx[dir] ] *
+                     ( real_c( stencil::cx[ dir ] ) * boundaryVelocity[0] +
+                       real_c( stencil::cy[ dir ] ) * boundaryVelocity[1] +
+                       real_c( stencil::cz[ dir ] ) * boundaryVelocity[2] );
+   }
+   else
+   {
+       pdf_new -= real_c(3.0) * alpha * LatticeModel_T::w[ Stencil_T::idx[dir] ] *
+                     ( real_c( stencil::cx[ dir ] ) * boundaryVelocity[0] +
+                       real_c( stencil::cy[ dir ] ) * boundaryVelocity[1] +
+                       real_c( stencil::cz[ dir ] ) * boundaryVelocity[2] );
+   }
 
    // carry out the boundary handling
    pdfField_->get( nx, ny, nz, Stencil_T::invDirIdx(dir) ) = pdf_new;
diff --git a/tests/pe/Marshalling.cpp b/tests/pe/Marshalling.cpp
index aedbc3566f44d99ccc0cc5e7c8fec198b3211c90..57411b0d279cf2e5524170450d82688314fc6697 100644
--- a/tests/pe/Marshalling.cpp
+++ b/tests/pe/Marshalling.cpp
@@ -23,8 +23,10 @@
 #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/rigidbody/Squirmer.h"
 #include "pe/communication/DynamicMarshalling.h"
 #include "pe/rigidbody/SetBodyTypeIDs.h"
 #include "pe/Materials.h"
@@ -39,7 +41,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 +110,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 +178,7 @@ int main( int argc, char** argv )
    testBox();
    testCapsule();
    testUnion();
+   testSquirmer();
 
    return EXIT_SUCCESS;
 }
diff --git a/tests/pe_coupling/CMakeLists.txt b/tests/pe_coupling/CMakeLists.txt
index 24b8d2a53791b812d2aa51b363e60408357103dc..e31ded9e826492cb2926ada9d7272dd4647bccd6 100644
--- a/tests/pe_coupling/CMakeLists.txt
+++ b/tests/pe_coupling/CMakeLists.txt
@@ -123,3 +123,6 @@ waLBerla_execute_test( NAME TorqueSphereMEMMRFuncTest      COMMAND $<TARGET_FILE
 waLBerla_execute_test( NAME TorqueSphereMEMMRSingleTest    COMMAND $<TARGET_FILE:TorqueSphereMEM> --MO_MR              PROCESSES 1 LABELS longrun     CONFIGURATIONS Release RelWithDbgInfo )
 waLBerla_execute_test( NAME TorqueSphereMEMMRParallelTest  COMMAND $<TARGET_FILE:TorqueSphereMEM> --MO_MR              PROCESSES 8 LABELS verylongrun CONFIGURATIONS Release RelWithDbgInfo )  
 
+waLBerla_compile_test( FILES momentum_exchange_method/SquirmerTest.cpp DEPENDS blockforest pe timeloop )
+waLBerla_execute_test( NAME SquirmerShortTest COMMAND $<TARGET_FILE:SquirmerTest> --shortrun PROCESSES 1 )
+waLBerla_execute_test( NAME SquirmerTest COMMAND $<TARGET_FILE:SquirmerTest> PROCESSES 1 CONFIGURATIONS Release RelWithDbgInfo  )
diff --git a/tests/pe_coupling/momentum_exchange_method/SquirmerTest.cpp b/tests/pe_coupling/momentum_exchange_method/SquirmerTest.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..235fcc2a62739f4d5cb559e164a34cf143db82a3
--- /dev/null
+++ b/tests/pe_coupling/momentum_exchange_method/SquirmerTest.cpp
@@ -0,0 +1,474 @@
+//======================================================================================================================
+//
+//  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 SquirmerTest.cpp
+//! \ingroup pe_coupling
+//! \author Christian Burkard <buch@icp.uni-stuttgart.de>
+//! \author Michael Kuron <mkuron@icp.uni-stuttgart.de>
+//
+//======================================================================================================================
+
+#include "blockforest/Initialization.h"
+#include "blockforest/communication/UniformBufferedScheme.h"
+
+#include "boundary/all.h"
+
+#include "core/DataTypes.h"
+#include "core/Environment.h"
+#include "core/debug/Debug.h"
+#include "core/debug/TestSubsystem.h"
+#include "core/logging/Logging.h"
+#include "core/math/all.h"
+#include "core/SharedFunctor.h"
+#include "core/timing/RemainingTimeLogger.h"
+#include "core/mpi/MPIManager.h"
+
+#include "domain_decomposition/SharedSweep.h"
+
+#include "field/AddToStorage.h"
+#include "field/communication/PackInfo.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/lattice_model/ForceModel.h"
+#include "lbm/sweeps/CellwiseSweep.h"
+#include "lbm/sweeps/SweepWrappers.h"
+
+#include "stencil/D3Q7.h"
+
+#include "timeloop/SweepTimeloop.h"
+
+#include "pe/rigidbody/SquirmerFactory.h"
+#include "pe/communication/rigidbody/Squirmer.h"
+#include "pe/basic.h"
+
+#include "pe_coupling/mapping/all.h"
+#include "pe_coupling/momentum_exchange_method/all.h"
+#include "pe_coupling/utility/all.h"
+
+#include "field/vtk/FlagFieldCellFilter.h"
+#include "field/vtk/VTKWriter.h"
+#include "lbm/vtk/Velocity.h"
+#include "vtk/VTKOutput.h"
+
+///////////
+// USING //
+///////////
+
+using namespace walberla;
+using walberla::uint_t;
+
+// PDF field, flag field & body field
+typedef lbm::force_model::None ForceModel_T;
+typedef lbm::D3Q19< lbm::collision_model::TRT, false, ForceModel_T >  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 >  PeBodyField_T;
+
+const uint_t FieldGhostLayers = 1;
+
+// boundary handling
+typedef pe_coupling::SimpleBB       < LatticeModel_T, FlagField_T >  MO_BB_T;
+
+typedef boost::tuples::tuple< MO_BB_T >               BoundaryConditions_T;
+typedef BoundaryHandling< FlagField_T, Stencil_T, BoundaryConditions_T > BoundaryHandling_T;
+
+typedef boost::tuple<pe::Squirmer> BodyTypeTuple ;
+
+///////////
+// FLAGS //
+///////////
+
+const FlagUID Fluid_Flag   ( "fluid" );
+const FlagUID MO_BB_Flag   ( "moving obstacle BB" );
+const FlagUID FormerMO_BB_Flag   ( "former moving obstacle BB" );
+
+/////////////////////////////////////
+// BOUNDARY HANDLING CUSTOMIZATION //
+/////////////////////////////////////
+
+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
+
+////////////////////////////////
+// VTK OUTPUT HELPER FUNCTION //
+////////////////////////////////
+
+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_ );
+   PeBodyField_T * bodyField     = block->getData< PeBodyField_T >( bodyFieldID_ );
+
+   const auto fluid = flagField->flagExists( Fluid_Flag ) ? flagField->getFlag( Fluid_Flag ) : flagField->registerFlag( Fluid_Flag );
+
+   BoundaryHandling_T * handling = new BoundaryHandling_T( "fixed obstacle boundary handling", flagField, fluid,
+                                                           boost::tuples::make_tuple( MO_BB_T (  "MO_BB",  MO_BB_Flag, pdfField, flagField, bodyField, fluid, *storage, *block ) ) );
+
+   handling->fillWithDomain( FieldGhostLayers );
+
+   return handling;
+}
+
+shared_ptr< vtk::VTKOutput> createFluidFieldVTKWriter( shared_ptr< StructuredBlockForest > & blocks,
+                                                       const BlockDataID & pdfFieldId, const BlockDataID & flagFieldId,
+                                                       std::string fname )
+{
+   auto pdfFieldVTKWriter = vtk::createVTKOutput_BlockData( blocks, fname, uint_c(1), uint_c(0) );
+
+   field::FlagFieldCellFilter< FlagField_T > fluidFilter( flagFieldId );
+   fluidFilter.addFlag( Fluid_Flag );
+   pdfFieldVTKWriter->addCellInclusionFilter( fluidFilter );
+
+   auto velocityWriter = make_shared< lbm::VelocityVTKWriter< LatticeModel_T, float > >( pdfFieldId, "VelocityFromPDF" );
+   pdfFieldVTKWriter->addCellDataWriter( velocityWriter );
+
+   return pdfFieldVTKWriter;
+}
+
+class ResetPosition
+{
+public:
+   
+   explicit ResetPosition( const shared_ptr<StructuredBlockStorage> & blockStorage,
+                           const BlockDataID & bodyStorageID,
+                           real_t dt )
+   : blockStorage_( blockStorage )
+   , bodyStorageID_( bodyStorageID )
+   , dt_( dt )
+   {}
+   
+   void operator()()
+   {
+      for (auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt)
+      {
+         for( auto bodyIt = pe::BodyIterator::begin(*blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+         {
+            bodyIt->setPosition( bodyIt->getPosition() - bodyIt->getLinearVel() * dt_ );
+         }
+      }
+   }
+   
+private:
+   
+   shared_ptr<StructuredBlockStorage> blockStorage_;
+   const BlockDataID &  bodyStorageID_;
+   
+   real_t dt_;
+   
+}; // class ResetPosition
+
+
+Vector3<real_t> analytical_velocity(Vector3<real_t> r, real_t b1, real_t b2, Vector3<real_t> e, real_t radius)
+{
+   auto r_len = r.length();
+   auto rs = r.getNormalized();
+   auto frac = radius/r_len;
+   auto frac2 = frac*frac;
+   auto frac3 = frac2*frac;
+   auto frac4 = frac3*frac;
+
+   return b1 * frac3 * ((e*rs) * rs - e / real_c(3.0)) + ( frac4 - frac2 ) * b2 * 0.5 * ( 3 * (e * rs) * (e * rs) - 1 ) * rs + frac4 * b2 * (e*rs) * ( (e*rs) * rs - e);
+}
+
+
+int main( int argc, char **argv )
+{
+   debug::enterTestMode();
+
+   mpi::Environment env( argc, argv );
+
+   auto processes = MPIManager::instance()->numProcesses();
+   
+   if( processes != 1 && processes != 2 && processes != 4 && processes != 8)
+   {
+      std::cerr << "Number of processes must be equal to either 1, 2, 4, or 8!" << std::endl;
+      return EXIT_FAILURE;
+   }
+   
+   bool shortrun = false;
+   bool writevtk = false;
+   for( int i = 1; i < argc; ++i )
+   {
+      if( std::strcmp( argv[i], "--shortrun" ) == 0 ) { shortrun       = true; continue; }
+      if( std::strcmp( argv[i], "--vtk" )      == 0 ) { writevtk       = true; continue; }
+   }
+   
+   ///////////////////////////
+   // SIMULATION PROPERTIES //
+   ///////////////////////////
+   
+   real_t R = 3;
+   uint_t L = 24;
+   uint_t T = 1000;
+   if (shortrun)
+   {
+       L = 10;
+       T = 10;
+   }
+   const real_t dx = real_c(1.0);
+   const real_t dt = real_c(1.0);
+   const real_t omega = lbm::collision_model::omegaFromViscosity( real_c(1.0/6.0) );
+   const real_t v0 = real_c(0.01);
+   const real_t beta = real_c(5.0);
+   
+   ///////////////////////////
+   // BLOCK STRUCTURE SETUP //
+   ///////////////////////////
+   const uint_t XBlocks = (processes >= 2) ? uint_t( 2 ) : uint_t( 1 );
+   const uint_t YBlocks = (processes >= 4) ? uint_t( 2 ) : uint_t( 1 );
+   const uint_t ZBlocks = (processes == 8) ? uint_t( 2 ) : uint_t( 1 );
+   const uint_t XCells = L / XBlocks;
+   const uint_t YCells = L / YBlocks;
+   const uint_t ZCells = L / ZBlocks;
+
+   // create fully periodic domain
+   auto blocks = blockforest::createUniformBlockGrid( XBlocks, YBlocks, ZBlocks, XCells, YCells, ZCells, dx, true,
+                                                      true, true, true );
+   
+   ////////
+   // PE //
+   ////////
+
+   auto globalBodyStorage = make_shared<pe::BodyStorage>();
+   pe::SetBodyTypeIDs<BodyTypeTuple>::execute();
+   auto bodyStorageID = blocks->addBlockData(pe::createStorageDataHandling<BodyTypeTuple>(), "pe Body Storage");
+   auto ccdID         = blocks->addBlockData(pe::ccd::createHashGridsDataHandling( globalBodyStorage, bodyStorageID ), "CCD");
+   auto fcdID         = blocks->addBlockData(pe::fcd::createSimpleFCDDataHandling<BodyTypeTuple>(), "FCD");
+   
+   pe::cr::HCSITS cr(globalBodyStorage, blocks->getBlockStoragePointer(), bodyStorageID, ccdID, fcdID);
+   
+   /////////////////
+   // PE COUPLING //
+   /////////////////
+
+   // connect to pe
+   const real_t overlap = real_c( 1.5 ) * dx;
+
+   if( R > real_c( L ) * real_c(0.5) - overlap )
+   {
+      std::cerr << "Periodic sphere is too large and would lead to invalid PE state!" << std::endl;
+      return EXIT_FAILURE;
+   }
+   
+   boost::function<void(void)> syncCall = boost::bind( pe::syncShadowOwners<BodyTypeTuple>, boost::ref(blocks->getBlockForest()), bodyStorageID, static_cast<WcTimingTree*>(NULL), overlap, false );
+   
+   const auto myMat = pe::createMaterial( "myMat", real_c(1), real_t(0), real_t(1), real_t(1), real_t(0), real_t(1), real_t(1), real_t(0), real_t(0) );
+
+   // create the squirmer in the middle of the domain
+   const Vector3<real_t> position (real_c(L)*real_c(0.5), real_c(L)*real_c(0.5), real_c(L)*real_c(0.5));
+   const Vector3<real_t> up(0.0, 0.0, 1.0);
+   const Vector3<real_t> orientation(1.0, 0.0, 0.0);
+   const auto w = std::acos(up*orientation);
+   math::Quaternion<real_t> q(up % orientation, w);
+   auto squirmer = pe::createSquirmer(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0, position, R, v0, beta, myMat);
+   if (squirmer)
+      squirmer->setOrientation(q);
+
+   syncCall();
+
+   ///////////////////////
+   // ADD DATA TO BLOCKS //
+   ////////////////////////
+
+   // create the lattice model
+   LatticeModel_T latticeModel = LatticeModel_T( lbm::collision_model::TRT::constructWithMagicNumber( omega ) );
+
+   // add PDF field ( uInit = <0,0,0>, rhoInit = 1 )
+   BlockDataID pdfFieldID = lbm::addPdfFieldToStorage< LatticeModel_T >( blocks, "pdf field (zyxf)", latticeModel, FieldGhostLayers );
+
+   // add flag field
+   BlockDataID flagFieldID = field::addFlagFieldToStorage< FlagField_T >( blocks, "flag field" );
+
+   // add body field
+   BlockDataID bodyFieldID = field::addToStorage< PeBodyField_T >( blocks, "body field", NULL, field::zyxf );
+
+   // add boundary handling & initialize outer domain boundaries
+   BlockDataID boundaryHandlingID = blocks->addStructuredBlockData< BoundaryHandling_T >(
+                                    MyBoundaryHandling( flagFieldID, pdfFieldID, bodyFieldID ), "boundary handling" );
+
+   ///////////////
+   // TIME LOOP //
+   ///////////////
+
+   // create the timeloop
+   SweepTimeloop timeloop( blocks->getBlockStorage(), T );
+   
+   // sweep for updating the pe body mapping into the LBM simulation
+   timeloop.add()
+   << Sweep( pe_coupling::BodyMapping< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MO_BB_Flag, FormerMO_BB_Flag ),
+            "Body Mapping" );
+   
+   // sweep for restoring PDFs in cells previously occupied by pe bodies
+   typedef pe_coupling::EquilibriumReconstructor< LatticeModel_T, BoundaryHandling_T > Reconstructor_T;
+   Reconstructor_T reconstructor( blocks, boundaryHandlingID, pdfFieldID, bodyFieldID );
+   timeloop.add()
+   << Sweep( pe_coupling::PDFReconstruction< LatticeModel_T, BoundaryHandling_T, Reconstructor_T >( blocks, boundaryHandlingID, bodyStorageID, bodyFieldID,
+                                                                                                   reconstructor, FormerMO_BB_Flag, Fluid_Flag  ), "PDF Restore" );
+
+   // setup of the LBM communication for synchronizing the pdf field between neighboring blocks
+   boost::function< void () > commFunction;
+
+   blockforest::communication::UniformBufferedScheme< stencil::D3Q27 > scheme( blocks );
+   scheme.addPackInfo( make_shared< field::communication::PackInfo< PdfField_T > >( pdfFieldID ) );
+   commFunction = scheme;
+
+   // uses standard bounce back boundary conditions
+   pe_coupling::mapMovingBodies< BoundaryHandling_T >( blocks, boundaryHandlingID, bodyStorageID, bodyFieldID, MO_BB_Flag );
+   
+   auto sweep = lbm::makeCellwiseSweep< LatticeModel_T, FlagField_T >( pdfFieldID, flagFieldID, Fluid_Flag );
+   
+   // collision sweep
+   timeloop.add() << Sweep( lbm::makeCollideSweep( sweep ), "cell-wise LB sweep (collide)" );
+
+   // add LBM communication function and boundary handling sweep (does the hydro force calculations and the no-slip treatment)
+   timeloop.add() << BeforeFunction( commFunction, "LBM Communication" )
+                  << Sweep( BoundaryHandling_T::getBlockSweep( boundaryHandlingID ), "Boundary Handling" );
+
+   // streaming
+   timeloop.add() << Sweep( lbm::makeStreamSweep( sweep ), "cell-wise LB sweep (stream)" );
+   
+   // add pe timesteps
+   timeloop.addFuncAfterTimeStep( pe_coupling::TimeStep( blocks, bodyStorageID, cr, syncCall, dt ), "pe Time Step" );
+   timeloop.addFuncAfterTimeStep( ResetPosition( blocks, bodyStorageID, dt ), "Reset pe positions" );
+
+   timeloop.addFuncAfterTimeStep( RemainingTimeLogger( timeloop.getNrOfTimeSteps() ), "Remaining Time Logger" );
+
+   ////////////////////////
+   // EXECUTE SIMULATION //
+   ////////////////////////
+   
+   WALBERLA_LOG_INFO_ON_ROOT("Starting test...");
+   WcTimingPool timeloopTiming;
+   timeloop.run( timeloopTiming );
+   timeloopTiming.logResultOnRoot();
+
+   if (writevtk)
+   {
+      WALBERLA_LOG_INFO_ON_ROOT("Writing out VTK file...");
+      auto pdfFieldVTKWriter = createFluidFieldVTKWriter( blocks, pdfFieldID, flagFieldID, "fluid_field" );
+      auto wf = vtk::writeFiles( pdfFieldVTKWriter );
+      wf();
+      auto flagFieldVTK = vtk::createVTKOutput_BlockData( blocks, "flag_field", 1, 0 );
+      flagFieldVTK->addCellDataWriter( make_shared< field::VTKWriter< FlagField_T > >( flagFieldID, "FlagField" ) );
+      vtk::writeFiles( flagFieldVTK )();
+   }
+   
+   if (squirmer) // this object only exists on the node that contains position
+   {
+      WALBERLA_ASSERT_FLOAT_EQUAL(squirmer->getSquirmerVelocity(), v0);
+      WALBERLA_ASSERT_FLOAT_EQUAL(squirmer->getSquirmerBeta(), beta);
+      WALBERLA_ASSERT_FLOAT_EQUAL(squirmer->getRadius(), R);
+      WALBERLA_ASSERT_FLOAT_EQUAL(squirmer->getPosition(), position);
+      WALBERLA_ASSERT_FLOAT_EQUAL(squirmer->getQuaternion(), q);
+   }
+   
+   WALBERLA_LOG_INFO_ON_ROOT("Checking result...");
+   auto b1 = real_c(1.5) * v0;
+   auto b2 = beta * b1;
+   auto e = q.rotate(up).getNormalized();
+   auto radius = R;
+   auto squirmer_pos = position;
+
+   real_t abs_tolerance = real_c(0.0026);
+   real_t rel_tolerance = real_c(0.10);
+
+   for ( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt)
+   {
+      IBlock & block = *blockIt;
+      auto src = block.getData< PdfField_T >(pdfFieldID);
+
+      auto flagField = block.getData<FlagField_T> (flagFieldID);
+      auto flag = flagField->getFlag(Fluid_Flag);
+
+      WALBERLA_FOR_ALL_CELLS_XYZ(src, {
+         Vector3<real_t> pos;
+         blocks->getBlockLocalCellCenter( block, Cell(x,y,z), pos[0], pos[1], pos[2] );
+
+         if( flagField->isFlagSet(x, y, z, flag) ) //check for fluid/non-boundary
+         {
+            auto r = pos - squirmer_pos;
+            
+            auto v_ana = analytical_velocity(r, b1, b2, e, radius);
+            
+            // Superposition with the nearest neighbours (faces)
+            for( auto s = stencil::D3Q7::beginNoCenter(); s != stencil::D3Q7::end(); ++s )
+            {
+               Vector3<real_t> imidpoint(squirmer_pos);
+               imidpoint[0] += real_c(s.cx() * int_c(L));
+               imidpoint[1] += real_c(s.cy() * int_c(L));
+               imidpoint[2] += real_c(s.cz() * int_c(L));
+               v_ana += analytical_velocity( pos - imidpoint, b1, b2, e, radius );
+            }
+            
+            if( !shortrun && r.length() > 2*radius ) // Exclude the volume around the squirmer, as the analytical solution is a far field solution only
+            {
+               auto v_data = src->getVelocity(x, y, z);
+               auto diff = v_data - v_ana;
+               
+               for( uint_t d = 0; d < 3; ++d )
+               {
+                  if( std::abs( diff[d] / v_ana[d] ) > rel_tolerance && std::abs( diff[d] ) > abs_tolerance )
+                  {
+                     WALBERLA_LOG_DEVEL("Difference too large in " << Cell(x,y,z) << " (" << r.length() << "). Expected " << v_ana << ", got " << v_data << ", relative " << std::abs(diff[d] / v_ana[d]) << ", absolute " << std::abs(diff[d]));
+                  }
+               }
+            }
+            
+            if( writevtk )
+            {
+               // store the value into the PDF field so we can write it to VTK later
+               src->setToEquilibrium( x, y, z, v_ana );
+            }
+         }
+      });
+   }
+   
+   if (writevtk)
+   {
+      WALBERLA_LOG_INFO_ON_ROOT("Writing out reference VTK file...");
+      auto pdfFieldVTKWriter = createFluidFieldVTKWriter( blocks, pdfFieldID, flagFieldID, "fluid_field_ref" );
+      auto wf = vtk::writeFiles( pdfFieldVTKWriter );
+      wf();
+   }
+   
+   WALBERLA_LOG_INFO_ON_ROOT("Completed test.");
+   
+   return 0;
+}