From 33ba96209cb3f84f04649549a93a37f2a6ba5061 Mon Sep 17 00:00:00 2001
From: Sebastian Eibl <sebastian.eibl@fau.de>
Date: Fri, 10 Nov 2017 13:17:53 +0100
Subject: [PATCH] added ellipsoid class

---
 src/pe/Types.h                               |   5 +
 src/pe/communication/DynamicMarshalling.h    |   1 +
 src/pe/communication/rigidbody/Ellipsoid.cpp |  59 ++++
 src/pe/communication/rigidbody/Ellipsoid.h   |  78 +++++
 src/pe/rigidbody/Ellipsoid.cpp               | 264 ++++++++++++++
 src/pe/rigidbody/Ellipsoid.h                 | 348 +++++++++++++++++++
 src/pe/rigidbody/EllipsoidFactory.cpp        |  83 +++++
 src/pe/rigidbody/EllipsoidFactory.h          |  69 ++++
 8 files changed, 907 insertions(+)
 create mode 100644 src/pe/communication/rigidbody/Ellipsoid.cpp
 create mode 100644 src/pe/communication/rigidbody/Ellipsoid.h
 create mode 100644 src/pe/rigidbody/Ellipsoid.cpp
 create mode 100644 src/pe/rigidbody/Ellipsoid.h
 create mode 100644 src/pe/rigidbody/EllipsoidFactory.cpp
 create mode 100644 src/pe/rigidbody/EllipsoidFactory.h

diff --git a/src/pe/Types.h b/src/pe/Types.h
index 35d559f4f..764bf4bb7 100644
--- a/src/pe/Types.h
+++ b/src/pe/Types.h
@@ -65,6 +65,7 @@ class Box;
 class Capsule;
 class Contact;
 class Cylinder;
+class Ellipsoid;
 class CylindricalBoundary;
 class FixedJoint;
 class ForceGenerator;
@@ -111,6 +112,10 @@ typedef CylindricalBoundary        CylindricalBoundaryType;        //!< Type of
 typedef CylindricalBoundary*       CylindricalBoundaryID;          //!< Handle for a cylindrical boundary primitive.
 typedef const CylindricalBoundary* ConstCylindricalBoundaryID;     //!< Handle for a constant cylindrical boundary primitive.
 
+typedef Ellipsoid             EllipsoidType;       //!< Type of the ellipsoid geometric primitive.
+typedef Ellipsoid*            EllipsoidID;         //!< Handle for a ellipsoid primitive.
+typedef const Ellipsoid*      ConstEllipsoidID;    //!< Handle for a constant ellipsoid primitive.
+
 typedef Plane                 PlaneType;           //!< Type of the plane geometric primitive.
 typedef Plane*                PlaneID;             //!< Handle for a plane primitive.
 typedef const Plane*          ConstPlaneID;        //!< Handle for a constant plane primitive.
diff --git a/src/pe/communication/DynamicMarshalling.h b/src/pe/communication/DynamicMarshalling.h
index 73f222c15..e2d2c4a6c 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/communication/rigidbody/Ellipsoid.h"
 #include "pe/utility/BodyCast.h"
 
 #include "core/Abort.h"
diff --git a/src/pe/communication/rigidbody/Ellipsoid.cpp b/src/pe/communication/rigidbody/Ellipsoid.cpp
new file mode 100644
index 000000000..d7e7199a4
--- /dev/null
+++ b/src/pe/communication/rigidbody/Ellipsoid.cpp
@@ -0,0 +1,59 @@
+//======================================================================================================================
+//
+//  This file is part of waLBerla. waLBerla is free software: you can
+//  redistribute it and/or modify it under the terms of the GNU General Public
+//  License as published by the Free Software Foundation, either version 3 of
+//  the License, or (at your option) any later version.
+//
+//  waLBerla is distributed in the hope that it will be useful, but WITHOUT
+//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+//  for more details.
+//
+//  You should have received a copy of the GNU General Public License along
+//  with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file Ellipsoid.cpp
+//! \author Tobias Preclik
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//! \brief Marshalling of objects for data transmission or storage.
+//
+//======================================================================================================================
+
+#include "Ellipsoid.h"
+
+namespace walberla {
+namespace pe {
+namespace communication {
+
+//*************************************************************************************************
+/*!\brief Marshalling a Ellipsoid primitive.
+ *
+ * \param buffer The buffer to be filled.
+ * \param obj The object to be marshalled.
+ * \return void
+ */
+void marshal( mpi::SendBuffer& buffer, const Ellipsoid& obj ) {
+   marshal( buffer, static_cast<const GeomPrimitive&>( obj ) );
+   buffer << obj.getSemiAxes();
+}
+//*************************************************************************************************
+
+
+//*************************************************************************************************
+/*!\brief Unmarshalling a Ellipsoid primitive.
+ *
+ * \param buffer The buffer from where to read.
+ * \param objparam The object to be reconstructed.
+ * \param hasSuperBody False if body is not part of a union. Passed on to rigid body unmarshalling.
+ * \return void
+ */
+void unmarshal( mpi::RecvBuffer& buffer, EllipsoidParameters& objparam ) {
+   unmarshal( buffer, static_cast<GeomPrimitiveParameters&>( objparam ) );
+   buffer >> objparam.semiAxes_;
+}
+//*************************************************************************************************
+
+}  // namespace communication
+}  // namespace pe
+}  // namespace walberla
diff --git a/src/pe/communication/rigidbody/Ellipsoid.h b/src/pe/communication/rigidbody/Ellipsoid.h
new file mode 100644
index 000000000..57916d268
--- /dev/null
+++ b/src/pe/communication/rigidbody/Ellipsoid.h
@@ -0,0 +1,78 @@
+//======================================================================================================================
+//
+//  This file is part of waLBerla. waLBerla is free software: you can
+//  redistribute it and/or modify it under the terms of the GNU General Public
+//  License as published by the Free Software Foundation, either version 3 of
+//  the License, or (at your option) any later version.
+//
+//  waLBerla is distributed in the hope that it will be useful, but WITHOUT
+//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+//  for more details.
+//
+//  You should have received a copy of the GNU General Public License along
+//  with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file Ellipsoid.h
+//! \author Tobias Preclik
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//! \brief Marshalling of objects for data transmission or storage.
+//
+//======================================================================================================================
+
+#pragma once
+
+
+//*************************************************************************************************
+// Includes
+//*************************************************************************************************
+
+#include "pe/communication/Instantiate.h"
+#include "pe/communication/Marshalling.h"
+#include "pe/rigidbody/Ellipsoid.h"
+
+namespace walberla {
+namespace pe {
+namespace communication {
+
+struct EllipsoidParameters : public GeomPrimitiveParameters {
+   Vec3 semiAxes_;
+};
+
+//*************************************************************************************************
+/*!\brief Marshalling a Ellipsoid primitive.
+ *
+ * \param buffer The buffer to be filled.
+ * \param obj The object to be marshalled.
+ * \return void
+ */
+void marshal( mpi::SendBuffer& buffer, const Ellipsoid& obj );
+//*************************************************************************************************
+
+//*************************************************************************************************
+/*!\brief Unmarshalling a Ellipsoid primitive.
+ *
+ * \param buffer The buffer from where to read.
+ * \param objparam The object to be reconstructed.
+ * \param hasSuperBody False if body is not part of a union. Passed on to rigid body unmarshalling.
+ * \return void
+ */
+void unmarshal( mpi::RecvBuffer& buffer, EllipsoidParameters& objparam );
+//*************************************************************************************************
+
+
+inline EllipsoidID instantiate( mpi::RecvBuffer& buffer, const math::AABB& domain, const math::AABB& block, EllipsoidID& newBody )
+{
+   EllipsoidParameters subobjparam;
+   unmarshal( buffer, subobjparam );
+   correctBodyPosition(domain, block.center(), subobjparam.gpos_);
+   newBody = new Ellipsoid( subobjparam.sid_, subobjparam.uid_, subobjparam.gpos_, subobjparam.rpos_, subobjparam.q_, subobjparam.semiAxes_, subobjparam.material_, false, subobjparam.communicating_, subobjparam.infiniteMass_ );
+   newBody->setLinearVel( subobjparam.v_ );
+   newBody->setAngularVel( subobjparam.w_ );
+   newBody->MPITrait.setOwner( subobjparam.mpiTrait_.owner_ );
+   return newBody;
+}
+
+}  // namespace communication
+}  // namespace pe
+}  // namespace walberla
diff --git a/src/pe/rigidbody/Ellipsoid.cpp b/src/pe/rigidbody/Ellipsoid.cpp
new file mode 100644
index 000000000..c3a0b6bb8
--- /dev/null
+++ b/src/pe/rigidbody/Ellipsoid.cpp
@@ -0,0 +1,264 @@
+//======================================================================================================================
+//
+//  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 Ellipsoid.cpp
+//! \author Klaus Iglberger
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#include "Ellipsoid.h"
+
+//*************************************************************************************************
+// Includes
+//*************************************************************************************************
+
+#include <iomanip>
+#include <iostream>
+#include <stdexcept>
+#include <cmath>
+#include <pe/Materials.h>
+#include <core/math/Matrix3.h>
+#include <core/debug/Debug.h>
+
+namespace walberla {
+namespace pe {
+
+//=================================================================================================
+//
+//  CONSTRUCTOR
+//
+//=================================================================================================
+
+
+//*************************************************************************************************
+//*************************************************************************************************
+/*!\brief Constructor for the Ellipsoid class.
+ *
+ * \param sid Unique system-specific ID for the Ellipsoid.
+ * \param uid User-specific ID for the Ellipsoid.
+ * \param gpos Global geometric center of the Ellipsoid.
+ * \param rpos The relative position within the body frame of a superordinate body.
+ * \param q The orientation of the Ellipsoid's body frame in the global world frame.
+ * \param radius The radius of the Ellipsoid \f$ (0..\infty) \f$.
+ * \param material The material of the Ellipsoid.
+ * \param global specifies if the Ellipsoid should be created in the global storage
+ * \param communicating specifies if the Ellipsoid should take part in synchronization (syncNextNeighbour, syncShadowOwner)
+ * \param infiniteMass specifies if the Ellipsoid has infinite mass and will be treated as an obstacle
+ */
+Ellipsoid::Ellipsoid( id_t sid, id_t uid, const Vec3& gpos, const Vec3& rpos, const Quat& q,
+                const Vec3& semiAxes, MaterialID material,
+                const bool global, const bool communicating, const bool infiniteMass )
+   : Ellipsoid::Ellipsoid( getStaticTypeID(), sid, uid, gpos, rpos, q, semiAxes, material, global, communicating, infiniteMass )
+{}
+Ellipsoid::Ellipsoid( id_t const typeId, id_t sid, id_t uid, const Vec3& gpos, const Vec3& rpos, const Quat& q,
+                const Vec3& semiAxes, MaterialID material,
+                const bool global, const bool communicating, const bool infiniteMass )
+   : GeomPrimitive( typeId, sid, uid, material )  // Initialization of the parent class
+   , semiAxes_(semiAxes)
+{
+   // Checking the radius
+   // Since the Ellipsoid constructor is never directly called but only used in a small number
+   // of functions that already check the Ellipsoid arguments, only asserts are used here to
+   // double check the arguments.
+   WALBERLA_ASSERT( semiAxes_[0] > real_c(0), "Invalid Ellipsoid radius" );
+   WALBERLA_ASSERT( semiAxes_[1] > real_c(0), "Invalid Ellipsoid radius" );
+   WALBERLA_ASSERT( semiAxes_[2] > real_c(0), "Invalid Ellipsoid radius" );
+   // Setting the center of the Ellipsoid
+   gpos_ = gpos;
+
+   // Initializing the instantiated Ellipsoid
+   rpos_   = rpos;                   // Setting the relative position
+   q_      = q;                      // Setting the orientation
+   R_      = q_.toRotationMatrix();  // Setting the rotation matrix
+
+   // Calculating the Ellipsoid mass
+   mass_ = calcMass(semiAxes_, Material::getDensity( material ));
+   invMass_ = real_c(1) / mass_;
+
+   // Calculating the moment of inertia
+   calcInertia();
+
+   setGlobal( global );
+   setMass( infiniteMass );
+   setCommunicating( communicating );
+   setFinite( true );
+
+   // Setting the axis-aligned bounding box
+   Ellipsoid::calcBoundingBox();
+}
+//*************************************************************************************************
+
+
+
+
+//=================================================================================================
+//
+//  DESTRUCTOR
+//
+//=================================================================================================
+
+//*************************************************************************************************
+/*!\brief Destructor for the Ellipsoid class.
+ */
+Ellipsoid::~Ellipsoid()
+{
+   // Logging the destruction of the Ellipsoid
+   WALBERLA_LOG_DETAIL( "Destroyed Ellipsoid " << sid_ );
+}
+//*************************************************************************************************
+
+
+
+
+//=================================================================================================
+//
+//  UTILITY FUNCTIONS
+//
+//=================================================================================================
+
+//*************************************************************************************************
+/*!\brief Checks, whether a point in body relative coordinates lies inside the Ellipsoid.
+ *
+ * \param px The x-component of the relative coordinate.
+ * \param py The y-component of the relative coordinate.
+ * \param pz The z-component of the relative coordinate.
+ * \return \a true if the point lies inside the Ellipsoid, \a false if not.
+ */
+bool Ellipsoid::containsRelPointImpl( real_t px, real_t py, real_t pz ) const
+{
+return ( (px * px)/(semiAxes_[0] * semiAxes_[0]) + (py * py)/(semiAxes_[1] * semiAxes_[1]) 
+		+ (pz * pz)/(semiAxes_[2] * semiAxes_[2]) <= real_t(1.0) );
+}
+//*************************************************************************************************
+
+
+//*************************************************************************************************
+/*!\brief Checks, whether a point in body relative coordinates lies on the surface of the Ellipsoid.
+ *
+ * \param px The x-component of the relative coordinate.
+ * \param py The y-component of the relative coordinate.
+ * \param pz The z-component of the relative coordinate.
+ * \return \a true if the point lies on the surface of the Ellipsoid, \a false if not.
+ *
+ * The (relative) tolerance level of the check is pe::surfaceThreshold.
+ */
+bool Ellipsoid::isSurfaceRelPointImpl( real_t px, real_t py, real_t pz ) const
+{
+   return floatIsEqual( (px * px)/(semiAxes_[0] * semiAxes_[0]) + (py * py)/(semiAxes_[1] * semiAxes_[1]) 
+		+ (pz * pz)/(semiAxes_[2] * semiAxes_[2]), real_t(1.0), pe::surfaceThreshold);
+}
+//*************************************************************************************************
+
+
+
+
+//=================================================================================================
+//
+//  OUTPUT FUNCTIONS
+//
+//=================================================================================================
+
+//*************************************************************************************************
+/*!\brief Output of the current state of a Ellipsoid.
+ *
+ * \param os Reference to the output stream.
+ * \param tab Indentation in front of every line of the Ellipsoid output.
+ * \return void
+ */
+void Ellipsoid::print( std::ostream& os, const char* tab ) const
+{
+   using std::setw;
+
+   os << tab << " Ellipsoid " << uid_ << " with semi-axis " << semiAxes_ << "\n";
+
+   os << tab << "   Fixed: " << isFixed() << " , sleeping: " << !isAwake() << "\n";
+
+   os << tab << "   System ID         = " << getSystemID() << "\n"
+      << tab << "   Total mass        = " << getMass() << "\n"
+      << tab << "   Material          = " << Material::getName( material_ ) << "\n"
+      << tab << "   Owner             = " << MPITrait.getOwner() << "\n"
+      << tab << "   Global position   = " << getPosition() << "\n"
+      << tab << "   Relative position = " << getRelPosition() << "\n"
+      << tab << "   Linear velocity   = " << getLinearVel() << "\n"
+      << tab << "   Angular velocity  = " << getAngularVel() << "\n";
+
+//   if( verboseMode )
+//   {
+      os << tab << "   Bounding box      = " << getAABB() << "\n"
+         << tab << "   Quaternion        = " << getQuaternion() << "\n"
+         << tab << "   Rotation matrix   = ( " << setw(9) << R_[0] << " , " << setw(9) << R_[1] << " , " << setw(9) << R_[2] << " )\n"
+         << tab << "                       ( " << setw(9) << R_[3] << " , " << setw(9) << R_[4] << " , " << setw(9) << R_[5] << " )\n"
+         << tab << "                       ( " << setw(9) << R_[6] << " , " << setw(9) << R_[7] << " , " << setw(9) << R_[8] << " )\n";
+
+      os << std::setiosflags(std::ios::right)
+         << tab << "   Moment of inertia = ( " << setw(9) << I_[0] << " , " << setw(9) << I_[1] << " , " << setw(9) << I_[2] << " )\n"
+         << tab << "                       ( " << setw(9) << I_[3] << " , " << setw(9) << I_[4] << " , " << setw(9) << I_[5] << " )\n"
+         << tab << "                       ( " << setw(9) << I_[6] << " , " << setw(9) << I_[7] << " , " << setw(9) << I_[8] << " )\n"
+         << std::resetiosflags(std::ios::right);
+//   }
+}
+//*************************************************************************************************
+
+
+
+
+//=================================================================================================
+//
+//  GLOBAL OPERATORS
+//
+//=================================================================================================
+
+//*************************************************************************************************
+/*!\brief Global output operator for Ellipsoids.
+ *
+ * \param os Reference to the output stream.
+ * \param s Reference to a constant Ellipsoid object.
+ * \return Reference to the output stream.
+ */
+std::ostream& operator<<( std::ostream& os, const Ellipsoid& s )
+{
+   os << "--" << "Ellipsoid PARAMETERS"
+      << "-------------------------------------------------------------\n";
+   s.print( os, "" );
+   os << "--------------------------------------------------------------------------------\n"
+      << std::endl;
+   return os;
+}
+//*************************************************************************************************
+
+
+//*************************************************************************************************
+/*!\brief Global output operator for Ellipsoid handles.
+ *
+ * \param os Reference to the output stream.
+ * \param s Constant Ellipsoid handle.
+ * \return Reference to the output stream.
+ */
+std::ostream& operator<<( std::ostream& os, ConstEllipsoidID s )
+{
+   os << "--" << "Ellipsoid PARAMETERS"
+      << "-------------------------------------------------------------\n";
+   s->print( os, "" );
+   os << "--------------------------------------------------------------------------------\n"
+      << std::endl;
+   return os;
+}
+//*************************************************************************************************
+
+id_t Ellipsoid::staticTypeID_ = std::numeric_limits<id_t>::max();
+
+} // namespace pe
+}  // namespace walberla
diff --git a/src/pe/rigidbody/Ellipsoid.h b/src/pe/rigidbody/Ellipsoid.h
new file mode 100644
index 000000000..f2baad173
--- /dev/null
+++ b/src/pe/rigidbody/Ellipsoid.h
@@ -0,0 +1,348 @@
+//======================================================================================================================
+//
+//  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 Ellipsoid.h
+//! \author Klaus Iglberger
+//! \author Sebastian Eibl <sebastian.eibl@fau.de
+//! \author Tobias Leemann <tobias.leemann@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+
+//*************************************************************************************************
+// Includes
+//*************************************************************************************************
+
+#include <pe/rigidbody/GeomPrimitive.h>
+#include <pe/Types.h>
+#include <core/math/Constants.h>
+#include <core/math/Matrix3.h>
+#include <core/math/Vector3.h>
+#include <core/DataTypes.h>
+#include <core/logging/Logging.h>
+#include <core/math/Constants.h>
+#include <core/math/Limits.h>
+#include <core/math/Utility.h>
+#include <pe/Config.h>
+
+
+namespace walberla {
+namespace pe {
+
+//=================================================================================================
+//
+//  CLASS DEFINITION
+//
+//=================================================================================================
+
+//*************************************************************************************************
+/**
+ * \ingroup pe
+ * \brief Base class for the Ellipsoid geometry.
+ *
+ * The Ellipsoid class represents the base class for the Ellipsoid geometry. It provides
+ * the basic functionality of a Ellipsoid. An Ellipsoid is obtained from a sphere by deforming it by means 
+ * of directional scalings. Its three semi-axes corresponding to the x, y, z direction are labeled 
+ * a, b, c.
+ */
+class Ellipsoid : public GeomPrimitive
+{
+public:
+   //**Constructors********************************************************************************
+   /*!\name Constructors */
+   //@{
+   explicit Ellipsoid( id_t sid, id_t uid, const Vec3& gpos, const Vec3& rpos, const Quat& q,
+                    const Vec3& semiAxes, MaterialID material,
+                    const bool global, const bool communicating, const bool infiniteMass );
+   explicit Ellipsoid( id_t const typeID, id_t sid, id_t uid, const Vec3& gpos, const Vec3& rpos, const Quat& q,
+                    const Vec3& semiAxes, MaterialID material,
+                    const bool global, const bool communicating, const bool infiniteMass );
+   //@}
+   //**********************************************************************************************
+
+   //**Destructor**********************************************************************************
+   /*!\name Destructor */
+   //@{
+   virtual ~Ellipsoid();
+   //@}
+   //**********************************************************************************************
+   //**********************************************************************************************
+
+public:
+   //**Get functions*******************************************************************************
+   /*!\name Get functions */
+   //@{
+   inline const Vec3& getSemiAxes() const;
+   virtual inline real_t getVolume()         const;
+   //@}
+   //**********************************************************************************************
+
+   //**Utility functions***************************************************************************
+   static inline id_t getStaticTypeID();
+   //@}
+   //**********************************************************************************************
+
+   //**Output functions****************************************************************************
+   /*!\name Output functions */
+   //@{
+   virtual void print( std::ostream& os, const char* tab ) const;
+   //@}
+   //**********************************************************************************************
+
+   //**Utility functions***************************************************************************
+   /*!\name Utility functions */
+   //@{
+   inline virtual Vec3 support( const Vec3& d ) const;
+   //@}
+   //**********************************************************************************************
+
+   //**Volume, mass and density functions**********************************************************
+   /*!\name Volume, mass and density functions */
+   //@{
+   static inline real_t calcVolume( const Vec3& semiAxes );
+   static inline real_t calcMass( const Vec3& semiAxes, real_t density );
+   static inline real_t calcDensity( const Vec3& semiAxes, real_t mass );
+   //@}
+   //**********************************************************************************************
+
+protected:
+   //**Utility functions***************************************************************************
+   /*!\name Utility functions */
+   //@{
+   virtual bool containsRelPointImpl ( real_t px, real_t py, real_t pz ) const;
+   virtual bool isSurfaceRelPointImpl( real_t px, real_t py, real_t pz ) const;
+   //@}
+   //**********************************************************************************************
+
+   //**Utility functions***************************************************************************
+   /*!\name Utility functions */
+   //@{
+   inline virtual void calcBoundingBox();  // Calculation of the axis-aligned bounding box
+   inline         void calcInertia();      // Calculation of the moment of inertia
+   //@}
+   //**********************************************************************************************
+
+   //**Member variables****************************************************************************
+   /*!\name Member variables */
+   //@{
+   Vec3 semiAxes_;  //!< Semi-axes of the Ellipsoid.
+                  /*!< The radius is constrained to values larger than 0.0. */
+   //@}
+   //**********************************************************************************************
+private:
+   static id_t staticTypeID_;  //< type id of Ellipsoid, will be set by SetBodyTypeIDs
+   static void setStaticTypeID(id_t typeID) {staticTypeID_ = typeID;}
+
+   //** friend declaration
+   /// needed to be able to set static type ids with setStaticTypeID
+   template <class T>
+   friend struct SetBodyTypeIDs;
+};
+//*************************************************************************************************
+
+
+
+
+//=================================================================================================
+//
+//  GET FUNCTIONS
+//
+//=================================================================================================
+
+//*************************************************************************************************
+/*!\brief Returns the semi-axes of the Ellipsoid.
+ *
+ * \return The semi-axes of the Ellipsoid.
+ */
+inline const Vec3& Ellipsoid::getSemiAxes() const
+{
+   return semiAxes_;
+}
+//*************************************************************************************************
+
+
+//*************************************************************************************************
+/*!\brief Returns the volume of the Ellipsoid.
+ *
+ * \return The volume of the Ellipsoid.
+ */
+inline real_t Ellipsoid::getVolume() const
+{
+   return Ellipsoid::calcVolume(getSemiAxes());
+}
+//*************************************************************************************************
+
+
+
+
+//=================================================================================================
+//
+//  VOLUME, MASS AND DENSITY FUNCTIONS
+//
+//=================================================================================================
+
+//*************************************************************************************************
+/*!\brief Calculates the volume of a Ellipsoid for a given vector of semi-axes.
+ *
+ * \param semiAxes The vector of semi-axes of the Ellipsoid.
+ * \return The volume of the Ellipsoid.
+ */
+inline real_t Ellipsoid::calcVolume(const Vec3& semiAxes ) 
+{
+   return real_c(4.0)/real_c(3.0) * math::M_PI * semiAxes[0] * semiAxes[1] * semiAxes[2];
+}
+//*************************************************************************************************
+
+
+//*************************************************************************************************
+/*!\brief Calculates the mass of a Ellipsoid for a given a given vector of semi-axes and density.
+ *
+ * \param semiAxes The vector of semi-axes of the Ellipsoid.
+ * \param density The density of the Ellipsoid.
+ * \return The total mass of the Ellipsoid.
+ */
+inline real_t Ellipsoid::calcMass(const Vec3& semiAxes, real_t density )
+{
+   return real_c(4.0)/real_c(3.0) * math::M_PI * semiAxes[0] * semiAxes[1] * semiAxes[2] * density;
+}
+//*************************************************************************************************
+
+
+//*************************************************************************************************
+/*!\brief Calculates the density of a Ellipsoid for a given vector of semi-axes and mass.
+ *
+ * \param semiAxes The vector of semi-axes of the Ellipsoid.
+ * \param mass The total mass of the Ellipsoid.
+ * \return The density of the Ellipsoid.
+ */
+inline real_t Ellipsoid::calcDensity(const Vec3& semiAxes, real_t mass )
+{
+   return real_c(0.75) * mass / ( math::M_PI * semiAxes[0] * semiAxes[1] * semiAxes[2] );
+}
+//*************************************************************************************************
+
+
+
+
+//=================================================================================================
+//
+//  UTILITY FUNCTIONS
+//
+//=================================================================================================
+
+//*************************************************************************************************
+/*!\brief Calculation of the bounding box of the Ellipsoid.
+ *
+ * \return void
+ *
+ * This function updates the axis-aligned bounding box of the Ellipsoid primitive according to the
+ * current position and orientation of the Ellipsoid. Note that the bounding box is increased in
+ * all dimensions by pe::contactThreshold to guarantee that rigid bodies in close proximity of
+ * the Ellipsoid are also considered during the collision detection process.
+ * Algorithm: Use a non-axes-aligned bounding box of the ellipse (box
+ * with sides twice the semi-axes long) and calc its AABB.
+ */
+inline void Ellipsoid::calcBoundingBox()
+{	
+   using std::fabs;
+
+   const real_t xlength( fabs(R_[0]*semiAxes_[0]) + fabs(R_[1]*semiAxes_[1]) + fabs(R_[2]*semiAxes_[2])  + contactThreshold );
+   const real_t ylength( fabs(R_[3]*semiAxes_[0]) + fabs(R_[4]*semiAxes_[1]) + fabs(R_[5]*semiAxes_[2])  + contactThreshold );
+   const real_t zlength( fabs(R_[6]*semiAxes_[0]) + fabs(R_[7]*semiAxes_[1]) + fabs(R_[8]*semiAxes_[2])  + contactThreshold );
+   aabb_ = math::AABB(
+         gpos_[0] - xlength,
+         gpos_[1] - ylength,
+         gpos_[2] - zlength,
+         gpos_[0] + xlength,
+         gpos_[1] + ylength,
+         gpos_[2] + zlength
+         );
+	//   WALBERLA_ASSERT( aabb_.isValid()        , "Invalid bounding box detected" );
+   WALBERLA_ASSERT( aabb_.contains( gpos_ ), "Invalid bounding box detected("<< getSystemID() <<")\n" << "pos: " << gpos_ << "\nlengths: " << xlength << "," << ylength << ", " << zlength<< "\nvel: " << getLinearVel() << "\nbox: " << aabb_ );
+}
+//*************************************************************************************************
+
+
+//*************************************************************************************************
+/*!\brief Calculation of the moment of inertia in reference to the body frame of the Ellipsoid.
+ *
+ * \return void
+ */
+inline void Ellipsoid::calcInertia()
+{
+   I_[0] = real_c(0.2) * mass_ * (semiAxes_[1] * semiAxes_[1] + semiAxes_[2] * semiAxes_[2]);
+   I_[4] = real_c(0.2) * mass_ * (semiAxes_[2] * semiAxes_[2] + semiAxes_[0] * semiAxes_[0]);
+   I_[8] = real_c(0.2) * mass_ * (semiAxes_[0] * semiAxes_[0] + semiAxes_[1] * semiAxes_[1]);
+   Iinv_ = I_.getInverse();
+}
+//*************************************************************************************************
+
+
+//*************************************************************************************************
+/*!\brief Estimates the point which is farthest in direction \a d.
+ *
+ * \param d The normalized search direction in world-frame coordinates.
+ * \return The support point in world-frame coordinates in direction a\ d.
+ */
+inline Vec3 Ellipsoid::support( const Vec3& d ) const
+{
+   auto len = d.sqrLength();
+   if (!math::equal(len, real_t(0)))
+   {
+	  Vec3 d_loc = vectorFromWFtoBF(d);
+	  Vec3 norm_vec(d_loc[0] * semiAxes_[0], d_loc[1] * semiAxes_[1], d_loc[2] * semiAxes_[2]);
+	  real_t norm_length = norm_vec.length();
+	  Vec3 local_support = (real_t(1.0) / norm_length) * Vec3(semiAxes_[0] * semiAxes_[0] * d_loc[0], 
+			semiAxes_[1] * semiAxes_[1] * d_loc[1], semiAxes_[2] * semiAxes_[2] * d_loc[2]);
+      return pointFromBFtoWF(local_support);
+   } else
+   {
+      return Vec3(0,0,0);
+   }
+}
+//*************************************************************************************************
+
+
+//*************************************************************************************************
+/*!\brief Returns unique type id of this type
+ *
+ * \return geometry specific type id
+ */
+inline id_t Ellipsoid::getStaticTypeID()
+{
+   return staticTypeID_;
+}
+//*************************************************************************************************
+
+
+//=================================================================================================
+//
+//  GLOBAL OPERATORS
+//
+//=================================================================================================
+
+//*************************************************************************************************
+/*!\name Ellipsoid operators */
+//@{
+std::ostream& operator<<( std::ostream& os, const Ellipsoid& s );
+std::ostream& operator<<( std::ostream& os, ConstEllipsoidID s );
+//@}
+//*************************************************************************************************
+
+
+} // namespace pe
+}  // namespace walberla
diff --git a/src/pe/rigidbody/EllipsoidFactory.cpp b/src/pe/rigidbody/EllipsoidFactory.cpp
new file mode 100644
index 000000000..e848dd285
--- /dev/null
+++ b/src/pe/rigidbody/EllipsoidFactory.cpp
@@ -0,0 +1,83 @@
+//======================================================================================================================
+//
+//  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 EllipsoidFactory.cpp
+//! \author Tobias Leemann <tobias.leemann@fau.de>
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#include "EllipsoidFactory.h"
+
+#include "pe/Materials.h"
+#include "pe/rigidbody/BodyStorage.h"
+#include "pe/rigidbody/Ellipsoid.h"
+
+namespace walberla {
+namespace pe {
+
+EllipsoidID createEllipsoid( BodyStorage& globalStorage, BlockStorage& blocks, BlockDataID storageID,
+                       id_t uid, const Vec3& gpos, const Vec3& semiAxes,
+                       MaterialID material,
+                       bool global, bool communicating, bool infiniteMass )
+{
+   WALBERLA_ASSERT_UNEQUAL( Ellipsoid::getStaticTypeID(), std::numeric_limits<id_t>::max(), "Ellipsoid TypeID not initalized!");
+   // Checking the semiAxes
+   if( semiAxes[0] <= real_c(0) || semiAxes[1] <= real_c(0) || semiAxes[2] <= real_c(0) )
+      throw std::invalid_argument( "Invalid Ellipsoid semi-axes" );
+
+   EllipsoidID ellipsoid = NULL;
+
+   if (global)
+   {
+      const id_t sid = UniqueID<RigidBody>::createGlobal();
+      WALBERLA_ASSERT_EQUAL(communicating, false);
+      WALBERLA_ASSERT_EQUAL(infiniteMass, true);
+      ellipsoid = new Ellipsoid(sid, uid, gpos, Vec3(0,0,0), Quat(), semiAxes, material, global, false, true);
+      globalStorage.add(ellipsoid);
+   } else
+   {
+      for (auto it = blocks.begin(); it != blocks.end(); ++it){
+         IBlock* block = (&(*it));
+         if (block->getAABB().contains(gpos))
+         {
+            const id_t sid( UniqueID<RigidBody>::create() );
+
+            Storage* bs = block->getData<Storage>(storageID);
+            ellipsoid = new Ellipsoid(sid, uid, gpos, Vec3(0,0,0), Quat(), semiAxes, material, global, communicating, infiniteMass);
+            ellipsoid->MPITrait.setOwner(Owner(MPIManager::instance()->rank(), block->getId().getID()));
+            (*bs)[0].add(ellipsoid);
+         }
+      }
+   }
+
+   if (ellipsoid != NULL)
+   {
+      // Logging the successful creation of the Ellipsoid
+      WALBERLA_LOG_DETAIL(
+                "Created Ellipsoid " << ellipsoid->getSystemID() << "\n"
+             << "   User-ID         = " << uid << "\n"
+             << "   Global position = " << gpos << "\n"
+             << "   Semi-axes       = " << semiAxes << "\n"
+             << "   LinVel          = " << ellipsoid->getLinearVel() << "\n"
+             << "   Material        = " << Material::getName( material )
+               );
+   }
+
+   return ellipsoid;
+}
+
+}  // namespace pe
+}  // namespace walberla
diff --git a/src/pe/rigidbody/EllipsoidFactory.h b/src/pe/rigidbody/EllipsoidFactory.h
new file mode 100644
index 000000000..78b335920
--- /dev/null
+++ b/src/pe/rigidbody/EllipsoidFactory.h
@@ -0,0 +1,69 @@
+//======================================================================================================================
+//
+//  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 EllipsoidFactory.h
+//! \author Tobias Leemann
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "pe/Materials.h"
+#include "pe/Types.h"
+
+#include "blockforest/BlockForest.h"
+#include "core/debug/Debug.h"
+
+namespace walberla {
+namespace pe {
+
+//=================================================================================================
+//
+//  Ellipsoid SETUP FUNCTIONS
+//
+//=================================================================================================
+
+//*************************************************************************************************
+/**
+ * \ingroup pe
+ * \brief Setup of a new Ellipsoid.
+ *
+ * \param globalStorage process local global storage
+ * \param blocks storage of all the blocks on this process
+ * \param storageID BlockDataID of the BlockStorage block datum
+ * \param uid The user-specific ID of the Ellipsoid.
+ * \param gpos The global position of the center of the Ellipsoid.
+ * \param semiAxes The semiAxes of the Ellipsoid \f$ (0..\infty) \f$.
+ * \param material The material of the Ellipsoid.
+ * \param global specifies if the Ellipsoid should be created in the global storage
+ * \param communicating specifies if the Ellipsoid should take part in synchronization (syncNextNeighbour, syncShadowOwner)
+ * \param infiniteMass specifies if the Ellipsoid has infinite mass and will be treated as an obstacle
+ * \return Handle for the new Ellipsoid.
+ * \exception std::invalid_argument Invalid Ellipsoid semi-axes.
+ * \exception std::invalid_argument Invalid global Ellipsoid position.
+ *
+ * This function creates a Ellipsoid primitive in the \b pe simulation system. The Ellipsoid with
+ * user-specific ID \a uid is placed at the global position \a gpos, has the semi-axes \a semiAxes,
+ * and consists of the material \a material.
+ */
+EllipsoidID createEllipsoid( BodyStorage& globalStorage, BlockStorage& blocks, BlockDataID storageID,
+                       id_t uid, const Vec3& gpos, const Vec3& semiAxes,
+                       MaterialID material = Material::find("iron"),
+                       bool global = false, bool communicating = true, bool infiniteMass = false );
+//*************************************************************************************************
+
+}  // namespace pe
+}  // namespace walberla
-- 
GitLab