diff --git a/src/pe/Types.h b/src/pe/Types.h
index 35d559f4f7c19d81a0ba36d04482a9d7639c7c43..764bf4bb70610b43ec54215b60c85a95219ae25f 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 73f222c15c55051d33474390cc6df2512193f2f6..e2d2c4a6ceb2b4e35d1aad180c9b4700514bc03f 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 0000000000000000000000000000000000000000..d7e7199a4e4186b467734ac36b778f276eb223bb
--- /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 0000000000000000000000000000000000000000..57916d268a30dbfea187071a18b4270ee63c83f3
--- /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 0000000000000000000000000000000000000000..c3a0b6bb8c4c48b745b1f7db34511d96fb7af0eb
--- /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 0000000000000000000000000000000000000000..f2baad173f8f64029055721511af7c220a06e400
--- /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 0000000000000000000000000000000000000000..e848dd28527d9c763c55bab2b7e02467c82b05a3
--- /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 0000000000000000000000000000000000000000..78b335920d0eb92390fff93fdd35f6ea0e08c5f0
--- /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