diff --git a/src/pe/collision/EPA.h b/src/pe/collision/EPA.h
index 6abe12a8ea75ae958fe745406fcc33d033e19629..c73183941b878f3f229c596dd3b2964e351f2915 100644
--- a/src/pe/collision/EPA.h
+++ b/src/pe/collision/EPA.h
@@ -90,6 +90,20 @@ public:
    //@}
    //**********************************************************************************************
 
+   //**Getter/Setter functions*****************************************************************************
+   /*!\name Getter and Setter functions */
+   //@{
+   inline void setMaxSupportPoints( size_t maxSupportPoints) {maxSupportPoints_ = maxSupportPoints;}
+
+   inline size_t getMaxSupportPoints() {return maxSupportPoints_;}
+
+   inline void setMaxTriangles( size_t maxTriangles) {maxTriangles_ = maxTriangles;}
+
+   inline size_t getMaxTriangles() {return maxTriangles_;}
+
+   //@}
+   //**********************************************************************************************
+
 private:
    //**Utility functions***************************************************************************
    /*!\name Utility functions */
@@ -128,8 +142,8 @@ private:
 
 private:
    //EPA constants
-   static const size_t maxSupportPoints_ = 100;
-   static const size_t maxTriangles_     = 200;
+   size_t maxSupportPoints_ = 100;
+   size_t maxTriangles_     = 200;
 };
 //*************************************************************************************************
 
diff --git a/src/pe/fcd/IterativeFCD.h b/src/pe/fcd/IterativeFCD.h
index 995faf7ea3861adccd29098e78d2018e6d1e49ce..5634f8bb6d291ef288e6c1fd15a43e4fa32051a0 100644
--- a/src/pe/fcd/IterativeFCD.h
+++ b/src/pe/fcd/IterativeFCD.h
@@ -109,7 +109,7 @@ private:
    * \param bd1 The first body.
    * \param bd2 The second body.
    *
-   * This function is called by generateContacts() on all normal
+   * This function is called by generateContactsPair() on all normal (single)
    * body pairs and on each single bodys of the union.
    * It checks if one of the bodies is a plane and calls the corresponding
    * detection function. If collision is detected a Contact is generated and
@@ -152,6 +152,38 @@ private:
       }
    }
 
+   /* Resolve pairs with unions and call generateContacts single for each pair of non-united
+    * particles.
+    */
+   void generateContactsPair(BodyID bd1, BodyID bd2){
+      bool hasSubbodies1 = bd1->hasSubBodies();
+      if( hasSubbodies1 ){
+         UnionGenericType* u1 = nullptr;
+         u1= static_cast<UnionGenericType*>(bd1);
+         //Loop over possible subbodies of first rigid body
+         auto it1 = u1->begin();
+         while(it1 != u1->end()){
+            generateContactsPair(*it1, bd2);
+            it1++;
+         }
+      }else{
+         //bd1 is a single body
+         bool hasSubbodies2 = bd2->hasSubBodies();
+         if( hasSubbodies2 ){
+            UnionGenericType* u2 = nullptr;
+            u2= static_cast<UnionGenericType*>(bd2);
+            auto it2 = u2->begin();
+            while(it2 != u2->end()){
+               generateContactsPair(bd1, *it2);
+               it2++;
+            }
+         }else{
+            //bd1 and bd2 are single bodies
+            generateContactsSingle(bd1, bd2);
+         }
+      }
+   }
+
 public:
    /*!\brief Checks each body pair in the given vector and generates
    * a vector of contacts.
@@ -164,25 +196,7 @@ public:
       contacts_.clear();
       for (auto it = possibleContacts.begin(); it != possibleContacts.end(); ++it)
       {
-         bool hasSubbodies1 = it->first->hasSubBodies();
-         UnionGenericType* u1 = nullptr;
-         if( hasSubbodies1 ){
-            u1= static_cast<UnionGenericType*>(it->first);
-         }
-         //Loop over possible subbodies of first rigid body
-         for(size_t bodynum1 = 0; bodynum1 < (hasSubbodies1 ? u1->size() : 1); bodynum1++){
-            BodyID bd1 = (hasSubbodies1 ? u1->getBody(bodynum1) : it->first);
-            bool hasSubbodies2 = it->second->hasSubBodies();
-            UnionGenericType* u2 = nullptr;
-            if( hasSubbodies2 ){
-               u2= static_cast<UnionGenericType*>(it->second);
-            }
-            //Loop over possible subbodies of second rigid body
-            for(size_t bodynum2 = 0; bodynum2 < (hasSubbodies2 ? u2->size() : 1); bodynum2++){
-               BodyID bd2 = (hasSubbodies2 ? u2->getBody(bodynum2) : it->second);
-               generateContactsSingle(bd1, bd2);
-            }
-         }
+         generateContactsPair(it->first, it->second);
       }
       return contacts_;
    }
diff --git a/tests/pe/Marshalling.cpp b/tests/pe/Marshalling.cpp
index 57411b0d279cf2e5524170450d82688314fc6697..30da9d3be10fbad20a2e545dbfa915a85a988ac0 100644
--- a/tests/pe/Marshalling.cpp
+++ b/tests/pe/Marshalling.cpp
@@ -26,6 +26,7 @@
 #include "pe/rigidbody/Squirmer.h"
 #include "pe/rigidbody/UnionFactory.h"
 #include "pe/rigidbody/Union.h"
+#include "pe/rigidbody/Ellipsoid.h"
 #include "pe/communication/rigidbody/Squirmer.h"
 #include "pe/communication/DynamicMarshalling.h"
 #include "pe/rigidbody/SetBodyTypeIDs.h"
@@ -41,7 +42,7 @@ typedef boost::tuple<Sphere>       UnionTypeTuple;
 typedef Union< UnionTypeTuple >    UnionT;
 typedef UnionT*                    UnionID;
 
-typedef boost::tuple<Box, Capsule, Sphere, Squirmer, UnionT> BodyTuple ;
+typedef boost::tuple<Box, Capsule, Sphere, Squirmer, UnionT, Ellipsoid> BodyTuple ;
 
 void testBox()
 {
@@ -128,6 +129,28 @@ void testSquirmer()
    WALBERLA_CHECK_FLOAT_EQUAL(s1.getSquirmerBeta(), s2->getSquirmerBeta());
 }
 
+void testEllipsoid()
+{
+   MaterialID iron = Material::find("iron");
+
+   Ellipsoid e1(759847, 1234795, Vec3(real_c(1), real_c(2), real_c(3)), Vec3(0,0,0), Quat(), Vec3(3,1,5), iron, false, false, false);
+   e1.setLinearVel(Vec3(real_c(5.2), real_c(6.3), real_c(7.4)));
+   e1.setAngularVel(Vec3(real_c(1.2), real_c(2.3), real_c(3.4)));
+
+   mpi::SendBuffer sb;
+   MarshalDynamically<BodyTuple>::execute(sb, e1);
+   mpi::RecvBuffer rb(sb);
+
+   EllipsoidID e2 = static_cast<EllipsoidID> (UnmarshalDynamically<BodyTuple>::execute(rb, Ellipsoid::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(e1.getPosition(), e2->getPosition());
+   WALBERLA_CHECK_FLOAT_EQUAL(e1.getLinearVel(), e2->getLinearVel());
+   WALBERLA_CHECK_FLOAT_EQUAL(e1.getAngularVel(), e2->getAngularVel());
+   WALBERLA_CHECK_FLOAT_EQUAL(e1.getSemiAxes(), e2->getSemiAxes());
+   WALBERLA_CHECK_EQUAL(e1.getID(), e2->getID());
+   WALBERLA_CHECK_EQUAL(e1.getSystemID(), e2->getSystemID());
+}
+
 void testUnion()
 {
    UnionT u1(159, 423, Vec3(real_c(1), real_c(2), real_c(3)), Vec3(0,0,0), Quat(), false, false, false);
@@ -179,6 +202,7 @@ int main( int argc, char** argv )
    testCapsule();
    testUnion();
    testSquirmer();
+   testEllipsoid();
 
    return EXIT_SUCCESS;
 }