From c45f3a3efebad801b9a7e9bab513814d71e42b93 Mon Sep 17 00:00:00 2001
From: Sebastian Eibl <sebastian.eibl@fau.de>
Date: Tue, 12 Jun 2018 13:55:51 +0200
Subject: [PATCH] [API] rebased Union on BodyStorage

---
 src/pe/communication/rigidbody/Union.h  |  12 +-
 src/pe/fcd/AnalyticCollisionDetection.h |   4 +-
 src/pe/fcd/GJKEPACollideFunctor.h       |   6 +-
 src/pe/rigidbody/Union.h                | 226 +++++++++---------------
 src/pe/rigidbody/UnionFactory.h         |  20 +--
 tests/pe/Collision.cpp                  |   8 +-
 tests/pe/CollisionTobiasGJK.cpp         |  32 ++--
 tests/pe/Marshalling.cpp                |   4 +-
 tests/pe/ShadowCopy.cpp                 |  16 +-
 tests/pe/Union.cpp                      |  21 ++-
 10 files changed, 143 insertions(+), 206 deletions(-)

diff --git a/src/pe/communication/rigidbody/Union.h b/src/pe/communication/rigidbody/Union.h
index 69c2d6bb1..596f54118 100644
--- a/src/pe/communication/rigidbody/Union.h
+++ b/src/pe/communication/rigidbody/Union.h
@@ -72,12 +72,12 @@ void marshal( mpi::SendBuffer& buffer, const Union<BodyTypeTuple>& obj )
    buffer << static_cast<size_t> (obj.size());                  // Encoding the number of contained bodies
 
    // Encoding the contained primitives
-   const typename Union<BodyTypeTuple>::ConstIterator begin( obj.begin() );
-   const typename Union<BodyTypeTuple>::ConstIterator end  ( obj.end()   );
-   for(  typename Union<BodyTypeTuple>::ConstIterator body=begin; body!=end; ++body )
+   const typename Union<BodyTypeTuple>::const_iterator begin( obj.begin() );
+   const typename Union<BodyTypeTuple>::const_iterator end  ( obj.end()   );
+   for(  typename Union<BodyTypeTuple>::const_iterator body=begin; body!=end; ++body )
    {
       buffer << body->getTypeID();
-      MarshalDynamically<BodyTypeTuple>::execute( buffer, **body );
+      MarshalDynamically<BodyTypeTuple>::execute( buffer, *body );
    }
 }
 
@@ -132,9 +132,9 @@ inline Union<BodyTypeTuple>* instantiate( mpi::RecvBuffer& buffer, const math::A
    {
       decltype ( static_cast<BodyID>(nullptr)->getTypeID() ) type;
       buffer >> type;
-      BodyID obj = UnmarshalDynamically<BodyTypeTuple>::execute(buffer, type, domain, block);
+      BodyPtr obj( UnmarshalDynamically<BodyTypeTuple>::execute(buffer, type, domain, block) );
       obj->setRemote( true );
-      newBody->add(obj);
+      newBody->add(std::move(obj));
    }
 
    return newBody;
diff --git a/src/pe/fcd/AnalyticCollisionDetection.h b/src/pe/fcd/AnalyticCollisionDetection.h
index 8b6fa2853..5318eaae9 100644
--- a/src/pe/fcd/AnalyticCollisionDetection.h
+++ b/src/pe/fcd/AnalyticCollisionDetection.h
@@ -2124,7 +2124,7 @@ bool collide( Union<BodyTypeTuple>* bd1, BodyB* bd2, Container& container )
    bool collision = false;
    for( auto it=bd1->begin(); it!=bd1->end(); ++it )
    {
-      collision |= SingleCast<BodyTypeTuple, AnalyticSingleCollideFunctor<BodyB, Container>, bool>::execute(*it, func);
+      collision |= SingleCast<BodyTypeTuple, AnalyticSingleCollideFunctor<BodyB, Container>, bool>::execute(it.getBodyID(), func);
    }
    return collision;
 }
@@ -2146,7 +2146,7 @@ bool collide( Union<BodyTypeTupleA>* bd1, Union<BodyTypeTupleB>* bd2, Container&
    {
       for( auto it2=bd2->begin(); it2!=bd2->end(); ++it2 )
       {
-         collision |= DoubleCast<BodyTypeTupleA, BodyTypeTupleB, AnalyticCollideFunctor<Container>, bool>::execute(*it1, *it2, func);
+         collision |= DoubleCast<BodyTypeTupleA, BodyTypeTupleB, AnalyticCollideFunctor<Container>, bool>::execute(it1.getBodyID(), it2.getBodyID(), func);
       }
    }
    return collision;
diff --git a/src/pe/fcd/GJKEPACollideFunctor.h b/src/pe/fcd/GJKEPACollideFunctor.h
index 35a72fcac..051d5bcbc 100644
--- a/src/pe/fcd/GJKEPACollideFunctor.h
+++ b/src/pe/fcd/GJKEPACollideFunctor.h
@@ -169,7 +169,7 @@ namespace gjkepa{
       bool collision = false;
       for( auto it=a->begin(); it!=a->end(); ++it )
       {
-         collision |= SingleCast<BodyTupleA, GJKEPASingleCollideFunctor<BodyB, Container>, bool>::execute(*it, func);
+         collision |= SingleCast<BodyTupleA, GJKEPASingleCollideFunctor<BodyB, Container>, bool>::execute(it.getBodyID(), func);
       }
       return collision;
    }
@@ -187,7 +187,7 @@ namespace gjkepa{
       {
          for( auto it2=b->begin(); it2!=b->end(); ++it2 )
          {
-            collision |= DoubleCast<BodyTupleA, BodyTupleB, GJKEPACollideFunctor<Container>, bool>::execute(*it1, *it2, func);
+            collision |= DoubleCast<BodyTupleA, BodyTupleB, GJKEPACollideFunctor<Container>, bool>::execute(it1.getBodyID(), it2.getBodyID(), func);
          }
       }
       return collision;
@@ -200,7 +200,7 @@ namespace gjkepa{
       bool collision = false;
       for( auto it=a->begin(); it!=a->end(); ++it )
       {
-         collision |= SingleCast<BodyTupleA, GJKEPASingleCollideFunctor<Plane, Container>, bool>::execute(*it, func);
+         collision |= SingleCast<BodyTupleA, GJKEPASingleCollideFunctor<Plane, Container>, bool>::execute(it.getBodyID(), func);
       }
       return collision;
    }
diff --git a/src/pe/rigidbody/Union.h b/src/pe/rigidbody/Union.h
index 81c8128f6..fe6e0ce94 100644
--- a/src/pe/rigidbody/Union.h
+++ b/src/pe/rigidbody/Union.h
@@ -27,7 +27,10 @@
 //*************************************************************************************************
 
 #include <pe/Config.h>
+#include <pe/rigidbody/BodyStorage.h>
 #include <pe/rigidbody/RigidBody.h>
+#include <pe/rigidbody/RigidBodyCastIterator.h>
+#include <pe/rigidbody/RigidBodyIterator.h>
 #include <pe/Types.h>
 
 #include <core/debug/Debug.h>
@@ -65,23 +68,24 @@ class Union : public RigidBody
 {
 public:
    //**Type definitions****************************************************************************
-   typedef BodyTypeTuple                     BodyTypeTupleT;
+   using BodyTypeTupleT        = BodyTypeTuple;
 
-   typedef PtrVector<RigidBody, PtrDelete>   Bodies;             //!< Iterator over the contained rigid bodies.
-   typedef Bodies::Iterator                  Iterator;           //!< Iterator over the contained rigid bodies.
-   typedef Bodies::ConstIterator             ConstIterator;      //!< Constant iterator over the contained rigid bodies.
-   //**++++++++++++++++****************************************************************************
+   //**********************************************************************************************
 
-   //**Forward declarations for nested classes*****************************************************
-//   template< typename C > class CastIterator;
-//   template< typename C > class ConstCastIterator;
+   using size_type             = BodyStorage::size_type;           //!< Size type of the body storage.
+   using iterator              = BodyStorage::iterator;            //!< Iterator over non-const bodies.
+   using const_iterator        = BodyStorage::const_iterator;      //!< Iterator over constant bodies.
+   template <typename C> //cast type
+   using cast_iterator         = BodyStorage::cast_iterator<C>;
+   template <typename C> //cast type
+   using const_cast_iterator   = BodyStorage::const_cast_iterator<C>;
    //**********************************************************************************************
 
    //**Constructors********************************************************************************
    /*!\name Constructors */
    //@{
    explicit Union( id_t sid, id_t uid, const Vec3& gpos, const Vec3& rpos, const Quat& q,
-                    const bool global, const bool communicating, const bool infiniteMass );
+                   const bool global, const bool communicating, const bool infiniteMass );
    //@}
    //**********************************************************************************************
 
@@ -95,28 +99,35 @@ public:
 
 public:
    //**Get functions*******************************************************************************
-   /*!\name Get functions */
+   /*!\name BodyStorage functions */
    //@{
-                          inline bool                 isEmpty() const;
-                          inline size_t               size()    const;
-   template< typename C > inline size_t               size()    const {return bodies_.size<C>();}
-                          inline BodyID               getBody( size_t index );
-                          inline ConstBodyID          getBody( size_t index ) const;
-
-                          inline Iterator                     begin()       {return bodies_.begin();}
-                          inline ConstIterator                begin() const {return bodies_.begin();}
-   template< typename C > inline Bodies::CastIterator<C>      begin()       {return bodies_.begin<C>();}
-   template< typename C > inline Bodies::ConstCastIterator<C> begin() const {return bodies_.begin<C>();}
-
-                          inline Iterator                     end()         {return bodies_.end();}
-                          inline ConstIterator                end() const   {return bodies_.end();}
-   template< typename C > inline Bodies::CastIterator<C>      end()         {return bodies_.end<C>();}
-   template< typename C > inline Bodies::ConstCastIterator<C> end() const   {return bodies_.end<C>();}
-
-   virtual inline real_t getVolume()         const;
+   inline bool                 isEmpty() const {return bodies_.isEmpty();}
+   inline size_t               size()    const {return bodies_.size();}
+
+   inline iterator       begin   ()       {return bodies_.begin();}
+   inline const_iterator begin   () const {return bodies_.begin();}
+   inline const_iterator cbegin  () const {return bodies_.cbegin();}
+   inline iterator       end     ()       {return bodies_.end();}
+   inline const_iterator end     () const {return bodies_.end();}
+   inline const_iterator cend    () const {return bodies_.cend();}
+
+   template< typename C >
+   inline cast_iterator<C> begin()              {return bodies_.begin<C>();}
+   template< typename C >
+   inline const_cast_iterator<C> begin() const  {return bodies_.begin<C>();}
+   template< typename C >
+   inline const_cast_iterator<C> cbegin() const {return bodies_.cbegin<C>();}
+   template< typename C >
+   inline cast_iterator<C> end()                {return bodies_.end<C>();}
+   template< typename C >
+   inline const_cast_iterator<C> end() const    {return bodies_.end<C>();}
+   template< typename C >
+   inline const_cast_iterator<C> cend() const   {return bodies_.cend<C>();}
    //@}
    //**********************************************************************************************
 
+   virtual inline real_t getVolume()         const;
+
    //**Set functions*******************************************************************************
    /*!\name Set functions */
    //@{
@@ -146,7 +157,7 @@ public:
    //**Rigid body manager functions****************************************************************
    /*!\name Rigid body manager functions */
    //@{
-   void add   ( BodyID body );
+   inline RigidBody& add( std::unique_ptr<RigidBody>&& body );
    //@}
    //**********************************************************************************************
 
@@ -183,18 +194,13 @@ protected:
    /*!\name Utility functions */
    //@{
    inline virtual void calcBoundingBox() WALBERLA_OVERRIDE;  // Calculation of the axis-aligned bounding box
-                  void calcCenterOfMass(); ///< updates the center of mass (gpos)
+   void calcCenterOfMass(); ///< updates the center of mass (gpos)
    inline         void calcInertia();      // Calculation of the moment of inertia
    //@}
    //**********************************************************************************************
 
-   //**Member variables****************************************************************************
-   /*!\name Member variables */
-   //@{
-   PtrVector<RigidBody, PtrDelete> bodies_;  //!< Rigid bodies contained in the union.
-   //@}
-   //**********************************************************************************************
 private:
+   BodyStorage bodies_;  //!< Rigid bodies contained in the union.
    std::vector<id_t> containedTypeIDs_;
 
    static id_t staticTypeID_;  //< type id of sphere, will be set by SetBodyTypeIDs
@@ -228,7 +234,7 @@ private:
  */
 template <typename BodyTypeTuple>
 Union<BodyTypeTuple>::Union( id_t sid, id_t uid, const Vec3& gpos, const Vec3& rpos, const Quat& q,
-                const bool global, const bool communicating, const bool /*infiniteMass*/ )
+                             const bool global, const bool communicating, const bool /*infiniteMass*/ )
    : RigidBody( getStaticTypeID(), sid, uid )  // Initialization of the parent class
 {
    // Initializing the instantiated union
@@ -241,7 +247,7 @@ Union<BodyTypeTuple>::Union( id_t sid, id_t uid, const Vec3& gpos, const Vec3& r
    calcInertia();
 
    setGlobal( global );
-//   setMass( infiniteMass );
+   //   setMass( infiniteMass );
    setCommunicating( communicating );
    setFinite( true );
 }
@@ -322,7 +328,7 @@ void Union<BodyTypeTuple>::calcBoundingBox()
    // Setting the bounding box of an empty union
    if( bodies_.isEmpty() ) {
       aabb_ = math::AABB(
-            gpos_[0] - real_t(0.01),
+                 gpos_[0] - real_t(0.01),
             gpos_[1] - real_t(0.01),
             gpos_[2] - real_t(0.01),
             gpos_[0] + real_t(0.01),
@@ -334,7 +340,7 @@ void Union<BodyTypeTuple>::calcBoundingBox()
    // Using the bounding box of the first contained bodies as initial bounding box
    // and merging it with the bounding boxes of all other bodies
    else {
-      aabb_ = bodies_[0]->getAABB();
+      aabb_ = bodies_.begin()->getAABB();
       for( auto b=bodies_.begin()+1; b!=bodies_.end(); ++b )
          aabb_.merge( b->getAABB() );
    }
@@ -366,7 +372,7 @@ void Union<BodyTypeTuple>::calcCenterOfMass()
    // Calculating the center of mass of a single body
    if( bodies_.size() == 1 )
    {
-      const BodyID body( bodies_[0] );
+      const BodyID body( bodies_.begin().getBodyID() );
       gpos_ = body->getPosition();
       mass_ = body->getMass();
       if( !isFixed() && mass_ > real_t(0) )
@@ -512,8 +518,8 @@ void Union<BodyTypeTuple>::setPositionImpl( real_t px, real_t py, real_t pz )
    gpos_ = gpos;
 
    // Updating the contained bodies
-   for( auto b=bodies_.begin(); b!=bodies_.end(); ++b )
-      (*b)->update( dp );
+   for( auto& b : bodies_ )
+      b.update( dp );
 
    Union<BodyTypeTuple>::calcBoundingBox();    // Setting the axis-aligned bounding box
    wake();               // Waking the union from sleep mode
@@ -548,8 +554,8 @@ void Union<BodyTypeTuple>::setOrientationImpl( real_t r, real_t i, real_t j, rea
    R_ = q_.toRotationMatrix();
 
    // Updating the contained bodies
-   for( auto b=bodies_.begin(); b!=bodies_.end(); ++b )
-      (*b)->update( dq );
+   for( auto& b : bodies_ )
+      b.update( dq );
 
    Union<BodyTypeTuple>::calcBoundingBox();  // Setting the axis-aligned bounding box
    wake();             // Waking the union from sleep mode
@@ -586,8 +592,8 @@ void Union<BodyTypeTuple>::update( const Vec3& dp )
    gpos_ += dp;
 
    // Updating the contained bodies
-   for( auto b=bodies_.begin(); b!=bodies_.end(); ++b )
-      (*b)->update( dp );
+   for( auto& b : bodies_ )
+      b.update( dp );
 
    // Setting the axis-aligned bounding box
    Union<BodyTypeTuple>::calcBoundingBox();
@@ -623,8 +629,8 @@ void Union<BodyTypeTuple>::update( const Quat& dq )
    R_ = q_.toRotationMatrix();
 
    // Updating the contained bodies
-   for( auto b=bodies_.begin(); b!=bodies_.end(); ++b )
-      (*b)->update( dq );
+   for( auto& b : bodies_ )
+      b.update( dq );
 
    // Setting the axis-aligned bounding box
    Union<BodyTypeTuple>::calcBoundingBox();
@@ -694,38 +700,38 @@ void Union<BodyTypeTuple>::update( const Quat& dq )
  *    union (to make the entire union (in-)visible.
  */
 template <typename BodyTypeTuple>
-void Union<BodyTypeTuple>::add( BodyID body )
+RigidBody& Union<BodyTypeTuple>::add( std::unique_ptr<RigidBody>&& body )
 {
    // Checking for "self-assignment"
-   if( body == BodyID( this ) ) return;
+   if( body.get() == BodyID( this ) ) return *this;
 
    // Checking the global flags of the body and the union in MPI parallel simulations
    if( body->isGlobal() ^ global_ )
       throw std::logic_error( "Global flags of body and union do not match" );
 
+   // Registering the rigid body
+   auto& bd = bodies_.add( std::move(body) );
+
    Vec3 oldCenterOfMass = getPosition();
    Vec3 oldImpulse      = getLinearVel() * getMass();
 
-   Vec3 bodyCenterOfMass = body->getPosition();
-   Vec3 bodyImpulse      = body->getLinearVel() * body->getMass();
+   Vec3 bodyCenterOfMass = bd.getPosition();
+   Vec3 bodyImpulse      = bd.getLinearVel() * bd.getMass();
 
-
-   // Registering the rigid body
-   body->setSB(this);
-   bodies_.pushBack( body );
+   bd.setSB(this); //having a superbody will forward all getVel calls to superbody!!!
 
    // Updating the axis-aligned bounding box
    if( bodies_.size() == 1 )
-      aabb_ = body->getAABB();
+      aabb_ = bd.getAABB();
    else
-      aabb_.merge( body->getAABB() );
+      aabb_.merge( bd.getAABB() );
 
    // Setting the union's total mass and center of mass
    calcCenterOfMass();
 
    // Setting the contained primitives' relative position in reference to the center of mass
-   for( auto b=bodies_.begin(); b!=bodies_.end(); ++b )
-      ( *b )->calcRelPosition();
+   for( auto& b : bodies_ )
+      b.calcRelPosition();
 
    // Setting the moment of inertia
    calcInertia();
@@ -735,11 +741,13 @@ void Union<BodyTypeTuple>::add( BodyID body )
    setAngularVel( Vec3(0,0,0) );
    addImpulseAtPos( oldImpulse, oldCenterOfMass);
    addImpulseAtPos( bodyImpulse, bodyCenterOfMass);
-   body->setLinearVel( Vec3(0,0,0) );
-   body->setAngularVel( Vec3(0,0,0) );
+   bd.setLinearVel( Vec3(0,0,0) );
+   bd.setAngularVel( Vec3(0,0,0) );
 
    // Signaling the internal modification to the superordinate body
    signalModification();
+
+   return bodies_.back();
 }
 //*************************************************************************************************
 
@@ -771,8 +779,8 @@ void Union<BodyTypeTuple>::translateImpl( real_t dx, real_t dy, real_t dz )
    gpos_ += dp;
 
    // Updating the contained bodies
-   for( auto b=bodies_.begin(); b!=bodies_.end(); ++b )
-      (*b)->update( dp );
+   for( auto& b : bodies_ )
+      b.update( dp );
 
    Union<BodyTypeTuple>::calcBoundingBox();    // Setting the axis-aligned bounding box
    wake();               // Waking the union from sleep mode
@@ -811,8 +819,8 @@ void Union<BodyTypeTuple>::rotateImpl( const Quat& dq )
    R_ = q_.toRotationMatrix();  // Updating the rotation of the union
 
    // Updating the contained bodies
-   for( auto b=bodies_.begin(); b!=bodies_.end(); ++b )
-      (*b)->update( dq );
+   for( auto& b : bodies_ )
+      b.update( dq );
 
    Union<BodyTypeTuple>::calcBoundingBox();  // Setting the axis-aligned bounding box
    wake();             // Waking the union from sleep mode
@@ -843,8 +851,8 @@ void Union<BodyTypeTuple>::rotateAroundOriginImpl( const Quat& dq )
    gpos_ = dq.rotate( gpos_ );     // Updating the global position of the union
 
    // Updating the contained bodies
-   for( auto b=bodies_.begin(); b!=bodies_.end(); ++b )
-      (*b)->update( dq );
+   for( auto& b : bodies_ )
+      b.update( dq );
 
    Union<BodyTypeTuple>::calcBoundingBox();    // Setting the axis-aligned bounding box
    wake();               // Waking the union from sleep mode
@@ -876,8 +884,8 @@ void Union<BodyTypeTuple>::rotateAroundPointImpl( const Vec3& point, const Quat
    gpos_ = point + dq.rotate( dp );  // Updating the global position of the union
 
    // Updating the contained bodies
-   for( auto b=bodies_.begin(); b!=bodies_.end(); ++b )
-      (*b)->update( dq );
+   for( auto& b : bodies_ )
+      b.update( dq );
 
    Union<BodyTypeTuple>::calcBoundingBox();    // Setting the axis-aligned bounding box
    wake();               // Waking the union from sleep mode
@@ -906,8 +914,8 @@ template <typename BodyTypeTuple>
 bool Union<BodyTypeTuple>::containsRelPointImpl( real_t px, real_t py, real_t pz ) const
 {
    const Vec3 gpos( pointFromBFtoWF( px, py, pz ) );
-   for( auto b=bodies_.begin(); b!=bodies_.end(); ++b )
-      if( b->containsPoint( gpos ) ) return true;
+   for( auto& b : bodies_ )
+      if( b.containsPoint( gpos ) ) return true;
    return false;
 }
 //*************************************************************************************************
@@ -929,9 +937,10 @@ bool Union<BodyTypeTuple>::isSurfaceRelPointImpl( real_t px, real_t py, real_t p
    bool surface( false );
    const Vec3 gpos( pointFromBFtoWF( px, py, pz ) );
 
-   for( auto b=bodies_.begin(); b!=bodies_.end(); ++b ) {
-      if( b->containsPoint( gpos ) ) return false;
-      else if( b->isSurfacePoint( gpos ) ) surface = true;
+   for( auto& b : bodies_ )
+   {
+      if( b.containsPoint( gpos ) ) return false;
+      else if( b.isSurfacePoint( gpos ) ) surface = true;
    }
 
    return surface;
@@ -974,8 +983,8 @@ void Union<BodyTypeTuple>::handleModification()
    calcCenterOfMass();
 
    // Setting the contained primitives' relative position in reference to the center of mass
-   for( auto b=bodies_.begin(); b!=bodies_.end(); ++b )
-      ( *b )->calcRelPosition();
+   for( auto& b : bodies_ )
+      b.calcRelPosition();
 
    // Setting the moment of inertia
    calcInertia();
@@ -1004,8 +1013,8 @@ void Union<BodyTypeTuple>::handleTranslation()
    calcCenterOfMass();
 
    // Setting the contained bodies' relative position in reference to the center of mass
-   for( auto b=bodies_.begin(); b!=bodies_.end(); ++b )
-      ( *b )->calcRelPosition();
+   for( auto& b : bodies_ )
+      b.calcRelPosition();
 
    // Setting the moment of inertia
    calcInertia();
@@ -1155,65 +1164,6 @@ std::ostream& operator<<( std::ostream& os, Union<BodyTypeTuple> const * u )
 //
 //=================================================================================================
 
-//*************************************************************************************************
-/*!\brief Returns whether the union contains any bodies or not.
- *
- * \return \a true if the union is empty, \a false if it not.
- */
-template <typename BodyTypeTuple>
-inline bool Union<BodyTypeTuple>::isEmpty() const
-{
-   return bodies_.size() == 0;
-}
-//*************************************************************************************************
-
-
-//*************************************************************************************************
-/*!\brief Returns the number of rigid bodies contained in the union.
- *
- * \return The number of rigid bodies.
- */
-template <typename BodyTypeTuple>
-inline size_t Union<BodyTypeTuple>::size() const
-{
-   return bodies_.size();
-}
-//*************************************************************************************************
-
-
-//*************************************************************************************************
-/*!\brief Returns a handle to the indexed rigid body.
- *
- * \param index Access index. The index has to be in the range \f$[0..size-1]\f$.
- * \return Handle to the accessed body.
- *
- * \b Note: No runtime check is performed to insure the validity of the access index.
- */
-template <typename BodyTypeTuple>
-inline BodyID Union<BodyTypeTuple>::getBody( size_t index )
-{
-   WALBERLA_ASSERT_LESS( index, bodies_.size(), "Invalid access index" );
-   return bodies_[index];
-}
-//*************************************************************************************************
-
-
-//*************************************************************************************************
-/*!\brief Returns a constant handle to the indexed rigid body.
- *
- * \param index Access index. The index has to be in the range \f$[0..size-1]\f$.
- * \return Constant handle to the accessed body.
- *
- * \b Note: No runtime check is performed to insure the validity of the access index.
- */
-template <typename BodyTypeTuple>
-inline ConstBodyID Union<BodyTypeTuple>::getBody( size_t index ) const
-{
-   WALBERLA_ASSERT_LESS( index, bodies_.size(), "Invalid access index" );
-   return bodies_[index];
-}
-//*************************************************************************************************
-
 
 //*************************************************************************************************
 /*!\brief Returns unique type id of this type
diff --git a/src/pe/rigidbody/UnionFactory.h b/src/pe/rigidbody/UnionFactory.h
index 4cb41e9a3..45c01d4a4 100644
--- a/src/pe/rigidbody/UnionFactory.h
+++ b/src/pe/rigidbody/UnionFactory.h
@@ -145,7 +145,6 @@ BoxID createBox( Union<BodyTypeTuple>* un,
    if( lengths[0] <= real_t(0) || lengths[1] <= real_t(0) || lengths[2] <= real_t(0) )
       throw std::invalid_argument( "Invalid side length" );
 
-   BoxID box = NULL;
    id_t  sid = 0;
 
    if (global)
@@ -158,9 +157,8 @@ BoxID createBox( Union<BodyTypeTuple>* un,
       sid = UniqueID<RigidBody>::create();
    }
 
-   box = new Box(sid, uid, gpos, Vec3(0,0,0), Quat(), lengths, material, global, communicating, infiniteMass);
+   std::unique_ptr<Box> box = std::make_unique<Box>(sid, uid, gpos, Vec3(0,0,0), Quat(), lengths, material, global, communicating, infiniteMass);
    box->MPITrait.setOwner( un->MPITrait.getOwner() );
-   un->add(box);
 
    if (box != NULL)
    {
@@ -175,7 +173,7 @@ BoxID createBox( Union<BodyTypeTuple>* un,
                );
    }
 
-   return box;
+   return static_cast<BoxID> (&(un->add(std::move(box))));
 }
 
 //*************************************************************************************************
@@ -217,7 +215,6 @@ CapsuleID createCapsule( Union<BodyTypeTuple>* un,
    if( length <= real_c(0) )
       throw std::invalid_argument( "Invalid capsule length" );
 
-   CapsuleID capsule = NULL;
    id_t      sid     = 0;
 
    if (global)
@@ -230,16 +227,15 @@ CapsuleID createCapsule( Union<BodyTypeTuple>* un,
       sid = UniqueID<RigidBody>::create();
    }
 
-   capsule = new Capsule(sid, uid, gpos, Vec3(0,0,0), Quat(), radius, length, material, global, communicating, infiniteMass);
+   std::unique_ptr<Capsule> capsule = std::make_unique<Capsule>(sid, uid, gpos, Vec3(0,0,0), Quat(), radius, length, material, global, communicating, infiniteMass);
    capsule->MPITrait.setOwner( un->MPITrait.getOwner() );
-   un->add(capsule);
 
    if (capsule != NULL)
    {
-      WALBERLA_LOG_DETAIL("Created capsule " << capsule->getSystemID() << "\n" << capsule);
+      WALBERLA_LOG_DETAIL("Created capsule " << capsule->getSystemID() << "\n" << *capsule);
    }
 
-   return capsule;
+   return static_cast<CapsuleID>(&(un->add(std::move(capsule))));
 }
 
 //*************************************************************************************************
@@ -277,7 +273,6 @@ SphereID createSphere( Union<BodyTypeTuple>* un,
       throw std::invalid_argument( "Invalid sphere radius" );
 
    id_t sid(0);
-   SphereID sphere = NULL;
 
    if (global)
    {
@@ -290,9 +285,8 @@ SphereID createSphere( Union<BodyTypeTuple>* un,
       sid = UniqueID<RigidBody>::create();
    }
 
-   sphere = new Sphere(sid, uid, gpos, Vec3(0,0,0), Quat(), radius, material, global, communicating, infiniteMass);
+   std::unique_ptr<Sphere> sphere = std::make_unique<Sphere>(sid, uid, gpos, Vec3(0,0,0), Quat(), radius, material, global, communicating, infiniteMass);
    sphere->MPITrait.setOwner( un->MPITrait.getOwner() );
-   un->add( sphere );
 
    if (sphere != NULL)
    {
@@ -307,7 +301,7 @@ SphereID createSphere( Union<BodyTypeTuple>* un,
                );
    }
 
-   return sphere;
+   return static_cast<SphereID>(&(un->add( std::move(sphere) )));
 }
 
 }  // namespace pe
diff --git a/tests/pe/Collision.cpp b/tests/pe/Collision.cpp
index 98ce1fc4a..54d4f1125 100644
--- a/tests/pe/Collision.cpp
+++ b/tests/pe/Collision.cpp
@@ -31,6 +31,7 @@
 #include "pe/rigidbody/Sphere.h"
 #include "pe/rigidbody/Plane.h"
 #include "pe/rigidbody/Union.h"
+#include "pe/rigidbody/UnionFactory.h"
 
 #include "pe/rigidbody/SetBodyTypeIDs.h"
 #include "pe/Types.h"
@@ -273,13 +274,10 @@ void CapsuleTest2()
 void UnionTest()
 {
    typedef Union< boost::tuple<Sphere> > UnionT;
-   MaterialID iron = Material::find("iron");
    UnionT  un1(120, 0, Vec3(0,0,0), Vec3(0,0,0), Quat(), false, true, false);
    UnionT  un2(121, 0, Vec3(real_t(1.5),0,0), Vec3(0,0,0), Quat(), false, true, false);
-   SphereID sp1 = new Sphere(123, 1, Vec3(0,0,0), Vec3(0,0,0), Quat(), 1, iron, false, true, false);
-   un1.add(sp1);
-   SphereID sp2 = new Sphere(124, 2, Vec3(real_t(1.5),0,0), Vec3(0,0,0), Quat(), 1, iron, false, true, false);
-   un2.add(sp2);
+   auto sp1 = createSphere(&un1, 123, Vec3(0,0,0), 1);
+   auto sp2 = createSphere(&un2, 124, Vec3(real_t(1.5),0,0), 1);
 
    std::vector<Contact> contacts;
 
diff --git a/tests/pe/CollisionTobiasGJK.cpp b/tests/pe/CollisionTobiasGJK.cpp
index 2c3fa052a..f1710053b 100644
--- a/tests/pe/CollisionTobiasGJK.cpp
+++ b/tests/pe/CollisionTobiasGJK.cpp
@@ -32,6 +32,7 @@
 #include "pe/rigidbody/Sphere.h"
 #include "pe/rigidbody/Plane.h"
 #include "pe/rigidbody/Union.h"
+#include "pe/rigidbody/UnionFactory.h"
 #include "pe/rigidbody/Ellipsoid.h"
 
 #include "pe/rigidbody/SetBodyTypeIDs.h"
@@ -355,19 +356,16 @@ void UnionTest(){
    Box box(179, 179, Vec3(0,0,0), Vec3(0,0,0), Quat(), Vec3(real_t(10),real_t(2), real_t(10)), iron, false, true, false);
 
 
-   Union<boost::tuple<Sphere>> *unsub = new Union<boost::tuple<Sphere>>(192, 192, Vec3(0,real_t(3.8),0), Vec3(0,0,0), Quat(), false, true, false);
+   using UnionT = Union<boost::tuple<Sphere>>;
+   auto unsub = std::make_unique<UnionT>(192, 192, Vec3(0,real_t(3.8),0), Vec3(0,0,0), Quat(), false, true, false);
 
-   Sphere sp1( 180, 180, Vec3(-3,real_t(3.8),0), Vec3(0,0,0), Quat(), real_t(3.0)  , iron, false, true, false );
-   Sphere sp2( 181, 181, Vec3(3,real_t(3.8),0), Vec3(0,0,0), Quat(), real_t(3.0), iron, false, true, false );
-
-   Sphere sp3( 182, 182, Vec3(0,real_t(6),0), Vec3(0,0,0), Quat(), real_t(3.0), iron, false, true, false );
-   unsub->add(&sp1);
-   unsub->add(&sp2);
+   auto sp1 = createSphere(unsub.get(), 180, Vec3(-3,real_t(3.8),0), real_t(3.0));
+   auto sp2 = createSphere(unsub.get(), 181, Vec3(3,real_t(3.8),0), real_t(3.0));
 
    //Create another union, and add sub union
    Union<boost::tuple<Sphere, Union<boost::tuple<Sphere>>>> *un = new Union<boost::tuple<Sphere, Union<boost::tuple<Sphere>>>>(193, 193, Vec3(0, 0, 0), Vec3(0,0,0), Quat(), false, true, false);
-   un->add(&sp3);
-   un->add(unsub);
+   createSphere(un, 182, Vec3(0,real_t(6),0), real_t(3.0));
+   un->add(std::move(unsub));
 
 
    PossibleContacts pcs;
@@ -378,9 +376,9 @@ void UnionTest(){
    Contact &c = container.back();
    WALBERLA_LOG_DEVEL( c.getDistance() << " " << c.getNormal() << " " << c.getPosition() );
    if(c.getBody1()->getID() == 181) {
-      checkContact( c, Contact(&sp2, &box,  Vec3(real_t(3), real_t(0.9), 0), Vec3(0, 1, 0), real_t(-0.2)), Vec3(0,0,0));
+      checkContact( c, Contact(sp2, &box,  Vec3(real_t(3), real_t(0.9), 0), Vec3(0, 1, 0), real_t(-0.2)), Vec3(0,0,0));
    } else if (c.getBody1()->getID() == 179) {
-      checkContact( c, Contact(&box, &sp2,  Vec3(real_t(3), real_t(0.9), 0), Vec3(0, -1, 0), real_t(-0.2)), Vec3(0,0,0));
+      checkContact( c, Contact(&box, sp2,  Vec3(real_t(3), real_t(0.9), 0), Vec3(0, -1, 0), real_t(-0.2)), Vec3(0,0,0));
    } else {
       WALBERLA_ABORT("Unknown ID!");
    }
@@ -390,9 +388,9 @@ void UnionTest(){
    c = container.back();
    WALBERLA_LOG_DEVEL( c.getDistance() << " " << c.getNormal() << " " << c.getPosition() );
    if(c.getBody1()->getID() == 180) {
-      checkContact( c, Contact(&sp1, &box,  Vec3(real_t(-3), real_t(0.9), 0), Vec3(0, 1, 0), real_t(-0.2)), Vec3(0,0,0));
+      checkContact( c, Contact(sp1, &box,  Vec3(real_t(-3), real_t(0.9), 0), Vec3(0, 1, 0), real_t(-0.2)), Vec3(0,0,0));
    } else if (c.getBody1()->getID() == 179) {
-      checkContact( c, Contact(&box, &sp1,  Vec3(real_t(-3), real_t(0.9), 0), Vec3(0, -1, 0), real_t(-0.2)), Vec3(0,0,0));
+      checkContact( c, Contact(&box, sp1,  Vec3(real_t(-3), real_t(0.9), 0), Vec3(0, -1, 0), real_t(-0.2)), Vec3(0,0,0));
    } else {
       WALBERLA_ABORT("Unknown ID!");
    }
@@ -406,9 +404,9 @@ void UnionTest(){
    c = container.back();
    WALBERLA_LOG_DEVEL( c.getDistance() << " " << c.getNormal() << " " << c.getPosition() );
    if(c.getBody1()->getID() == 181) {
-      checkContact( c, Contact(&sp2, &box,  Vec3(real_t(3), real_t(0.9), 0), Vec3(0, 1, 0), real_t(-0.2)), Vec3(0,0,0));
+      checkContact( c, Contact(sp2, &box,  Vec3(real_t(3), real_t(0.9), 0), Vec3(0, 1, 0), real_t(-0.2)), Vec3(0,0,0));
    } else if (c.getBody1()->getID() == 179) {
-      checkContact( c, Contact(&box, &sp2,  Vec3(real_t(3), real_t(0.9), 0), Vec3(0, -1, 0), real_t(-0.2)), Vec3(0,0,0));
+      checkContact( c, Contact(&box, sp2,  Vec3(real_t(3), real_t(0.9), 0), Vec3(0, -1, 0), real_t(-0.2)), Vec3(0,0,0));
    } else {
       WALBERLA_ABORT("Unknown ID!");
    }
@@ -417,9 +415,9 @@ void UnionTest(){
    c = container.back();
    WALBERLA_LOG_DEVEL( c.getDistance() << " " << c.getNormal() << " " << c.getPosition() );
    if(c.getBody1()->getID() == 180) {
-      checkContact( c, Contact(&sp1, &box,  Vec3(real_t(-3), real_t(0.9), 0), Vec3(0, 1, 0), real_t(-0.2)), Vec3(0,0,0));
+      checkContact( c, Contact(sp1, &box,  Vec3(real_t(-3), real_t(0.9), 0), Vec3(0, 1, 0), real_t(-0.2)), Vec3(0,0,0));
    } else if (c.getBody1()->getID() == 179) {
-      checkContact( c, Contact(&box, &sp1,  Vec3(real_t(-3), real_t(0.9), 0), Vec3(0, -1, 0), real_t(-0.2)), Vec3(0,0,0));
+      checkContact( c, Contact(&box, sp1,  Vec3(real_t(-3), real_t(0.9), 0), Vec3(0, -1, 0), real_t(-0.2)), Vec3(0,0,0));
    } else {
       WALBERLA_ABORT("Unknown ID!");
    }
diff --git a/tests/pe/Marshalling.cpp b/tests/pe/Marshalling.cpp
index 30da9d3be..b815a5d96 100644
--- a/tests/pe/Marshalling.cpp
+++ b/tests/pe/Marshalling.cpp
@@ -170,8 +170,8 @@ void testUnion()
    WALBERLA_CHECK_EQUAL(u1.size(), u2->size());
 
    //getting spheres of second union
-   SphereID s12 = static_cast<SphereID> (*(u2->begin()));
-   SphereID s22 = static_cast<SphereID> (*(++(u2->begin())));
+   SphereID s12 = static_cast<SphereID> (u2->begin().getBodyID());
+   SphereID s22 = static_cast<SphereID> ((++(u2->begin())).getBodyID());
    WALBERLA_CHECK_UNEQUAL( s12, s22 );
 
    WALBERLA_CHECK_FLOAT_EQUAL( s11->getPosition(),    s12->getPosition());
diff --git a/tests/pe/ShadowCopy.cpp b/tests/pe/ShadowCopy.cpp
index da6c31cef..8c5dba7c4 100644
--- a/tests/pe/ShadowCopy.cpp
+++ b/tests/pe/ShadowCopy.cpp
@@ -19,6 +19,7 @@
 //======================================================================================================================
 
 #include "pe/basic.h"
+#include "pe/rigidbody/UnionFactory.h"
 #include "pe/utility/GetBody.h"
 #include "pe/utility/DestroyBody.h"
 
@@ -95,12 +96,9 @@ int main( int argc, char** argv )
    destroyBodyBySID( *globalBodyStorage, forest->getBlockStorage(), storageID, sid );
 
    WALBERLA_LOG_PROGRESS_ON_ROOT( " *** UNION *** ");
-   MaterialID iron = Material::find("iron");
    UnionT* un   = createUnion< boost::tuple<Sphere> >( *globalBodyStorage, forest->getBlockStorage(), storageID, 0, Vec3(2,2,2) );
-   SphereID sp1 = new Sphere( 10, 0, Vec3(real_t(4.9),2,2), Vec3(0,0,0), Quat(), real_t(1)  , iron, false, true, false );
-   SphereID sp2 = new Sphere( 11, 0, Vec3(3,2,2)          , Vec3(0,0,0), Quat(), real_t(1.5), iron, false, true, false );
-   un->add(sp1);
-   un->add(sp2);
+   auto sp1 = createSphere(un, 10, Vec3(real_t(4.9),2,2), real_t(1));
+   auto sp2 = createSphere(un, 11, Vec3(3,2,2), real_t(1.5));
    un->setPosition( Vec3( real_t(4.9), 2, 2) );
    auto relPosSp1 = sp1->getRelPosition();
    auto relPosSp2 = sp2->getRelPosition();
@@ -113,8 +111,8 @@ int main( int argc, char** argv )
    syncCall();
 
    un  = static_cast<UnionT*> (getBody( *globalBodyStorage, forest->getBlockStorage(), storageID, sid, StorageSelect::LOCAL ));
-   sp1 = static_cast<SphereID> (*(un->begin()));
-   sp2 = static_cast<SphereID> (*(++(un->begin())));
+   sp1 = static_cast<SphereID> (un->begin().getBodyID());
+   sp2 = static_cast<SphereID> ((++(un->begin())).getBodyID());
    WALBERLA_CHECK_NOT_NULLPTR(sp1);
    WALBERLA_CHECK_NOT_NULLPTR(sp2);
    WALBERLA_CHECK_EQUAL( sp1->getTypeID(), Sphere::getStaticTypeID() );
@@ -152,8 +150,8 @@ int main( int argc, char** argv )
 
    posUnion = Vec3(real_t(0.9),2,2);
    un  = static_cast<UnionT*> (getBody( *globalBodyStorage, forest->getBlockStorage(), storageID, sid, StorageSelect::LOCAL ));
-   sp1 = static_cast<SphereID> (*(un->begin()));
-   sp2 = static_cast<SphereID> (*(++(un->begin())));
+   sp1 = static_cast<SphereID> (un->begin().getBodyID());
+   sp2 = static_cast<SphereID> ((++(un->begin())).getBodyID());
    WALBERLA_CHECK_NOT_NULLPTR(sp1);
    WALBERLA_CHECK_NOT_NULLPTR(sp2);
    WALBERLA_CHECK_EQUAL( sp1->getTypeID(), Sphere::getStaticTypeID() );
diff --git a/tests/pe/Union.cpp b/tests/pe/Union.cpp
index 6cd0aec53..1a814c28a 100644
--- a/tests/pe/Union.cpp
+++ b/tests/pe/Union.cpp
@@ -25,6 +25,7 @@
 
 #include "pe/basic.h"
 #include "pe/rigidbody/Union.h"
+#include "pe/rigidbody/UnionFactory.h"
 #include "pe/ccd/SimpleCCDDataHandling.h"
 #include "pe/synchronization/SyncNextNeighbors.h"
 #include "pe/vtk/BodyVtkOutput.h"
@@ -75,12 +76,9 @@ void SnowManFallingOnPlane()
 
    createPlane( *globalBodyStorage, 0, Vec3(0,0,1), Vec3(0,0,0) );
 
-   MaterialID iron = Material::find("iron");
-   UnionType* un   = createUnion< boost::tuple<Sphere> >( *globalBodyStorage, forest->getBlockStorage(), storageID, 0, Vec3(5,5,5) );
-   SphereID sp1 = new Sphere( 10, 0, Vec3(5,5,1), Vec3(0,0,0), Quat(), real_t(1)  , iron, false, true, false );
-   SphereID sp2 = new Sphere( 11, 0, Vec3(real_t(6.7),5,real_t(1.2)), Vec3(0,0,0), Quat(), real_t(1.1), iron, false, true, false );
-   un->add(sp1);
-   un->add(sp2);
+   UnionType* un  = createUnion< boost::tuple<Sphere> >( *globalBodyStorage, forest->getBlockStorage(), storageID, 0, Vec3(5,5,5) );
+   auto sp1       = createSphere(un, 10, Vec3(5,5,1), real_t(1));
+   auto sp2       = createSphere(un, 11, Vec3(real_t(6.7),5,real_t(1.2)), real_t(1.1));
 
    auto distance = (sp1->getPosition() - sp2->getPosition()).length();
 
@@ -106,14 +104,15 @@ void ImpulsCarryover()
 {
    MaterialID iron = Material::find("iron");
 
-   UnionType* un   = new Union< boost::tuple<Sphere> >(12, 0, Vec3(0,0,0), Vec3(0,0,0), Quat(), false, true, false);
-   SphereID sp1    = new Sphere( 10, 0, Vec3( 1,0,0), Vec3(0,0,0), Quat(), real_t(1), iron, false, true, false );
-   SphereID sp2    = new Sphere( 11, 0, Vec3(-1,0,0), Vec3(0,0,0), Quat(), real_t(1), iron, false, true, false );
+   auto un  = std::make_unique<UnionType>(12, 0, Vec3(0,0,0), Vec3(0,0,0), Quat(), false, true, false);
+   auto sp1 = std::make_unique<Sphere>( 10, 0, Vec3( 1,0,0), Vec3(0,0,0), Quat(), real_t(1), iron, false, true, false );
+   auto sp2 = std::make_unique<Sphere>( 11, 0, Vec3(-1,0,0), Vec3(0,0,0), Quat(), real_t(1), iron, false, true, false );
+
    sp1->setLinearVel(Vec3(0,real_c(+1),0));
    sp2->setLinearVel(Vec3(0,real_c(-1),0));
 
-   un->add(sp1);
-   un->add(sp2);
+   un->add( std::move(sp1) );
+   un->add( std::move(sp2) );
 
    WALBERLA_CHECK_FLOAT_EQUAL( un->getPosition(),  Vec3(0,0,0) );
    WALBERLA_CHECK_FLOAT_EQUAL( un->getLinearVel(), Vec3(0,0,0) );
-- 
GitLab