diff --git a/apps/tutorials/pe/01_ConfinedGas.cpp b/apps/tutorials/pe/01_ConfinedGas.cpp
index 4f0190197ff14bc26b22daf7ad2afb35dc530889..5e4bc03e8a8ae31680773fdcab2d87f4f7e89cca 100644
--- a/apps/tutorials/pe/01_ConfinedGas.cpp
+++ b/apps/tutorials/pe/01_ConfinedGas.cpp
@@ -75,7 +75,7 @@ int main( int argc, char ** argv )
    //! [StorageDataHandling]
    //! [AdditionalBlockData]
    auto ccdID               = forest->addBlockData(ccd::createHashGridsDataHandling( globalBodyStorage, storageID ), "CCD");
-   auto fcdID               = forest->addBlockData(fcd::createSimpleFCDDataHandling<BodyTypeTuple>(), "FCD");
+   auto fcdID               = forest->addBlockData(fcd::createGenericFCDDataHandling<BodyTypeTuple, fcd::AnalyticCollideFunctor>(), "FCD");
    //! [AdditionalBlockData]
 
    WALBERLA_LOG_INFO("*** INTEGRATOR ***");
diff --git a/src/blockforest/BlockDataHandling.h b/src/blockforest/BlockDataHandling.h
index 099cb9834d83a7f616683c7cb229e14dc42b40c0..a9e612e84560bdf5593ef00d979c13d8af5bf5b0 100644
--- a/src/blockforest/BlockDataHandling.h
+++ b/src/blockforest/BlockDataHandling.h
@@ -69,7 +69,7 @@ public:
    T * deserialize( IBlock * const block ) { return this->initialize( block ); }
    T * deserializeCoarseToFine( Block * const block ) { return this->initialize( block ); }
    T * deserializeFineToCoarse( Block * const block ) { return this->initialize( block ); }
-   
+
    void deserialize( IBlock * const, const BlockDataID &, mpi::RecvBuffer & ) {}
    void deserializeCoarseToFine( Block * const, const BlockDataID &, mpi::RecvBuffer & ) {}
    void deserializeFineToCoarse( Block * const, const BlockDataID &, mpi::RecvBuffer &, const uint_t ) {}
@@ -77,6 +77,17 @@ public:
 
 
 
+template< typename T >
+class AlwaysCreateBlockDataHandling : public AlwaysInitializeBlockDataHandling<T>
+{
+public:
+   ~AlwaysCreateBlockDataHandling() {}
+
+   T * initialize( IBlock * const /*block*/ ) {return new T();}
+};
+
+
+
 namespace internal {
 
 
diff --git a/src/pe/collision/Collide.cpp b/src/pe/collision/Collide.cpp
deleted file mode 100644
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..0000000000000000000000000000000000000000
diff --git a/src/pe/communication/DynamicMarshalling.h b/src/pe/communication/DynamicMarshalling.h
index 121931250e48a482583ca3994f443d51bcb5489d..73f222c15c55051d33474390cc6df2512193f2f6 100644
--- a/src/pe/communication/DynamicMarshalling.h
+++ b/src/pe/communication/DynamicMarshalling.h
@@ -33,6 +33,7 @@
 #include "pe/communication/rigidbody/Capsule.h"
 #include "pe/communication/rigidbody/Sphere.h"
 #include "pe/communication/rigidbody/Union.h"
+#include "pe/utility/BodyCast.h"
 
 #include "core/Abort.h"
 
@@ -44,7 +45,19 @@ namespace pe {
 namespace communication {
 
 template < typename BodyTypeTuple >
-struct MarshalDynamically{
+class MarshalDynamically{
+private:
+   struct MarshalFunctor
+   {
+      mpi::SendBuffer& buffer_;
+
+      MarshalFunctor(mpi::SendBuffer& buffer) : buffer_(buffer) {}
+
+      template< typename BodyType >
+      void operator()( const BodyType* bd) { marshal( buffer_, *bd); }
+   };
+
+public:
    //*************************************************************************************************
    /*!\brief Marshalling rigid body parameters dynamically.
     *
@@ -57,29 +70,30 @@ struct MarshalDynamically{
     */
    static void execute(mpi::SendBuffer& buffer, const RigidBody& b)
    {
-      static_assert(boost::is_base_of<RigidBody, typename BodyTypeTuple::head_type>::value, "only derived types from RigidBody are allowed!");
-
-      if (BodyTypeTuple::head_type::getStaticTypeID() == b.getTypeID())
-      {
-         typedef typename BodyTypeTuple::head_type BodyT;
-         marshal( buffer, static_cast<const BodyT&>(b));
-      } else
-      {
-        MarshalDynamically<typename BodyTypeTuple::tail_type>::execute(buffer, b);
-      }
+      MarshalFunctor func(buffer);
+      return SingleCast<BodyTypeTuple, MarshalFunctor, void>::execute (&b, func);
    }
 };
 
-template < >
-struct MarshalDynamically< boost::tuples::null_type>{
-    static void execute(mpi::SendBuffer& /*buffer*/, const RigidBody& b)
-    {
-       WALBERLA_ABORT("MarshalDynamically: BodyTypeID " << b.getTypeID() << " could not be resolved!");
-    }
-};
-
 template < typename BodyTypeTuple >
-struct UnmarshalDynamically{
+class UnmarshalDynamically{
+private:
+   struct UnmarshalFunctor
+   {
+      mpi::RecvBuffer& buffer_;
+      const math::AABB& domain_;
+      const math::AABB& block_;
+
+      UnmarshalFunctor(mpi::RecvBuffer& buffer, const math::AABB& domain, const math::AABB& block)
+         : buffer_(buffer)
+         , domain_(domain)
+         , block_(block) {}
+
+      template< typename BodyType >
+      BodyID operator()( BodyType* bd) { return instantiate( buffer_, domain_, block_, bd ); }
+   };
+
+public:
    //*************************************************************************************************
    /*!\brief Marshalling rigid body parameters dynamically.
     *
@@ -92,26 +106,11 @@ struct UnmarshalDynamically{
     */
    static BodyID execute(mpi::RecvBuffer& buffer, const id_t typeID, const math::AABB& domain, const math::AABB& block)
    {
-      if (BodyTypeTuple::head_type::getStaticTypeID() == typeID)
-      {
-         typedef typename BodyTypeTuple::head_type BodyT;
-         BodyT* newBody;
-         return instantiate( buffer, domain, block, newBody );
-      } else
-      {
-         return UnmarshalDynamically<typename BodyTypeTuple::tail_type>::execute(buffer, typeID, domain, block);
-      }
+      UnmarshalFunctor func(buffer, domain, block);
+      return SingleCast<BodyTypeTuple, UnmarshalFunctor, BodyID>::execute (typeID, func);
    }
 };
 
-template < >
-struct UnmarshalDynamically< boost::tuples::null_type>{
-    static BodyID execute(mpi::RecvBuffer& /*buffer*/, const id_t typeID, const math::AABB& /*domain*/, const math::AABB& /*block*/)
-    {
-       WALBERLA_ABORT("UnmarshalDynamically: BodyTypeID " << typeID << " could not be resolved!");
-    }
-};
-
 }  // namespace communication
 }  // namespace pe
 }  // namespace walberla
diff --git a/src/pe/communication/rigidbody/Union.h b/src/pe/communication/rigidbody/Union.h
index 0850d0231cc604b4e7f5a5c9a3caec059b819adf..69c2d6bb172ac8589b56f68e7ecc2946ac0b471f 100644
--- a/src/pe/communication/rigidbody/Union.h
+++ b/src/pe/communication/rigidbody/Union.h
@@ -37,11 +37,11 @@ namespace communication {
 
 //forward declaration
 template < typename BodyTypeTuple >
-struct MarshalDynamically;
+class MarshalDynamically;
 
 //forward declaration
 template < typename BodyTypeTuple >
-struct UnmarshalDynamically;
+class UnmarshalDynamically;
 
 struct UnionParameters : public GeomPrimitiveParameters {
    real_t         m_;
diff --git a/src/pe/collision/Collide.h b/src/pe/fcd/AnalyticCollisionDetection.h
similarity index 95%
rename from src/pe/collision/Collide.h
rename to src/pe/fcd/AnalyticCollisionDetection.h
index 1baf317ccdfe314e18c8f01a275430de7acbbbe8..4777794366f6d6cd9eda691b45b187785ed1da01 100644
--- a/src/pe/collision/Collide.h
+++ b/src/pe/fcd/AnalyticCollisionDetection.h
@@ -23,14 +23,13 @@
 #pragma once
 
 #include "pe/Types.h"
+#include "pe/contact/Contact.h"
 #include "pe/rigidbody/Box.h"
 #include "pe/rigidbody/Capsule.h"
 #include "pe/rigidbody/Plane.h"
 #include "pe/rigidbody/Sphere.h"
 #include "pe/rigidbody/Union.h"
-#include "pe/contact/Contact.h"
-
-#include "GJKEPAHelper.h"
+#include "pe/utility/BodyCast.h"
 
 #include "core/debug/Debug.h"
 #include "core/math/RotationMatrix.h"
@@ -41,12 +40,107 @@
 
 namespace walberla {
 namespace pe {
-
 namespace fcd {
-//Forward Declaration
-template < typename BodyTypeTuple, typename BaseT >
-struct DoubleDispatch;
-}
+namespace analytic {
+
+template <typename Container>
+inline
+bool collide( GeomID bd1, GeomID bd2, Container& container );
+
+template <typename Container>
+inline
+bool collide( SphereID s1, SphereID s2, Container& container );
+
+template <typename Container>
+inline
+bool collide( SphereID s, PlaneID p, Container& container );
+
+template <typename Container>
+inline
+bool collide( PlaneID p, SphereID s, Container& container );
+template <typename Container>
+inline
+bool collide( SphereID s, BoxID b, Container& container );
+
+template <typename Container>
+inline
+bool collide( BoxID b, SphereID s, Container& container );
+template <typename Container>
+inline
+bool collide( BoxID b1, BoxID b2, Container& container );
+template <typename Container>
+inline
+bool collide( BoxID b, PlaneID p, Container& container );
+
+template <typename Container>
+inline
+bool collide( PlaneID p, BoxID b, Container& container );
+template <typename Container>
+inline
+bool collide( CapsuleID c1, CapsuleID c2, Container& container );
+template <typename Container>
+inline
+bool collide( CapsuleID c, PlaneID p, Container& container );
+
+template <typename Container>
+inline
+bool collide( PlaneID p, CapsuleID c, Container& container );
+
+template <typename Container>
+inline
+bool collide( SphereID s, CapsuleID c, Container& container );
+
+template <typename Container>
+inline
+bool collide( CapsuleID c, SphereID s, Container& container );
+
+
+template <typename Container>
+inline
+bool collide( BoxID b, CapsuleID c, Container& container );
+
+template <typename Container>
+inline
+bool collide( CapsuleID c, BoxID b, Container& container );
+
+template <typename BodyTypeTuple, typename BodyB, typename Container>
+inline
+bool collide( Union<BodyTypeTuple>* bd1, BodyB* bd2, Container& container );
+
+template <typename BodyA, typename BodyTypeTuple, typename Container>
+inline
+bool collide( BodyA* bd1, Union<BodyTypeTuple>* bd2, Container& container );
+
+template <typename BodyTypeTupleA, typename BodyTypeTupleB, typename Container>
+inline
+bool collide( Union<BodyTypeTupleA>* bd1, Union<BodyTypeTupleB>* bd2, Container& container );
+
+} //namespace analytic
+
+template <typename Container>
+struct AnalyticCollideFunctor
+{
+   Container& contacts_;
+
+   AnalyticCollideFunctor(Container& contacts) : contacts_(contacts) {}
+
+   template< typename BodyType1, typename BodyType2 >
+   bool operator()( BodyType1* bd1, BodyType2* bd2) { using namespace analytic; return collide( bd1, bd2, contacts_); }
+};
+
+template <typename BodyType1, typename Container>
+struct AnalyticSingleCollideFunctor
+{
+   BodyType1* bd1_;
+   Container& contacts_;
+
+   AnalyticSingleCollideFunctor(BodyType1* bd1, Container& contacts) : bd1_(bd1), contacts_(contacts) {}
+
+   template< typename BodyType2 >
+   bool operator()( BodyType2* bd2) { using namespace analytic; return collide( bd1_, bd2, contacts_); }
+};
+
+namespace analytic {
 
 //*************************************************************************************************
 /*!\brief Contact generation between two colliding rigid bodies.
@@ -60,21 +154,9 @@ struct DoubleDispatch;
  */
 template <typename Container>
 inline
-bool collide( GeomID bd1, GeomID bd2, Container& container )
+bool collide( GeomID /*bd1*/, GeomID /*bd2*/, Container& /*container*/ )
 {
    WALBERLA_ABORT("UNSUPPORTED COLLISION!");
-   Vec3   contactPoint;
-   Vec3   contactNormal;
-   real_t penetrationDepth;
-
-   bool collision = collideGJK(bd1, bd2, contactPoint, contactNormal, penetrationDepth);
-
-   if (collision)
-   {
-      container.push_back( Contact(bd1, bd2, contactPoint, contactNormal, penetrationDepth) );
-      return true;
-   }
-
    return false;
 }
 
@@ -2002,10 +2084,11 @@ template <typename BodyTypeTuple, typename BodyB, typename Container>
 inline
 bool collide( Union<BodyTypeTuple>* bd1, BodyB* bd2, Container& container )
 {
+   AnalyticSingleCollideFunctor<BodyB, Container> func(bd2, container);
    bool collision = false;
    for( auto it=bd1->begin(); it!=bd1->end(); ++it )
    {
-      collision |= fcd::DoubleDispatch< BodyTypeTuple, boost::tuple<BodyB> >::execute(*it, bd2, container);
+      collision |= SingleCast<BodyTypeTuple, AnalyticSingleCollideFunctor<BodyB, Container>, bool>::execute(*it, func);
    }
    return collision;
 }
@@ -2021,16 +2104,19 @@ template <typename BodyTypeTupleA, typename BodyTypeTupleB, typename Container>
 inline
 bool collide( Union<BodyTypeTupleA>* bd1, Union<BodyTypeTupleB>* bd2, Container& container )
 {
+   AnalyticCollideFunctor<Container> func(container);
    bool collision = false;
    for( auto it1=bd1->begin(); it1!=bd1->end(); ++it1 )
    {
       for( auto it2=bd2->begin(); it2!=bd2->end(); ++it2 )
       {
-         collision |= fcd::DoubleDispatch< BodyTypeTupleA, BodyTypeTupleB >::execute(*it1, *it2, container);
+         collision |= DoubleCast<BodyTypeTupleA, BodyTypeTupleB, AnalyticCollideFunctor<Container>, bool>::execute(*it1, *it2, func);
       }
    }
    return collision;
 }
 
-}
-}
+} //namespace analytic
+} //namespace fcd
+} //namespace pe
+} //namespace walberla
diff --git a/src/pe/fcd/GenericFCD.h b/src/pe/fcd/GenericFCD.h
new file mode 100644
index 0000000000000000000000000000000000000000..27a311f3c3bdb48552103bed3c47c18f7d05e0d0
--- /dev/null
+++ b/src/pe/fcd/GenericFCD.h
@@ -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 GenericFCD.h
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "IFCD.h"
+
+#include "pe/utility/BodyCast.h"
+
+#include "blockforest/BlockDataHandling.h"
+
+namespace walberla{
+namespace pe{
+namespace fcd {
+
+///
+/// \brief Uses CollideFunctor to call collide function without additional namespace inclusion.
+///
+template <typename BodyTypeTuple, template <typename Container> class CollisionFunctor >
+class GenericFCD : public IFCD{
+public:
+   virtual Contacts& generateContacts(PossibleContacts& possibleContacts)
+   {
+      contacts_.clear();
+      CollisionFunctor<decltype(contacts_)> func(contacts_);
+      for (auto it = possibleContacts.begin(); it != possibleContacts.end(); ++it)
+      {
+         DoubleCast<BodyTypeTuple, BodyTypeTuple, CollisionFunctor<decltype(contacts_)>, bool>::execute(it->first, it->second, func);
+      }
+      return contacts_;
+   }
+};
+
+template <typename BodyTypeTuple, template <typename Container> class CollisionFunctor>
+shared_ptr< blockforest::AlwaysCreateBlockDataHandling<GenericFCD<BodyTypeTuple, CollisionFunctor> > > createGenericFCDDataHandling()
+{
+   return make_shared< blockforest::AlwaysCreateBlockDataHandling<GenericFCD<BodyTypeTuple, CollisionFunctor> > >( );
+}
+
+}
+}
+}
diff --git a/src/pe/fcd/SimpleFCD.h b/src/pe/fcd/SimpleFCD.h
index 456d9b73b41d22b69f83ba5bcf586d5ee5163883..ba5e45316692dc4929823c1af10833b7b60287f9 100644
--- a/src/pe/fcd/SimpleFCD.h
+++ b/src/pe/fcd/SimpleFCD.h
@@ -20,9 +20,8 @@
 
 #pragma once
 
-#include "IFCD.h"
-
-#include "pe/collision/Collide.h"
+#include "AnalyticCollisionDetection.h"
+#include "GenericFCD.h"
 
 #include <boost/type_traits/is_base_of.hpp>
 #include <boost/tuple/tuple.hpp>
@@ -31,68 +30,10 @@ namespace walberla{
 namespace pe{
 namespace fcd {
 
-template < typename TypeA, typename TypeListB >
-struct SingleDispatch{
-   static bool execute(RigidBody* a, RigidBody* b, Contacts& contacts){
-      static_assert(boost::is_base_of<RigidBody, typename TypeListB::head_type>::value, "only downcasting allowed!");
-      if (TypeListB::head_type::getStaticTypeID() == b->getTypeID())
-      {
-         typedef typename TypeListB::head_type TypeB;
-
-         auto bd1 = static_cast<TypeA *>(a);
-         auto bd2 = static_cast<TypeB *>(b);
-
-         return collide(bd1, bd2, contacts);
-      } else
-      {
-         return SingleDispatch<TypeA, typename TypeListB::tail_type>::execute(a, b, contacts);
-      }
-   }
-};
-
-template < typename TypeA >
-struct SingleDispatch< TypeA, boost::tuples::null_type >{
-   static bool execute(RigidBody* /*a*/, RigidBody* b, Contacts& /*contacts*/){
-      WALBERLA_ABORT("SingleDispatch: Type of body " << b->getSystemID() << " could not be determined (" << b->getTypeID() << ")");
-      return false;
-   }
-};
-
-template < typename TypeListA, typename TypeListB = TypeListA>
-struct DoubleDispatch{
-   static bool execute(RigidBody* a, RigidBody* b, Contacts& contacts){
-      // Force a defined order of collision detection across processes
-      if( b->getSystemID() < a->getSystemID() )
-         std::swap( a, b );
-      static_assert(boost::is_base_of<RigidBody, typename TypeListA::head_type>::value, "only downcasting allowed!");
-      if (TypeListA::head_type::getStaticTypeID() == a->getTypeID()) {
-         return SingleDispatch<typename TypeListA::head_type, TypeListB>::execute(a, b, contacts);
-      } else {
-         return DoubleDispatch<typename TypeListA::tail_type, TypeListB>::execute(a, b, contacts);
-      }
-   }
-};
-
-template < typename TypeListB>
-struct DoubleDispatch< boost::tuples::null_type, TypeListB>{
-   static bool execute(RigidBody* /*a*/, RigidBody* /*b*/, Contacts& /*contacts*/){
-      WALBERLA_ABORT("DoubleDispatch: Type could not be determined");
-      return false;
-   }
-};
-
 template <typename BodyTypeTuple>
-class SimpleFCD : public IFCD{
+class SimpleFCD : public GenericFCD<BodyTypeTuple, AnalyticCollideFunctor>{
 public:
-   virtual Contacts& generateContacts(PossibleContacts& possibleContacts)
-   {
-      contacts_.clear();
-      for (auto it = possibleContacts.begin(); it != possibleContacts.end(); ++it)
-      {
-         DoubleDispatch<BodyTypeTuple>::execute(it->first, it->second, contacts_);
-      }
-      return contacts_;
-   }
+   virtual ~SimpleFCD(){}
 };
 
 }
diff --git a/src/pe/fcd/SimpleFCDDataHandling.h b/src/pe/fcd/SimpleFCDDataHandling.h
index 4f46eccf86c7b25a35a19231516d0cebe7b2cf21..1fce240e0b3db419f2bb4ede1048d93ec06a22e8 100644
--- a/src/pe/fcd/SimpleFCDDataHandling.h
+++ b/src/pe/fcd/SimpleFCDDataHandling.h
@@ -23,6 +23,7 @@
 #include "SimpleFCD.h"
 
 #include "blockforest/BlockDataHandling.h"
+#include "core/Deprecated.h"
 
 namespace walberla{
 namespace pe{
@@ -34,6 +35,10 @@ public:
     SimpleFCD<BodyTypeTuple> * initialize( IBlock * const /*block*/ ) {return new SimpleFCD<BodyTypeTuple>();}
 };
 
+/// \attention This function is deprecated. Use createGenericFCDDataHandling<BodyTypeTuple, AnalyticCollideFunctor>() instead!
+template <typename BodyTypeTuple>
+WALBERLA_DEPRECATED(shared_ptr<SimpleFCDDataHandling<BodyTypeTuple> > createSimpleFCDDataHandling());
+
 template <typename BodyTypeTuple>
 shared_ptr<SimpleFCDDataHandling<BodyTypeTuple> > createSimpleFCDDataHandling()
 {
diff --git a/src/pe/utility/BodyCast.h b/src/pe/utility/BodyCast.h
new file mode 100644
index 0000000000000000000000000000000000000000..93c14befb3c6d4d719296855e3b2024832690c31
--- /dev/null
+++ b/src/pe/utility/BodyCast.h
@@ -0,0 +1,162 @@
+//======================================================================================================================
+//
+//  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 BodyCast.h
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include <core/DataTypes.h>
+#include <pe/rigidbody/RigidBody.h>
+
+#include <boost/tuple/tuple.hpp>
+
+namespace walberla {
+namespace pe {
+
+template < typename TypeList, typename Functor, typename ReturnType >
+class SingleCast
+{
+public:
+   static ReturnType execute(const id_t typeID, Functor& func){
+      static_assert(boost::is_base_of<RigidBody, typename TypeList::head_type>::value, "only downcasting allowed!");
+      if (TypeList::head_type::getStaticTypeID() == typeID)
+      {
+         typedef typename TypeList::head_type CastBodyType;
+         CastBodyType* bd = NULL;
+         return func( static_cast<CastBodyType *>( bd ) );
+      } else
+      {
+         return SingleCast<typename TypeList::tail_type, Functor, ReturnType>::execute(typeID, func);
+      }
+   }
+
+   static ReturnType execute(RigidBody* bd, Functor& func){
+      static_assert(boost::is_base_of<RigidBody, typename TypeList::head_type>::value, "only downcasting allowed!");
+      if (TypeList::head_type::getStaticTypeID() == bd->getTypeID())
+      {
+         typedef typename TypeList::head_type CastBodyType;
+         return func( static_cast<CastBodyType *>(bd) );
+      } else
+      {
+         return SingleCast<typename TypeList::tail_type, Functor, ReturnType>::execute(bd, func);
+      }
+   }
+
+   static ReturnType execute(const RigidBody* bd, Functor& func){
+      static_assert(boost::is_base_of<RigidBody, typename TypeList::head_type>::value, "only downcasting allowed!");
+      if (TypeList::head_type::getStaticTypeID() == bd->getTypeID())
+      {
+         typedef typename TypeList::head_type CastBodyType;
+         return func( static_cast<const CastBodyType *>(bd) );
+      } else
+      {
+         return SingleCast<typename TypeList::tail_type, Functor, ReturnType>::execute(bd, func);
+      }
+   }
+};
+
+template < typename Functor, typename ReturnType >
+struct SingleCast< boost::tuples::null_type, Functor, ReturnType >{
+   static ReturnType execute(const id_t typeID, Functor& /*func*/)
+   {
+      WALBERLA_ABORT("SingleCast: BodyType could not be determined (" << typeID << ")");
+      return ReturnType();
+   }
+
+   static ReturnType execute(RigidBody* bd, Functor& /*func*/)
+   {
+      WALBERLA_ABORT("SingleCast: Type of body " << bd->getSystemID() << " could not be determined (" << bd->getTypeID() << ")");
+      return ReturnType();
+   }
+   static ReturnType execute(const RigidBody* bd, Functor& /*func*/)
+   {
+      WALBERLA_ABORT("SingleCast: Type of body " << bd->getSystemID() << " could not be determined (" << bd->getTypeID() << ")");
+      return ReturnType();
+   }
+};
+
+template < typename TypeListA, typename TypeListB, typename Functor, typename ReturnType >
+class DoubleCast
+{
+private:
+   template< typename BodyType1 >
+   struct SingleCastFunctor
+   {
+      BodyType1* a_;
+      Functor& func_;
+
+      SingleCastFunctor(BodyType1* a, Functor& func) : a_(a), func_(func) {}
+
+      template< typename BodyType2 >
+      ReturnType operator()( BodyType2* bd) { return func_( a_, bd); }
+   };
+   template< typename BodyType1 >
+   struct SingleCastConstFunctor
+   {
+      BodyType1 const * a_;
+      Functor& func_;
+
+      SingleCastConstFunctor(BodyType1 const * a, Functor& func) : a_(a), func_(func) {}
+
+      template< typename BodyType2 >
+      ReturnType operator()( BodyType2 const * bd) { return func_( a_, bd); }
+   };
+public:
+   static ReturnType execute(RigidBody* a, RigidBody* b, Functor& func){
+      static_assert(boost::is_base_of<RigidBody, typename TypeListA::head_type>::value, "only downcasting allowed!");
+      if (TypeListA::head_type::getStaticTypeID() == a->getTypeID())
+      {
+         typedef typename TypeListA::head_type CastBodyType;
+         SingleCastFunctor<CastBodyType> singleFunc( static_cast<CastBodyType *>(a), func);
+         return SingleCast<TypeListB, SingleCastFunctor<CastBodyType>, ReturnType>::execute(b, singleFunc );
+      } else
+      {
+         return DoubleCast<typename TypeListA::tail_type, TypeListB, Functor, ReturnType>::execute(a, b, func);
+      }
+   }
+
+   static ReturnType execute(const RigidBody* a, const RigidBody* b, Functor& func){
+      static_assert(boost::is_base_of<RigidBody, typename TypeListA::head_type>::value, "only downcasting allowed!");
+      if (TypeListA::head_type::getStaticTypeID() == a->getTypeID())
+      {
+         typedef typename TypeListA::head_type CastBodyType;
+         SingleCastConstFunctor<CastBodyType> singleFunc( static_cast<CastBodyType *>(a), func);
+         return SingleCast<TypeListB, SingleCastConstFunctor<CastBodyType>, ReturnType>::execute(b, singleFunc );
+      } else
+      {
+         return DoubleCast<typename TypeListA::tail_type, TypeListB, Functor, ReturnType>::execute(a, b, func);
+      }
+   }
+};
+
+template < typename TypeListB, typename Functor, typename ReturnType >
+struct DoubleCast< boost::tuples::null_type, TypeListB, Functor, ReturnType >{
+   static ReturnType execute(RigidBody* a, RigidBody* /*b*/, Functor& /*func*/)
+   {
+      WALBERLA_ABORT("DoubleCast: Type of body " << a->getSystemID() << " could not be determined (" << a->getTypeID() << ")");
+      return ReturnType();
+   }
+   static ReturnType execute(const RigidBody* a, const RigidBody* b, Functor& /*func*/)
+   {
+      WALBERLA_ABORT("DoubleCast: Type of body " << a->getSystemID() << " could not be determined (" << a->getTypeID() << ")");
+      return ReturnType();
+   }
+};
+
+}  // namespace pe
+}  // namespace walberla
diff --git a/tests/pe/BodyFlags.cpp b/tests/pe/BodyFlags.cpp
index 824b43e811f87bc2d809ebb1156cf0c720dcd6c8..901e85c71bf41bb6b0c4238fc8fff6aa6e550db1 100644
--- a/tests/pe/BodyFlags.cpp
+++ b/tests/pe/BodyFlags.cpp
@@ -64,7 +64,7 @@ int main( int argc, char ** argv )
 
    auto storageID           = forest->addBlockData(createStorageDataHandling<BodyTuple>(), "Storage");
    auto ccdID               = forest->addBlockData(ccd::createHashGridsDataHandling( globalStorage, storageID ), "CCD");
-   auto fcdID               = forest->addBlockData(fcd::createSimpleFCDDataHandling<BodyTuple>(), "FCD");
+   auto fcdID               = forest->addBlockData(fcd::createGenericFCDDataHandling<BodyTuple, fcd::AnalyticCollideFunctor>(), "FCD");
 
    //cr::DEM    cr(globalStorage, forest->getBlockStorage(), storageID, ccdID, fcdID, NULL );
    cr::HCSITS cr(globalStorage, forest->getBlockStoragePointer(), storageID, ccdID, fcdID, NULL );
diff --git a/tests/pe/Collision.cpp b/tests/pe/Collision.cpp
index faa5a5f420f7c2f6f259b03f79cf01012c4748f7..9abc2fb2e7080262dda47b15ca7b48c8f18dbd6c 100644
--- a/tests/pe/Collision.cpp
+++ b/tests/pe/Collision.cpp
@@ -18,8 +18,9 @@
 //
 //======================================================================================================================
 
-#include "pe/collision/Collide.h"
 #include "pe/collision/GJKEPAHelper.h"
+#include "pe/fcd/AnalyticCollisionDetection.h"
+#include "pe/utility/BodyCast.h"
 
 #include "pe/contact/Contact.h"
 #include "pe/fcd/SimpleFCD.h"
@@ -40,6 +41,7 @@
 
 using namespace walberla;
 using namespace walberla::pe;
+using walberla::pe::fcd::analytic::collide;
 
 void checkContact(const Contact& c1, const Contact& c2)
 {
@@ -59,11 +61,12 @@ void SphereTest()
    Plane  pl1(223, 1, Vec3(0,0,0), Vec3(1,1,1).getNormalized(), 0, iron);
 
    std::vector<Contact> contacts;
+   fcd::AnalyticCollideFunctor< std::vector<Contact> > collideFunc(contacts);
 
    // SPHERE <-> SPHERE
    WALBERLA_LOG_INFO("SPHERE <-> SPHERE");
-   WALBERLA_CHECK( !collide(&sp1, &sp3, contacts) );
-   WALBERLA_CHECK(  collide(&sp1, &sp2, contacts) );
+   WALBERLA_CHECK( !collideFunc(&sp1, &sp3) );
+   WALBERLA_CHECK(  collideFunc(&sp1, &sp2) );
    checkContact( contacts.at(0),
                  Contact( &sp1, &sp2, Vec3(real_t(0.75), 0, 0), Vec3(-1, 0, 0), real_t(-0.5)) );
 
@@ -71,17 +74,17 @@ void SphereTest()
    // SPHERE <-> PLANE
    WALBERLA_LOG_INFO("SPHERE <-> PLANE");
    contacts.clear();
-   WALBERLA_CHECK(  collide(&sp1, &pl1, contacts) );
-   WALBERLA_CHECK(  collide(&sp2, &pl1, contacts) );
+   WALBERLA_CHECK(  collideFunc(&sp1, &pl1) );
+   WALBERLA_CHECK(  collideFunc(&sp2, &pl1) );
    checkContact( contacts.at(1),
                  Contact( &sp2, &pl1, Vec3(1,real_t(-0.5),real_t(-0.5)), Vec3(1, 1, 1).getNormalized(), real_t(-0.133974596215561181)) );
    // ORDER INVARIANCE
    contacts.clear();
-   WALBERLA_CHECK(  collide(&pl1, &sp2, contacts) );
+   WALBERLA_CHECK(  collideFunc(&pl1, &sp2) );
    checkContact( contacts.at(0),
                  Contact( &sp2, &pl1, Vec3(1,real_t(-0.5),real_t(-0.5)), Vec3(1, 1, 1).getNormalized(), real_t(-0.133974596215561181)) );
 
-   WALBERLA_CHECK( !collide(&sp3, &pl1, contacts) );
+   WALBERLA_CHECK( !collideFunc(&sp3, &pl1) );
 }
 
 void BoxTest()
@@ -98,6 +101,7 @@ void BoxTest()
    b5.rotate( Vec3(1,0,0), real_t(math::PI * 0.25) );
 
    std::vector<Contact> contacts;
+   fcd::AnalyticCollideFunctor< std::vector<Contact> > collideFunc(contacts);
 
    real_t penetrationDepth = real_t(0);
    Vec3   contactPoint = Vec3();
@@ -108,12 +112,12 @@ void BoxTest()
    // BOX <-> BOX
    WALBERLA_LOG_INFO("BOX <-> BOX");
    WALBERLA_CHECK( !collideGJK(&b1, &b3, contactPoint, contactNormal, penetrationDepth) );
-   WALBERLA_CHECK( !collide(&b1, &b3, contacts) );
+   WALBERLA_CHECK( !collideFunc(&b1, &b3) );
 //   WALBERLA_LOG_WARNING("contactPoint    : " << contactPoint);
 //   WALBERLA_LOG_WARNING("contactNormal   : " << contactNormal);
 //   WALBERLA_LOG_WARNING("penetrationDepth: " << penetrationDepth);
    WALBERLA_CHECK(  collideGJK(&b1, &b2, contactPoint, contactNormal, penetrationDepth) );
-   WALBERLA_CHECK(  collide(&b1, &b2, contacts) );
+   WALBERLA_CHECK(  collideFunc(&b1, &b2) );
 //   WALBERLA_LOG_WARNING("contactPoint    : " << contactPoint);
 //   WALBERLA_LOG_WARNING("contactNormal   : " << contactNormal);
 //   WALBERLA_LOG_WARNING("penetrationDepth: " << penetrationDepth);
@@ -121,29 +125,29 @@ void BoxTest()
 
    b4.setPosition( (Vec3(0,0,1) * real_t(sqrt(3)) + Vec3(0,0,1)) * 0.999);
    WALBERLA_CHECK( collideGJK(&b1, &b4, contactPoint, contactNormal, penetrationDepth) );
-   WALBERLA_CHECK( collide(&b1, &b4, contacts) );
+   WALBERLA_CHECK( collideFunc(&b1, &b4) );
 //   WALBERLA_LOG_WARNING("contactPoint    : " << contacts.back().getPosition());
 //   WALBERLA_LOG_WARNING("contactNormal   : " << contacts.back().getNormal());
 //   WALBERLA_LOG_WARNING("penetrationDepth: " << contacts.back().getDistance());
 
    b4.setPosition( (Vec3(0,0,1) * real_t(sqrt(3)) + Vec3(0,0,1)) * 1.001);
    WALBERLA_CHECK( !collideGJK(&b1, &b4, contactPoint, contactNormal, penetrationDepth) );
-   WALBERLA_CHECK( !collide(&b1, &b4, contacts) );
+   WALBERLA_CHECK( !collideFunc(&b1, &b4) );
 
    b5.setPosition( (Vec3(0,0,1) * real_t(sqrt(3)) + Vec3(0,0,1)) * 0.99);
    WALBERLA_CHECK( collideGJK(&b1, &b5, contactPoint, contactNormal, penetrationDepth) );
-   WALBERLA_CHECK( collide(&b1, &b5, contacts) );
+   WALBERLA_CHECK( collideFunc(&b1, &b5) );
 //   WALBERLA_LOG_WARNING("contactPoint    : " << contacts.back().getPosition());
 //   WALBERLA_LOG_WARNING("contactNormal   : " << contacts.back().getNormal());
 //   WALBERLA_LOG_WARNING("penetrationDepth: " << contacts.back().getDistance());
 
    b5.setPosition( (Vec3(0,0,1) * real_t(sqrt(3)) + Vec3(0,0,1)) * 1.01);
    WALBERLA_CHECK( !collideGJK(&b1, &b5, contactPoint, contactNormal, penetrationDepth) );
-   WALBERLA_CHECK( !collide(&b1, &b5, contacts) );
+   WALBERLA_CHECK( !collideFunc(&b1, &b5) );
 
    Sphere s1(126, 0, Vec3(real_t(1.5), real_t(1.5), real_t(1.5)), Vec3(0,0,0), Quat(), 1, iron, false, true, false);
    WALBERLA_CHECK( collideGJK(&b1, &s1, contactPoint, contactNormal, penetrationDepth) );
-   WALBERLA_CHECK( collide(&b1, &s1, contacts) );
+   WALBERLA_CHECK( collideFunc(&b1, &s1) );
 //   WALBERLA_LOG_WARNING("contactPoint    : " << contactPoint);
 //   WALBERLA_LOG_WARNING("contactNormal   : " << contactNormal);
 //   WALBERLA_LOG_WARNING("penetrationDepth: " << penetrationDepth);
@@ -156,6 +160,7 @@ void CapsuleTest()
    Sphere sp1(123, 123, Vec3(0,0,0), Vec3(0,0,0), Quat(), 1, iron, false, true, false);
 
    std::vector<Contact> contacts;
+   fcd::AnalyticCollideFunctor< std::vector<Contact> > collideFunc(contacts);
 
    // CAPSULE <-> SPHERE
    WALBERLA_LOG_INFO("CAPSULE <-> SPHERE");
@@ -168,30 +173,30 @@ void CapsuleTest()
    WALBERLA_CHECK( !collideGJK(&c1, &sp1, contacts) );
 
    sp1.setPosition(0, real_t(1.9), 0);
-   WALBERLA_CHECK( collide(&c1, &sp1, contacts) );
+   WALBERLA_CHECK( collideFunc(&c1, &sp1) );
 //   WALBERLA_LOG_WARNING("contactPoint    : " << contacts.at(1).getPosition());
 //   WALBERLA_LOG_WARNING("contactNormal   : " << contacts.at(1).getNormal());
 //   WALBERLA_LOG_WARNING("penetrationDepth: " << contacts.at(1).getDistance());
    sp1.setPosition(0, real_t(2.1), 0);
-   WALBERLA_CHECK( !collide(&c1, &sp1, contacts) );
+   WALBERLA_CHECK( !collideFunc(&c1, &sp1) );
 
    // CAPSULE <-> PLANE
    WALBERLA_LOG_INFO("CAPSULE <-> PLANE");
    Plane pl1(124, 124, Vec3(0,0,real_t(-0.9)), Vec3(0,0,1), 0, iron);
-   WALBERLA_CHECK( collide(&c1, &pl1, contacts) );
+   WALBERLA_CHECK( collideFunc(&c1, &pl1) );
 //   WALBERLA_LOG_WARNING("contactPoint    : " << contacts.at(2).getPosition());
 //   WALBERLA_LOG_WARNING("contactNormal   : " << contacts.at(2).getNormal());
 //   WALBERLA_LOG_WARNING("penetrationDepth: " << contacts.at(2).getDistance());
    pl1.setPosition(0,0, real_t(-1.1));
-   WALBERLA_CHECK( !collide(&c1, &pl1, contacts) );
+   WALBERLA_CHECK( !collideFunc(&c1, &pl1) );
 
    Plane pl2(124, 124, Vec3(real_t(1.9),0,0), Vec3(-1,0,0), 0, iron);
-   WALBERLA_CHECK( collide(&c1, &pl2, contacts) );
+   WALBERLA_CHECK( collideFunc(&c1, &pl2) );
 //   WALBERLA_LOG_WARNING("contactPoint    : " << contacts.at(3).getPosition());
 //   WALBERLA_LOG_WARNING("contactNormal   : " << contacts.at(3).getNormal());
 //   WALBERLA_LOG_WARNING("penetrationDepth: " << contacts.at(3).getDistance());
    pl2.setPosition(real_t(2.1),0, 0);
-   WALBERLA_CHECK( !collide(&c1, &pl2, contacts) );
+   WALBERLA_CHECK( !collideFunc(&c1, &pl2) );
 
 }
 
@@ -256,9 +261,12 @@ void UnionTest()
 
    std::vector<Contact> contacts;
 
+   using namespace walberla::pe::fcd;
    // SPHERE <-> SPHERE
    WALBERLA_LOG_INFO("UNION <-> UNION");
-   fcd::DoubleDispatch< boost::tuple<UnionT, Sphere> >::execute(&un1, &un2, contacts);
+   AnalyticCollideFunctor< std::vector<Contact> > func(contacts);
+   func(&un1, &un2);
+
    checkContact( contacts.at(0),
                  Contact( sp1, sp2, Vec3(real_t(0.75), 0, 0), Vec3(-1, 0, 0), real_t(-0.5)) );
 }
diff --git a/tests/pe/HCSITS.cpp b/tests/pe/HCSITS.cpp
index 93faeca367767d043da284410a51fef9a9acad8e..ceda7017e1ce51f0a0ccddb08081e2480aa76dcc 100644
--- a/tests/pe/HCSITS.cpp
+++ b/tests/pe/HCSITS.cpp
@@ -116,7 +116,7 @@ int main( int argc, char** argv )
 
    auto storageID           = forest->addBlockData(createStorageDataHandling<BodyTuple>(), "Storage");
    auto hccdID              = forest->addBlockData(ccd::createHashGridsDataHandling( globalBodyStorage, storageID ), "HCCD");
-   auto fcdID               = forest->addBlockData(fcd::createSimpleFCDDataHandling<BodyTuple>(), "FCD");
+   auto fcdID               = forest->addBlockData(fcd::createGenericFCDDataHandling<BodyTuple, fcd::AnalyticCollideFunctor>(), "FCD");
    cr::HCSITS cr(globalBodyStorage, forest->getBlockStoragePointer(), storageID, hccdID, fcdID);
    cr.setMaxIterations( 10 );
    cr.setRelaxationParameter    ( real_t(0.7) );
diff --git a/tests/pe/HashGrids.cpp b/tests/pe/HashGrids.cpp
index 1a01e29c6fbe01d8d30a080bf05f1cfa5f43a24d..cc48c5b038b6f6034449bcf523ae4ab7591bce24 100644
--- a/tests/pe/HashGrids.cpp
+++ b/tests/pe/HashGrids.cpp
@@ -65,7 +65,7 @@ int main( int argc, char** argv )
     auto storageID           = forest->addBlockData(createStorageDataHandling<BodyTuple>(), "Storage");
     auto sccdID              = forest->addBlockData(ccd::createSimpleCCDDataHandling( globalBodyStorage, storageID ), "SCCD");
     auto hccdID              = forest->addBlockData(ccd::createHashGridsDataHandling( globalBodyStorage, storageID ), "HCCD");
-    auto fcdID               = forest->addBlockData(fcd::createSimpleFCDDataHandling<BodyTuple>(), "FCD");
+    auto fcdID               = forest->addBlockData(fcd::createGenericFCDDataHandling<BodyTuple, fcd::AnalyticCollideFunctor>(), "FCD");
     cr::PlainIntegrator cr(globalBodyStorage, forest->getBlockStoragePointer(), storageID, NULL);
 
     pe::createPlane( *globalBodyStorage, 0, Vec3(0, +1, 0), Vec3(5, 0,5), iron);
diff --git a/tests/pe/ParallelEquivalence.cpp b/tests/pe/ParallelEquivalence.cpp
index 0448c509ee0a370b07947bba1f3090c64f9c2b6d..6a211cd6739ad9613e682832b7269d9d5184531a 100644
--- a/tests/pe/ParallelEquivalence.cpp
+++ b/tests/pe/ParallelEquivalence.cpp
@@ -103,7 +103,7 @@ void sim(shared_ptr< StructuredBlockForest > forest, std::vector<BodyData>& res,
    {
       ccdID_                 = forest->addBlockData(ccd::createSimpleCCDDataHandling( globalStorage, storageID_ ), "CCD");
    }
-   auto fcdID_               = forest->addBlockData(fcd::createSimpleFCDDataHandling<BodyTuple>(), "FCD");
+   auto fcdID_               = forest->addBlockData(fcd::createGenericFCDDataHandling<BodyTuple, fcd::AnalyticCollideFunctor>(), "FCD");
 
    cr::DEM cr( globalStorage, forest->getBlockStoragePointer(), storageID_, ccdID_, fcdID_, NULL);
 
diff --git a/tests/pe/PeDocumentationSnippets.cpp b/tests/pe/PeDocumentationSnippets.cpp
index 109f23d7cc477bbda2b79777bdef7102b46ab89a..4cf46e2f8b0be31450f3900261a9e5bdcd53650d 100644
--- a/tests/pe/PeDocumentationSnippets.cpp
+++ b/tests/pe/PeDocumentationSnippets.cpp
@@ -79,7 +79,7 @@ int main( int argc, char ** argv )
 
    auto storageID           = forest->addBlockData(createStorageDataHandling<BodyTypeTuple>(), "Storage");
    auto ccdID               = forest->addBlockData(ccd::createHashGridsDataHandling( globalBodyStorage, storageID ), "CCD");
-   auto fcdID               = forest->addBlockData(fcd::createSimpleFCDDataHandling<BodyTypeTuple>(), "FCD");
+   auto fcdID               = forest->addBlockData(fcd::createGenericFCDDataHandling<BodyTypeTuple, fcd::AnalyticCollideFunctor>(), "FCD");
 
    //! [Setup HCSITS]
    cr::HCSITS hcsits(globalBodyStorage, forest->getBlockStoragePointer(), storageID, ccdID, fcdID);
diff --git a/tests/pe/Refinement.cpp b/tests/pe/Refinement.cpp
index cb8291d28632befc0e92a15ad47ec38f8258e8f6..cb64385b6b755302933c2c5a50d4988d66f1f5e9 100644
--- a/tests/pe/Refinement.cpp
+++ b/tests/pe/Refinement.cpp
@@ -138,7 +138,7 @@ int main( int argc, char ** argv )
 
    auto storageID           = forest->addBlockData(createStorageDataHandling<BodyTuple>(), "Storage");
    auto ccdID               = forest->addBlockData(ccd::createHashGridsDataHandling( globalStorage, storageID ), "CCD");
-   auto fcdID               = forest->addBlockData(fcd::createSimpleFCDDataHandling<BodyTuple>(), "FCD");
+   auto fcdID               = forest->addBlockData(fcd::createGenericFCDDataHandling<BodyTuple, fcd::AnalyticCollideFunctor>(), "FCD");
 
    //cr::DEM    cr(globalStorage, forest->getBlockStorage(), storageID, ccdID, fcdID, NULL );
    cr::HCSITS cr(globalStorage, forest->getBlockStoragePointer(), storageID, ccdID, fcdID, NULL );
diff --git a/tests/pe/SetBodyTypeIDs.cpp b/tests/pe/SetBodyTypeIDs.cpp
index dee15759db4c51d2f8c9ee88d66a2bfde8af081b..9e99059b24d2fb962f788ff9697130c836dbcf57 100644
--- a/tests/pe/SetBodyTypeIDs.cpp
+++ b/tests/pe/SetBodyTypeIDs.cpp
@@ -20,6 +20,8 @@
 
 #include "core/debug/TestSubsystem.h"
 
+#include "pe/Materials.h"
+
 #include "pe/rigidbody/SetBodyTypeIDs.h"
 
 #include "pe/rigidbody/Box.h"
@@ -28,6 +30,8 @@
 #include "pe/rigidbody/Sphere.h"
 #include "pe/rigidbody/Union.h"
 
+#include "pe/utility/BodyCast.h"
+
 namespace walberla {
 namespace pe {
 
@@ -85,6 +89,25 @@ id_t C::staticTypeID_ = 100;
 
 using namespace walberla::pe;
 
+struct SingleTypeFunctor
+{
+   std::vector<walberla::id_t> ids_;
+
+   void operator()( const BodyID /*bd*/) { }
+   void operator()( const SphereID /*bd*/) { ids_.push_back(Sphere::getStaticTypeID()); }
+   void operator()( const BoxID /*bd*/) { ids_.push_back(Box::getStaticTypeID()); }
+   void operator()( const CapsuleID /*bd*/) { ids_.push_back(Capsule::getStaticTypeID()); }
+};
+
+struct DoubleTypeFunctor
+{
+   std::vector<walberla::id_t> ids_;
+
+   void operator()( const BodyID /*bd1*/, const BodyID /*bd2*/) { }
+   void operator()( const BoxID /*bd1*/, const CapsuleID /*bd2*/) { ids_.push_back(Box::getStaticTypeID()); ids_.push_back(Capsule::getStaticTypeID()); }
+   void operator()( const CapsuleID /*bd1*/, const CapsuleID /*bd2*/) { ids_.push_back(5); }
+};
+
 int main( int argc, char** argv )
 {
    walberla::debug::enterTestMode();
@@ -108,5 +131,30 @@ int main( int argc, char** argv )
    WALBERLA_CHECK_UNEQUAL(Plane::getStaticTypeID(), Box::getStaticTypeID());
    WALBERLA_CHECK_UNEQUAL(Plane::getStaticTypeID(), Capsule::getStaticTypeID());
 
+   SingleTypeFunctor singleFunc;
+   Box     bx (0, 0, Vec3(0), Vec3(0), Quat(), Vec3(1), Material::find("iron"), false, false, false);
+   Capsule cap(0, 0, Vec3(0), Vec3(0), Quat(), 1, 1, Material::find("iron"), false, false, false);
+
+   SingleCast<BodyTuple2, SingleTypeFunctor, void>::execute(Box::getStaticTypeID(), singleFunc);
+   SingleCast<BodyTuple2, SingleTypeFunctor, void>::execute(Capsule::getStaticTypeID(), singleFunc);
+   WALBERLA_CHECK_EQUAL( singleFunc.ids_.size(), 2);
+   WALBERLA_CHECK_EQUAL( singleFunc.ids_[0], Box::getStaticTypeID());
+   WALBERLA_CHECK_EQUAL( singleFunc.ids_[1], Capsule::getStaticTypeID());
+
+   singleFunc.ids_.clear();
+   SingleCast<BodyTuple2, SingleTypeFunctor, void>::execute(&bx, singleFunc);
+   SingleCast<BodyTuple2, SingleTypeFunctor, void>::execute(&cap, singleFunc);
+   WALBERLA_CHECK_EQUAL( singleFunc.ids_.size(), 2);
+   WALBERLA_CHECK_EQUAL( singleFunc.ids_[0], Box::getStaticTypeID());
+   WALBERLA_CHECK_EQUAL( singleFunc.ids_[1], Capsule::getStaticTypeID());
+
+   DoubleTypeFunctor doubleFunc;
+   DoubleCast<BodyTuple2, BodyTuple2, DoubleTypeFunctor, void>::execute(&bx, &cap, doubleFunc);
+   DoubleCast<BodyTuple2, BodyTuple2, DoubleTypeFunctor, void>::execute(&cap, &cap, doubleFunc);
+   WALBERLA_CHECK_EQUAL( doubleFunc.ids_.size(), 3);
+   WALBERLA_CHECK_EQUAL( doubleFunc.ids_[0], Box::getStaticTypeID());
+   WALBERLA_CHECK_EQUAL( doubleFunc.ids_[1], Capsule::getStaticTypeID());;
+   WALBERLA_CHECK_EQUAL( doubleFunc.ids_[2], 5);;
+
    return EXIT_SUCCESS;
 }
diff --git a/tests/pe/ShadowCopy.cpp b/tests/pe/ShadowCopy.cpp
index 1661343557cd7a9018241c07c85d25637b13d928..d0d3c119fd47b8c7f3a2ccef15a08f86ab935870 100644
--- a/tests/pe/ShadowCopy.cpp
+++ b/tests/pe/ShadowCopy.cpp
@@ -54,7 +54,7 @@ int main( int argc, char** argv )
 
    auto storageID           = forest->addBlockData(createStorageDataHandling<BodyTuple>(), "Storage");
                               forest->addBlockData(ccd::createHashGridsDataHandling( globalBodyStorage, storageID ), "HCCD");
-                              forest->addBlockData(fcd::createSimpleFCDDataHandling<BodyTuple>(), "FCD");
+                              forest->addBlockData(fcd::createGenericFCDDataHandling<BodyTuple, fcd::AnalyticCollideFunctor>(), "FCD");
 
 
 //   logging::Logging::instance()->setStreamLogLevel(logging::Logging::DETAIL);
diff --git a/tests/pe/SyncEquivalence.cpp b/tests/pe/SyncEquivalence.cpp
index 728268de134ef3e58d74fbfcaf3c3756fbb02d1c..9d5cb83b52bfeb5e2c9a977e8a9f191215dc9827 100644
--- a/tests/pe/SyncEquivalence.cpp
+++ b/tests/pe/SyncEquivalence.cpp
@@ -113,7 +113,7 @@ void createSimulation(math::AABB& simulationDomain,
     // add block data
     info.storageID           = info.forest->addBlockData(createStorageDataHandling<BodyTuple>(), "Storage");
     info.ccdID               = info.forest->addBlockData(ccd::createHashGridsDataHandling( info.globalBodyStorage, info.storageID ), "CCD");
-    info.fcdID               = info.forest->addBlockData(fcd::createSimpleFCDDataHandling<BodyTuple>(), "FCD");
+    info.fcdID               = info.forest->addBlockData(fcd::createGenericFCDDataHandling<BodyTuple, fcd::AnalyticCollideFunctor>(), "FCD");
 
     info.cr = shared_ptr<cr::ICR>(new cr::HCSITS(info.globalBodyStorage, info.forest->getBlockStoragePointer(), info.storageID, info.ccdID, info.fcdID) );
 
diff --git a/tests/pe/Union.cpp b/tests/pe/Union.cpp
index b451b128841267d75fb3317116152b34c101f24e..427e3752eda06d13e335e9448bcd0b0ab56695a4 100644
--- a/tests/pe/Union.cpp
+++ b/tests/pe/Union.cpp
@@ -68,7 +68,7 @@ int main( int argc, char ** argv )
 
    auto storageID           = forest->addBlockData(createStorageDataHandling<BodyTuple>(), "Storage");
    auto ccdID               = forest->addBlockData(ccd::createHashGridsDataHandling( globalBodyStorage, storageID ), "CCD");
-   auto fcdID               = forest->addBlockData(fcd::createSimpleFCDDataHandling<BodyTuple>(), "FCD");
+   auto fcdID               = forest->addBlockData(fcd::createGenericFCDDataHandling<BodyTuple, fcd::AnalyticCollideFunctor>(), "FCD");
 
    cr::HCSITS cr(globalBodyStorage, forest->getBlockStoragePointer(), storageID, ccdID, fcdID);
    cr.setMaxIterations( uint_c(10) );
diff --git a/tests/pe_coupling/momentum_exchange_method/BodyAtBlockBoarderCheck.cpp b/tests/pe_coupling/momentum_exchange_method/BodyAtBlockBoarderCheck.cpp
index a8922a4813a8bd5273a1636c21527979e2367657..03cd4c395ae7b52c397c623b21cddf8aecb8f673 100644
--- a/tests/pe_coupling/momentum_exchange_method/BodyAtBlockBoarderCheck.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/BodyAtBlockBoarderCheck.cpp
@@ -304,7 +304,7 @@ int main( int argc, char **argv )
    shared_ptr<pe::BodyStorage> globalBodyStorage = make_shared<pe::BodyStorage>();
    auto bodyStorageID = blocks->addBlockData(pe::createStorageDataHandling<BodyTypeTuple>(), "Storage");
    auto ccdID         = blocks->addBlockData(pe::ccd::createHashGridsDataHandling( globalBodyStorage, bodyStorageID ), "CCD");
-   auto fcdID         = blocks->addBlockData(pe::fcd::createSimpleFCDDataHandling<BodyTypeTuple>(), "FCD");
+   auto fcdID         = blocks->addBlockData(pe::fcd::createGenericFCDDataHandling<BodyTypeTuple, pe::fcd::AnalyticCollideFunctor>(), "FCD");
    pe::cr::HCSITS cr(globalBodyStorage, blocks->getBlockStoragePointer(), bodyStorageID, ccdID, fcdID);
 
 
diff --git a/tests/pe_coupling/momentum_exchange_method/LubricationCorrectionMEM.cpp b/tests/pe_coupling/momentum_exchange_method/LubricationCorrectionMEM.cpp
index 64145d3dde11640b1d397c2039099e7f6c9748ff..d782e03fe8c8608734193588099e22189bf0a974 100644
--- a/tests/pe_coupling/momentum_exchange_method/LubricationCorrectionMEM.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/LubricationCorrectionMEM.cpp
@@ -783,7 +783,7 @@ int main( int argc, char **argv )
 
    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");
+   auto fcdID          = blocks->addBlockData(pe::fcd::createGenericFCDDataHandling<BodyTypeTuple, pe::fcd::AnalyticCollideFunctor>(), "FCD");
 
    // set up collision response, here DEM solver
    pe::cr::DEM cr(globalBodyStorage, blocks->getBlockStoragePointer(), bodyStorageID, ccdID, fcdID, NULL);
diff --git a/tests/pe_coupling/momentum_exchange_method/PeriodicParticleChannelMEM.cpp b/tests/pe_coupling/momentum_exchange_method/PeriodicParticleChannelMEM.cpp
index 2fda5409eeec90387b16940621dbb67c5da91ef5..51ae2623cf6413f54c9efa451cebcc94b7f64e18 100644
--- a/tests/pe_coupling/momentum_exchange_method/PeriodicParticleChannelMEM.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/PeriodicParticleChannelMEM.cpp
@@ -411,7 +411,7 @@ int main( int argc, char **argv )
 
    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");
+   auto fcdID          = blocks->addBlockData(pe::fcd::createGenericFCDDataHandling<BodyTypeTuple, pe::fcd::AnalyticCollideFunctor>(), "FCD");
 
    // set up collision response, here DEM solver
    pe::cr::DEM cr(globalBodyStorage, blocks->getBlockStoragePointer(), bodyStorageID, ccdID, fcdID, NULL);
diff --git a/tests/pe_coupling/momentum_exchange_method/SegreSilberbergMEM.cpp b/tests/pe_coupling/momentum_exchange_method/SegreSilberbergMEM.cpp
index 38295ecaae6bcf82a43b0f686ba91aecceb8440d..7bd4e73a2e892b2bf4aceb6b0caa3107adf909c4 100644
--- a/tests/pe_coupling/momentum_exchange_method/SegreSilberbergMEM.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/SegreSilberbergMEM.cpp
@@ -603,7 +603,7 @@ int main( int argc, char **argv )
 
    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");
+   auto fcdID          = blocks->addBlockData(pe::fcd::createGenericFCDDataHandling<BodyTypeTuple, pe::fcd::AnalyticCollideFunctor>(), "FCD");
 
    // set up collision response, here DEM solver
    pe::cr::DEM cr(globalBodyStorage, blocks->getBlockStoragePointer(), bodyStorageID, ccdID, fcdID, NULL);
diff --git a/tests/pe_coupling/momentum_exchange_method/SquirmerTest.cpp b/tests/pe_coupling/momentum_exchange_method/SquirmerTest.cpp
index 235fcc2a62739f4d5cb559e164a34cf143db82a3..f0c248c1419f2ebaf1a4f2ae33736273dbae9509 100644
--- a/tests/pe_coupling/momentum_exchange_method/SquirmerTest.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/SquirmerTest.cpp
@@ -270,7 +270,7 @@ int main( int argc, char **argv )
    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");
+   auto fcdID         = blocks->addBlockData(pe::fcd::createGenericFCDDataHandling<BodyTypeTuple, pe::fcd::AnalyticCollideFunctor>(), "FCD");
    
    pe::cr::HCSITS cr(globalBodyStorage, blocks->getBlockStoragePointer(), bodyStorageID, ccdID, fcdID);
    
diff --git a/tests/pe_coupling/partially_saturated_cells_method/SegreSilberbergPSM.cpp b/tests/pe_coupling/partially_saturated_cells_method/SegreSilberbergPSM.cpp
index 4d1cef74abc117df8fc3a16c6db23c5ddd28d2d7..ef3378bcace821843d098592ee5781e0739fb141 100644
--- a/tests/pe_coupling/partially_saturated_cells_method/SegreSilberbergPSM.cpp
+++ b/tests/pe_coupling/partially_saturated_cells_method/SegreSilberbergPSM.cpp
@@ -560,7 +560,7 @@ int main( int argc, char **argv )
 
    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");
+   auto fcdID          = blocks->addBlockData(pe::fcd::createGenericFCDDataHandling<BodyTypeTuple, pe::fcd::AnalyticCollideFunctor>(), "FCD");
 
    // set up collision response, here DEM solver
    pe::cr::DEM cr( globalBodyStorage, blocks->getBlockStoragePointer(), bodyStorageID, ccdID, fcdID, NULL );