From 018b6cd5ef02173b7bbab88a3e302365da535cbc Mon Sep 17 00:00:00 2001
From: Sebastian Eibl <sebastian.eibl@fau.de>
Date: Wed, 7 Nov 2018 14:06:28 +0100
Subject: [PATCH] added Rot3 class and various free functions for vectors
---
src/core/math/Quaternion.h | 4 +-
src/core/math/Rot3.h | 124 +++++++++++++++++++++++++++++++++++++
src/core/math/Vector3.h | 57 +++++++++++++++++
src/core/mpi/RecvBuffer.h | 2 +-
4 files changed, 185 insertions(+), 2 deletions(-)
create mode 100644 src/core/math/Rot3.h
diff --git a/src/core/math/Quaternion.h b/src/core/math/Quaternion.h
index 4e531ea73..f6822daf3 100644
--- a/src/core/math/Quaternion.h
+++ b/src/core/math/Quaternion.h
@@ -36,6 +36,7 @@
#include "core/mpi/SendBuffer.h"
#include "core/mpi/RecvBuffer.h"
+#include <core/logging/Logging.h>
#include <boost/type_traits/is_floating_point.hpp>
#include <boost/type_traits/is_const.hpp>
@@ -406,7 +407,8 @@ inline Type Quaternion<Type>::operator[]( size_t index ) const
template< typename Type > // Data type of the quaternion
inline Quaternion<Type>& Quaternion<Type>::set( Type r, Type i, Type j, Type k )
{
- WALBERLA_CHECK_FLOAT_EQUAL( std::fabs( r*r + i*i + j*j + k*k ), Type(1), "Invalid quaternion parameters" );
+ WALBERLA_CHECK_FLOAT_EQUAL( std::fabs( r*r + i*i + j*j + k*k ), Type(1),
+ "Invalid quaternion parameters: " << r << ", "<< i << ", "<< j << ", "<< k );
v_[0] = r;
v_[1] = i;
v_[2] = j;
diff --git a/src/core/math/Rot3.h b/src/core/math/Rot3.h
new file mode 100644
index 000000000..30b4952d6
--- /dev/null
+++ b/src/core/math/Rot3.h
@@ -0,0 +1,124 @@
+//======================================================================================================================
+//
+// 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 Rot3.h
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include <core/math/Matrix3.h>
+#include <core/math/Quaternion.h>
+
+#include <type_traits>
+
+namespace walberla {
+namespace math {
+
+/**
+ * Rotation class which merges quaternion and matrix representation.
+ *
+ * For numerical reasons the representation of a rotation via a quaternion
+ * is favourable. However application of the rotation to vectors and matrices
+ * is numerical more efficient with matrices. Therefore this class combines both
+ * representations and takes care that both are in sync.
+ */
+template <typename Type> // floating point type
+class Rot3
+{
+ //**Compile time checks*************************************************************************
+ /*! \cond internal */
+ static_assert(std::is_floating_point<Type>::value, "T has to be floating point!");
+ static_assert(!std::is_const<Type>::value, "T has to be non const!");
+ static_assert(!std::is_volatile<Type>::value, "T has to be non volatile!");
+ /*! \endcond */
+ //**********************************************************************************************
+public:
+ Rot3(const Quaternion<Type>& q);
+
+ const Quaternion<Type>& getQuaternion() const { return quat_; }
+ const Matrix3<Type>& getMatrix() const { return mat_; }
+
+ void rotate(const Quaternion<Type>& q);
+private:
+ Quaternion<Type> quat_;
+ Matrix3<Type> mat_;
+};
+
+template< typename Type > // floating point type
+inline Rot3<Type>::Rot3(const Quaternion<Type>& q)
+ : quat_(q)
+ , mat_ (q.toRotationMatrix())
+{
+ WALBERLA_ASSERT_FLOAT_EQUAL( mat_.getDeterminant(), real_t(1), "Corrupted rotation matrix determinant" );
+}
+
+template< typename Type > // floating point type
+inline void Rot3<Type>::rotate(const Quaternion<Type>& q)
+{
+ quat_ = q * quat_;
+ mat_ = quat_.toRotationMatrix();
+ WALBERLA_ASSERT_FLOAT_EQUAL( mat_.getDeterminant(), real_t(1), "Corrupted rotation matrix determinant" );
+}
+
+template< typename Type > // floating point type
+inline std::ostream& operator<<( std::ostream& os, const Rot3<Type>& r )
+{
+ os << r.getMatrix();
+ return os;
+}
+
+} // math
+} // walberla
+
+//======================================================================================================================
+//
+// Send/Recv Buffer Serialization Specialization
+//
+//======================================================================================================================
+
+namespace walberla {
+namespace mpi {
+
+template< typename T, // Element type of SendBuffer
+ typename G, // Growth policy of SendBuffer
+ typename V > // value type
+mpi::GenericSendBuffer<T,G>& operator<<( mpi::GenericSendBuffer<T,G> & buf, const math::Rot3<V>& obj )
+{
+ buf.addDebugMarker( "ro" );
+ buf << obj.getQuaternion();
+ return buf;
+}
+
+template< typename T, // Element type of RecvBuffer
+ typename V > // value type
+mpi::GenericRecvBuffer<T>& operator>>( mpi::GenericRecvBuffer<T> & buf, math::Rot3<V>& objparam )
+{
+ math::Quaternion<V> q;
+ buf.readDebugMarker( "ro" );
+ buf >> q;
+ objparam = math::Rot3<V>(q);
+ return buf;
+}
+
+template< typename V > // value type
+struct BufferSizeTrait< math::Rot3<V> > {
+ static const bool constantSize = true;
+ static const uint_t size = BufferSizeTrait< math::Quaternion<V> >::size + mpi::BUFFER_DEBUG_OVERHEAD;
+};
+
+} // mpi
+} // walberla
diff --git a/src/core/math/Vector3.h b/src/core/math/Vector3.h
index da90e4cfd..24485fef5 100644
--- a/src/core/math/Vector3.h
+++ b/src/core/math/Vector3.h
@@ -593,6 +593,21 @@ inline Vector3<HIGH> Vector3<Type>::operator%( const Vector3<Other>& rhs ) const
//**********************************************************************************************************************
+//**********************************************************************************************************************
+/*!\brief Cross product (outer product) of two vectors (\f$ \vec{a}=\vec{b}\times\vec{c} \f$).
+//
+// \param lhs The left-hand-side vector for the cross product.
+// \param rhs The right-hand-side vector for the cross product.
+// \return The cross product.
+*/
+template< typename Type >
+inline Vector3<Type> cross( const Vector3<Type>& lhs, const Vector3<Type>& rhs )
+{
+ return lhs % rhs;
+}
+//**********************************************************************************************************************
+
+
//**********************************************************************************************************************
/*!\fn Vector3<HIGH> Vector3<Type>::operator+( const Vector3<Other>& rhs ) const
// \brief Addition operator for the addition of two vectors (\f$ \vec{a}=\vec{b}+\vec{c} \f$).
@@ -699,6 +714,13 @@ inline Vector3<HIGH> Vector3<Type>::operator/( Other rhs ) const
}
//**********************************************************************************************************************
+template< typename Type, typename Other >
+inline Vector3<HIGH> operator/( Other lhs, const Vector3<Type>& rhs )
+{
+ return Vector3<HIGH>( lhs/rhs[0], lhs/rhs[1], lhs/rhs[2] );
+}
+//**********************************************************************************************************************
+
//======================================================================================================================
//
@@ -1743,6 +1765,41 @@ Vector3<T> & normalize( Vector3<T> & v )
}
//**********************************************************************************************************************
+//**********************************************************************************************************************
+/**
+// \brief Length of the vector.
+//
+// \return Length of the vector.
+*/
+template<typename T>
+real_t length( const Vector3<T> & v )
+{
+ return v.length();
+}
+//**********************************************************************************************************************
+//**********************************************************************************************************************
+/**
+// \brief Length of the vector squared.
+//
+// \return Length of the vector squared.
+*/
+template<typename T>
+real_t sqrLength( const Vector3<T> & v )
+{
+ return v.sqrLength();
+}
+//**********************************************************************************************************************
+//**********************************************************************************************************************
+/**
+// \brief Dot product of two vectors.
+*/
+template<typename T>
+real_t dot( const Vector3<T> & v1, const Vector3<T> & v2 )
+{
+ return v1*v2;
+}
+//**********************************************************************************************************************
+
//**********************************************************************************************************************
/**
diff --git a/src/core/mpi/RecvBuffer.h b/src/core/mpi/RecvBuffer.h
index 3caa764b7..19d45825e 100644
--- a/src/core/mpi/RecvBuffer.h
+++ b/src/core/mpi/RecvBuffer.h
@@ -655,7 +655,7 @@ inline void GenericRecvBuffer<T>::readDebugMarker( const char * expectedMarker )
valid = false;
}
- WALBERLA_ASSERT( valid , "Packed and unpacked datatype do not match." );
+ WALBERLA_ASSERT( valid , "Packed and unpacked (" << expectedMarker << ") datatype do not match." );
}
#else
--
GitLab